基因差异表达分析——基于RSEM对比,DESeq2操作实例
以下操作全部基于R中的DESeq2函数包
使用DESeq2的两点要求:
- DEseq2要求输入数据是由整数组成的矩阵。
- DESeq2要求矩阵是没有标准化的。
下面是基本流程:
一、准备工作:
确定工作路径,以确保上传文件时无误
getwd() #显示当前工作目录
setwd() #设置工作目录
操作实例:
> getwd()
[1] "C:/Users/Larkin/Documents"
> setwd("C:/Users/Larkin/Desktop")
> getwd()
[1] "C:/Users/Larkin/Desktop"
二、数据上传及初步求整和整理
database <- read.table(file = "文件路径", sep = "\t", header = T, row.names = 1) #上传文件时,文件路径一定用“/”
database <- round(as.matrix(database)) #使用“round”函数,对数据求整;as.matrix()可将table转为matrix格式
其中,read.table() 函数中
参数sep = “\t"代表上传文件中将每列分隔的分隔符为”;”
参数header = T 代表上传文件中第一行的列数比其余行列数少一(header = F 代表所有行的列数都一样)
参数row.names = 1 代表将每一行以第一列中元素命名。
操作实例:
> database <- read.table(file = "C:/Users/Larkin/Desktop/48h.matrix", sep = "\t", header = T, row.names = 1)
得到下表:
继续运行
> database <- round(as.matrix(database))
则得到下表:
备注:
read.table()函数详细用法,请参考:
https://www.jianshu.com/p/90e1d430c9ef
三、构建dds矩阵
- dds构建前准备:
condition <- factor(c(rep("control",i),rep("exp",j))) #构建condition因子,作为基因差异表达的自变量(即实验组和对照组的对比)。i和j分别为control和exp的个数
coldata <- data.frame(row.names =colnames(database),condition) #将database中各个实验名称加上condition标签,即说明是control还是exp
library(DESeq2) #加载DESeq2包
其中,
rep(“ ”,i)代表将“”里元素重复i次。
data.frame( , ,) 代表将以逗号分隔开的几个元素放在一个dataframe里面
操作实例:
> condition <- factor(c(rep("control",3),rep("exp",3))) #3个对照,3个实验
> coldata <- data.frame(row.names = colnames(database),condition)
> coldatacondition
mus.48h.control1.results.genes.results control
mus.48h.control2.results.genes.results control
mus.48h.control3.results.genes.results control
mus.48h.exp1.results.genes.results exp
mus.48h.exp2.results.genes.results exp
mus.48h.exp3.results.genes.results exp
> library(DESeq2)
- 构建dds及利用results函数得到最终结果:
dds <- DESeqDataSetFromMatrix(countData = database,colData = coldata,design = ~condition) #countData用于说明数据来源,colData用于说明不同组数据的实验操作类型,design用于声明自变量,即谁和谁进行对比
dds <- DESeq(dds) #用于对dds数据进行运算及分析
res <- results(dds2,contrast = c("condition","control","exp")) #生成结果文件
具体操作请注意大小写
操作实例:
> dds <- DESeqDataSetFromMatrix(countData = database_1,colData = coldata,design = ~condition)
converting counts to integer mode> dds2 <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> res <- results(dds2,contrast = c("condition","control","exp"))
> res
log2 fold change (MLE): condition control vs exp
Wald test p-value: condition control vs exp
DataFrame with 42202 rows and 6 columnsbaseMean log2FoldChange lfcSE stat<numeric> <numeric> <numeric> <numeric>
gene:ENSMUSG00000000001 4975.47463350681 0.40129412229032 0.158622834629876 2.52986351698155
gene:ENSMUSG00000000003 0 NA NA NA
gene:ENSMUSG00000000028 1354.30376755978 0.0279337153287204 0.145918838510525 0.191433235172753
gene:ENSMUSG00000000037 0.167499533540566 0.850236344411763 4.0804728567969 0.208367111913393
gene:ENSMUSG00000000049 0 NA NA NA
... ... ... ... ...
gene:ENSMUSG00000118636 0.342472921967151 -0.111544335017095 4.0804728567969 -0.0273361296427433
gene:ENSMUSG00000118637 0 NA NA NA
gene:ENSMUSG00000118638 0 NA NA NA
gene:ENSMUSG00000118639 0 NA NA NA
gene:ENSMUSG00000118640 24.0588804554878 1.73883481104997 0.503013829968203 3.45683300826916pvalue padj<numeric> <numeric>
gene:ENSMUSG00000000001 0.0114106903449377 0.0366334294881918
gene:ENSMUSG00000000003 NA NA
gene:ENSMUSG00000000028 0.848186183621707 0.912698954692459
gene:ENSMUSG00000000037 0.834942333626489 NA
gene:ENSMUSG00000000049 NA NA
... ... ...
gene:ENSMUSG00000118636 0.978191640340057 NA
gene:ENSMUSG00000118637 NA NA
gene:ENSMUSG00000118638 NA NA
gene:ENSMUSG00000118639 NA NA
gene:ENSMUSG00000118640 0.000546563433122328 0.00292742174196498
至此,分析过程结束
四、结果文件整理及输出
res <- res[order(res$padj),] #将结果文件按照res文件中的padj这一列进行降序排列,其中$符号表示res中的padj这一列
resdata <- merge(as.data.frame(res),as.data.frame(counts(dds,normalized = TRUE)),by = "row.names",sort =FALSE) #合并结果文件res和构建的数据帧文件dds
write.csv(resdata,file = "结果文件名.csv") #在工作目录中生成结果文件
操作实例:
> res <- res[order(res$padj),]
> resdata <- merge(as.data.frame(res),as.data.frame(counts(dds,normalized = TRUE)),by = "row.names",sort =FALSE)
> write.csv(resdata,file = "48h_results.csv")
备注:
write.csv函数详细用法,请参考:
https://blog.csdn.net/u013421629/article/details/72771241
以上是做基因差异表达分析的基本流程
完成!
基因差异表达分析——基于RSEM对比,DESeq2操作实例相关推荐
- ballgown包进行基因差异表达分析
ballgown包可以读入Stringtie 的转录组组装及定量数据,进行基因差异表达分析. 1. 数据读入 # if (!requireNamespace("BiocManager&quo ...
- edger多组差异性分析_R语言利用edgeR package进行基因差异表达分析 举例
R语言利用edgeR package进行基因差异表达分析 举例 实验数据: 同一组织,分为两组,control vs treat,每组7例sample.数据第一列为基因名,后14列为对应的count. ...
- mysql修改工资字段_基于Linux的MySQL操作实例(修改表结构,MySQL索引,MySQL数据引擎)...
基于Linux的MySQL操作实例(修改表结构,MySQL索引,MySQL数据引擎) 前言 本篇是基于Linux下针对MySQL表结构的修改,MySQL索引的操作以及MySQL数据引擎的配置和说明. ...
- 哈佛大学——差异表达分析(九)DESeq2步骤描述
文章目录 学习目标 DESeq2差异基因表达分析流程 第一步:估计大小因子 第二步:估计基因离散(gene-wise dispersion) 第三步:拟合曲线到基因的分散估计 第四步:将基因离散估计值 ...
- 差异表达分析之FDR
差异表达分析之FDR 随着测序成本的不断降低,转录组测序分析已逐渐成为一种很常用的分析手段.但对于转录组分析当中的一些概念,很多人还不是很清楚.今天,小编就来谈谈在转录组分析中,经常会遇到的一个概念F ...
- edger多组差异性分析_简单使用DESeq2/EdgeR做差异分析
DESeq2和EdgeR都可用于做基因差异表达分析,主要也是用于RNA-Seq数据,同样也可以处理类似的ChIP-Seq,shRNA以及质谱数据. 这两个都属于R包,其相同点在于都是对count da ...
- edger多组差异性分析_简单使用DESeq2/EdgeR做差异分析 – 生信笔记
DESeq2和EdgeR都可用于做基因差异表达分析,主要也是用于RNA-Seq数据,同样也可以处理类似的ChIP-Seq,shRNA以及质谱数据. 这两个都属于R包,其相同点在于都是对count da ...
- 哈佛大学——差异表达分析(七)设计公式(Design formulas)
文章目录 学习目标 利用DESeq2进行差异表达分析 运行DESeq2 设计公式(design formula) 复杂的设计 MOV10 差异表达分析 学习目标 使用DESeq2执行差异表达分析工作流 ...
- 「用一个更复杂的例子,来深入学习DESeq2差异表达分析后的小分析」
这篇文章,对Griffith Lab的DESeq2分析流程做一个解读. 理解数据 Griffith Lab所使用的基因表达量矩阵总共包含了54个sample,这些sample可以划分为1)normal ...
最新文章
- CNCF 2019 年度报告重磅发布 | 云原生生态周报 Vol. 41
- 汇编指令prefix rep:
- zuul 1.x 和gateway性能对比
- umask命令:设置文件的默认权限掩码
- 2.5 指数加权平均的偏差修正
- ORACLE expdp/impdp导出实例
- Python 紧急修复远程代码执行漏洞
- Linux系统编程(18)——正则表达式实用举例
- linux进程自动启动,linux 嵌入式 自启动 系统自动登录-自动启动程序或脚本
- python中md5方法返回值_python中的md5加密
- 熟悉Linux基本操作
- 移动端人脸识别活体检测,高效集成
- Java流处理之转换编码的转换流
- NIOS ii 实战篇 --- 按键控制LED
- rm -rf /* 数据恢复记录
- 微信小程序的家校通系统(家校联系)
- 视频监控的2017 有什么看头?
- 厄米高斯模式与拉盖尔高斯模式
- c语言小学生四则运算出题程序
- Could not set property ‘xxx‘ of ‘class ‘xxx‘ with value ‘xxx‘,Cause: java.lang.IllegalArgumentExcept
热门文章
- 三菱FX3U与东元Teco变频器N310通讯实战程序
- 【人脸识别实战一】系统架构设计
- LeetCode题解(0197):查询相较于昨天温度上升的日期(SQL)
- 推送url到百度php,PHPCMS自动推送URL到百度站长平台
- 【训练题】航线设计 | 使用最长上升子序列(LIS)长度的O(nlogn)算法优化
- UHD630核显驱动方法及驱动后闪屏严重问题解决记录
- Arduino开源硬件-达芬奇-Level_1 (Mixly编程)
- 驾照理论考试android应用
- 两个经纬度点之间计算距离【经纬度距离计算】
- 架构师成长之路 5 --如何获取知识(学习境界)(方法)