基因差异表达分析——基于RSEM对比,DESeq2操作实例

以下操作全部基于R中的DESeq2函数包

使用DESeq2的两点要求:

  1. DEseq2要求输入数据是由整数组成的矩阵。
  2. 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矩阵

  1. 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)
> coldata
                                       condition
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)
  1. 构建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 columns
                                 baseMean     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.45683300826916
                                      pvalue                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

以上是做基因差异表达分析的基本流程
完成!

  • 6
    点赞
  • 41
    收藏
    觉得还不错? 一键收藏
  • 5
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 5
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值