转录组入门7-用DESeq2进行差异表达分析

任务

这个步骤推荐在R里面做,载入表达矩阵,然后设置好分组信息,统一用DEseq2进行差异分析,当然也可以走走edgeR或者limma的voom流程。
基本任务是得到差异分析结果,进阶任务是比较多个差异分析结果的异同点。

软件介绍

差异基因表达分析通常都是在R软件中,用各种R包进行分析的,常用的有DESeq2,edgeR,limma等几种。这次主要介绍用DESeq2来进行差异表达分析

关于使用DESeq2进行差异基因分析的介绍:

bioconductor.org/packag

dwheelerau.com/2014/02/

bio-info-trainee.com/bi

推荐大家先把以上三个教程学完,至少要把bioconductor上的说明(第一个链接)看完。

安装DESeq2(在R中进行操作)

source("https://bioconductor.org/biocLite.R") #载入安装工具
biocLite("DESeq2")  #安装包
library("DESeq2") #测试是否安装成功

DESeq2包的安装经常会出现错误。我尝试在windows、ubuntu和centos三种系统下安装这个包,都没有一次成功的。而且安装前我都已经使用bioconda安装过r-essentials(biodonda可以自动配置各种包的依赖环境),但是依然遇到各种不同的错误。根据我的经验,建议大家在遇到安装R包的问题时可以从以下方面解决:

1.尽量科学上网。各种R包常常会不能正常的下载,科学上网以后,这个问题会得到一定程度的解决。

2.注意R中的报错信息。如果科学上网后仍然不能解决R包安装的问题,就要认真注意R的报错信息,确认一下是哪种问题。

3.利用搜索引擎,搜索报错信息和解决办法。这实际上也需要科学上网,或者用必应国际版,查找解决办法。

DESeq2对于输入数据的要求

1.DEseq2要求输入数据是由整数组成的矩阵。

2.DESeq2要求矩阵是没有标准化的。

DESeq2进行差异表达分析

DESeq2包分析差异表达基因简单来说只有三步:构建dds矩阵,标准化,以及进行差异分析。

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) #~在R里面用于构建公式对象,~左边为因变量,右边为自变量。
dds <- DESeq(dds) #标准化
res <- results(dds, contrast=c("condition","treated","control")) #差异分析结果

构建dds矩阵

构建dds矩阵的代码就是快速入门中的第一行代码:

 dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) #~在R里面用于构建公式对象,~左边为因变量,右边为自变量。

构建dds矩阵需要:

表达矩阵 即上述代码中的countData,就是我们前面通过read count计算后并融合生成的矩阵,行为各个基因,列为各个样品,中间为计算reads或者fragment得到的整数。我们后面要用的是这个文件(mouse_all_count.txt)

样品信息矩阵即上述代码中的colData,它的类型是一个dataframe(数据框),第一列是样品名称,第二列是样品的处理情况(对照还是处理等),即condition,condition的类型是一个factor。这些信息可以从mouse_all_count.txt中导出,也可以自己单独建立。

差异比较矩阵 即上述代码中的design。 差异比较矩阵就是告诉差异分析函数是要从要分析哪些变量间的差异,简单说就是说明哪些是对照哪些是处理。

正式构建dds矩阵

library(DESeq2)  #加载包
countData <- mouse_all_count
condition <- factor(c("control","control","KD","KD"))
colData <- data.frame(row.names=colnames(countData), condition)

我们可以看一下condition和colData的具体内容:

#正式构建dds矩阵
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )
head(dds)  #查看一下构建好的矩阵

可以看到rownames都是原来mouse_all_count.txt中的行号,我们也可以在构建dds矩阵的时候,把行号直接替换成gene_id,这样操作后生成的dds矩阵如下:

dds <- DESeq(dds)   #对原始dds进行normalize

相关说明: DESeq包含三步,estimation of size factors(estimateSizeFactors), estimation of dispersion(estimateDispersons), Negative Binomial GLM fitting and Wald statistics(nbinomWaldTest),可以分布运行,也可用一步到位,最后返回 results可用的DESeqDataSet对象。

# 查看结果的名称,本次实验中是 "Intercept","condition_akap95_vs_control"
resultsNames(dds2)
# 将结果用results()函数来获取,赋值给res变量
res <- results(dds2)
# summary一下,看一下结果的概要信息
summary(res)

提取差异分析结果

# 获取padj(p值经过多重校验校正后的值)小于0.05,表达倍数取以2为对数后大于1或者小于-1的差异表达基因。
> table(res$padj<0.05) #取P值小于0.05的结果
> res <- res[order(res$padj),]
> diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
> diff_gene_deseq2 <- row.names(diff_gene_deseq2)
> resdata <-  merge(as.data.frame(res),as.data.frame(counts(dds2,normalize=TRUE)),by="row.names",sort=FALSE)
# 得到csv格式的差异表达分析结果
> write.csv(resdata,file= "H:/baidu/Rworkdir/count_matrix2/control_vs_akap95.cvs",row.names = F)

保存好的csv文件可以用于进一步绘图展示分析结果,我们下一次再具体介绍。


感谢:

生信技能书的转录组入门活动:RNA-seq基础入门传送门-转录组-生信技能树

hoptop的《(伪)从零开始学转录组》

panda姐的《PANDA姐的转录组入门》

金建峰的《浙大植物学小白的转录组笔记》


PS:从上次更新完转录组入门6(2个月前),我就一直纠缠在安装DESeq2以及在不同的电脑上部署各种生信分析环境,时常感叹要是能够有类似一键Ghost那样的简单快捷的安装方法或者软件来布置常用的生信分析环境就好了。这个问题其实也有解,一种方法是docker。但是docker的安装需要管理权限,而且常常有很多docker无法正常下载。另外一种方法是 @hoptop 在微信公众号@生信媛 介绍过的,主要是通过conda批量安装需要的程序。有机会我也会总结一下这方面的学习心得。

发布于 2017-10-22

文章被以下专栏收录