简单使用DESeq2/EdgeR做差异分析

最新推荐文章于 2024-07-02 13:57:53 发布

刘永鑫Adam

最新推荐文章于 2024-07-02 13:57:53 发布

阅读量4.9k

点赞数 3

文章标签：机器学习大数据数据分析数据挖掘人工智能

原文链接：https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA#rd

版权

DESeq2和EdgeR都可用于做基因差异表达分析，主要也是用于RNA-Seq数据，同样也可以处理类似的ChIP-Seq,shRNA以及质谱数据。

这两个都属于R包，其相同点在于都是对count data数据进行处理，都是基于负二项分布模型。因此会发现，用两者处理同一组数据，最后在相同阈值下筛选出的大部分基因都是一样的，但是有一部分不同应该是由于其估计离散度的不同方法所导致的。 ### DESeq2的使用方法：

输入矩阵数据，行名为sample，列名为gene；DESeq2不支持无生物学重复的数据，因此我选择了2个样本，3个生物学重复的数据；并对count data取整（经大神指点，这里需要说明下，我的测试数据readcount是RSEM定量的结果，并不是常见的htseq-count的结果，所以count值会有小数点，而DESeq2包不支持count数有小数点，所以这里需要round取整）。
```
 database_all <- read.table(file = "readcount", sep = "\t", header = T, row.names = 1)
 database <- database_all[,1:6]
 type <- factor(c(rep("LC_1",3), rep("LC_2",3)))
 database <- round(as.matrix(database))
```

设置分组信息以及构建dds对象

 condition <- factor(c(rep("LC_1",3), rep("LC_2",3)))
 coldata <- data.frame(row.names = colnames(database), condition)
 dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)

使用DESeq函数进行估计离散度，然后进行标准的差异表达分析，得到res对象结果
```
 dds <- DESeq(dds)
 res <- results(dds)
```

最后设定阈值，筛选差异基因，导出数据

 table(res$padj <0.05)
 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 = "LC_1_vs_LC_2.csv")

EdgeR的使用方法：

跟DESeq2一样，EdgeR输入矩阵数据，行名为sample，列名为gene；DESeq2不支持无生物学重复的数据，因此我选择了2个样本，3个生物学重复的数据。

 exprSet_all <- read.table(file = "readcount", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
 exprSet <- exprSet_all[,1:6]
 group_list <- factor(c(rep("LC_1",3), rep("LC_2",3)))

设置分组信息，去除低表达量的gene以及做TMM标准化

 exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]
 exprSet <- DGEList(counts = exprSet, group = group_list)
 exprSet <- calcNormFactors(exprSet)

使用qCML（quantile-adjusted conditional maximum likelihood）估计离散度（只针对单因素实验设计）
```
 exprSet <- estimateCommonDisp(exprSet)
 exprSet <- estimateTagwiseDisp(exprSet)
```

寻找差异gene(这里的exactTest函数还是基于qCML并且只针对单因素实验设计)，然后按照阈值进行筛选即可

 et <- exactTest(exprSet)
 tTag <- topTags(et, n=nrow(exprSet))
 tTag <- as.data.frame(tTag)
 write.csv(tTag,file = "LC_1_vs_LC_2_edgeR.csv")

Summary

以上我主要针对单因素两两比较组进行差异分析，其实DESeq2和EdgeR两个R包都可以对多因素进行差异分析。

DESeq2修改以上代码的分组信息design参数以及在差异分析results函数中添加所选定的分组因素，其他代码基本一样，具体参照DESeq2手册
EdgeR则需要用Cox-Reid profile-adjusted likelihood (CR)方法来估算离散度，y <- estimateDisp(y, design)或者分别使用三个函数（y <- estimateGLMCommonDisp(y, design)，y <- estimateGLMTrendedDisp(y, design)，）y <- estimateGLMTagwiseDisp(y, design)；然后差异表达分析也跟单因素分析不同，主要使用generalized linear model (GLM) likelihood ratio test 或者 quasi-likelihood(QL) F-test，具体代码可以参照EdgeR手册。

http://www.bioinfo-scrounger.com

系列教程：微生物组入门 Biostar 微生物组宏基因组

专业技能：学术图表高分文章生信宝典不可或缺的人

一文读懂：宏基因组寄生虫益处进化树

必备技能：提问搜索 Endnote

文献阅读热心肠 SemanticScholar Geenmedical

扩增子分析：图表解读分析流程统计绘图

16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun

在线工具：16S预测培养基生信绘图

科研经验：云笔记云协作公众号

编程模板: Shell R Perl

生物科普: 肠道细菌人体上的生命生命大跃进细胞暗战人体奥秘

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。