简单的转录组基因差异表达分析 -- edgeR
经典的转录组差异分析通常会使用到三个工具limma/voom, edgeR和DESeq2， 今天我们同样使用一个小规模的转录组测序数据来演示edgeR的简单流程。
由于，edgeR和DESeq2都是使用基于负二项分布的广义线性回归模型（GLM）来对RNA-seq数据进行拟合和差异分析，所以我们都用同一个数据来分析。
文中使用的数据来自Standford 大学的一个拟南芥的small RNA-seq数据（https://bios221.stanford.edu/data/mobData.RData）
该数据包含6个样本：SL236，SL260，SL237，SL238，SL239，SL240， 并分成了三组，分别是：
MM="triple mutatnt shoot grafted onto triple mutant root"
WM="wild-type shoot grafted onto triple mutant root"
WW="wild-type shoot grafted onto wild-type root"
简而言之，
WW组可以认为是实验的对照组，而MM和WM则是两个处理组。
P.S. 本文的分析是基于有生物学重复的单因子差异分析，关于无生物学重复或者多因子的情况，以后有机会再做展开。
对于edgeR的分析流程而言，我们需要输入的数据包括：
- 表达矩阵（
counts） - 分组信息（
group） - 拟合信息（
design）：指明如何根据样本的分组进行建模
下面就以mobData 中的数据为例简单介绍edgeR 的分析流程
载入数据及生成DGEList
由于mobData 中的行名没有提供基因的ID，我们也不是为了探究生物学问题，就以mobData 的行数作为其ID
library(edgeR)
load("data/mobData.RData")
head(mobData)
## SL236 SL260 SL237 SL238 SL239 SL240
## [1,] 21 52 4 4 86 68
## [2,] 18 21 1 5 1 1
## [3,] 1 2 2 3 0 0
## [4,] 68 87 270 184 396 368
## [5,] 68 87 270 183 396 368
## [6,] 1 0 6 10 6 12
row.names(mobData) <- as.character(c(1:dim(mobData)[1]))
# MM="triple mutatnt shoot grafted onto triple mutant root"
# WM="wild-type shoot grafted onto triple mutant root"
# WW="wild-type shoot grafted onto wild-type root"
mobGroups <- c("MM", "MM", "WM", "WM", "WW", "WW")
degs <- DGEList(counts = mobData, group = mobGroups);degs
## An object of class "DGEList"
## $counts
## SL236 SL260 SL237 SL238 SL239 SL240
## 1 21 52 4 4 86 68
## 2 18 21 1 5 1 1
## 3 1 2 2 3 0 0
## 4 68 87 270 184 396 368
## 5 68 87 270 183 396 368
## 2995 more rows ...
##
## $samples
## group lib.size norm.factors
## SL236 MM 152461 1
## SL260 MM 309995 1
## SL237 WM 216924 1
## SL238 WM 208841 1
## SL239 WW 258404 1
## SL240 WW 276434 1
edgeR将数据存储在列表形式的DGEList对象中，需要指定的参数包括counts 和group . DGEList 将表达矩阵存储在$counts中，将样本的信息，例如分组情况和文库大小等存储在$samples中。
预处理
在进行差异分析之前，需要对样本数据的表达矩阵进行预处理，包括：
- 去除低表达量基因
- 探索样本分组信息 -- 有助于挖掘潜在的差异样本
这里我们根据CPM normalization后的基因表达量作为过滤低表达基因的指标
# cpm normalization
countsPerMillion <- cpm(degs)
countCheck <- countsPerMillion > 1
keep <- which(rowSums(countCheck) >= 2)
degs.keep <- degs[keep,]
dim(degs.keep)
## [1] 2861 6
edgeR默认使用 trimmed mean of M-values (TMM) 计算文库的scale factor进行normalization，以最大程度地缩小样本间基因表达量的log-fold change。这是因为TMM 法认为样本间大部分的基因都没有发生差异表达，而那些真正差异表达的基因并不会受到normalization的严重影响。如此一来，便将那些由于测序引起的差异表达基因的表达量给校正了，消除了一部分的假阳性。
degs.norm <- calcNormFactors(degs.keep, method = 'TMM')
plotMDS(degs.norm, col=as.numeric(degs.norm$samples$group))
legend("bottomleft",as.character(unique(degs.norm$samples$group)), col=1:3, pch=20)
这里使用plotMDS 查看样本的分组情况（通过logFC），各组都分得很开。plotMDS在多因子的情况下可以更好地观察各个样本组是否有良好的分组。
关于Normalization
在差异分析中，我们常常更关注的是相对表达量的变化，例如处理组的A基因表达量相对于对照组的而言是上调还是下调了。而基因表达量的定量准确性则在差异分析中不太重要，因此，在进行差异分析时，像RPKM/FPKM这种对转录本长度进行normalization方法是并不常用，也是没有必要的。
在常规的RNA-seq中，影响基因表达量更大的技术因素往往是测序深度以及有效文库大小（effective libraries size）。这也是一般的差异分析软件会进行normalize的部分。
差异分析
首先，我们构建出design矩阵，指明差异分析所要比较的关系
designMat <- model.matrix(~0+mobGroups);designMat
## mobGroupsMM mobGroupsWM mobGroupsWW
## 1 1 0 0
## 2 1 0 0
## 3 0 1 0
## 4 0 1 0
## 5 0 0 1
## 6 0 0 1
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$mobGroups
## [1] "contr.treatment"
然后，进行dispersion的估计
degs.norm <- estimateGLMCommonDisp(degs.norm,design=designMat)
degs.norm <- estimateGLMTrendedDisp(degs.norm, design=designMat)
degs.norm <- estimateGLMTagwiseDisp(degs.norm, design=designMat)
plotBCV(degs.norm)
good fitting
plotBCV 反映不同表达量的基因与模型拟合的情况，如果模型拟合得好则"Tagwise"点的分布会拟合到“Trend”这条曲线上，如上图所展示的情况。但也可以看到低表达量的基因有点离散。
关于dispersion estimation
一般而言，样本间的变异系数（coefficient of variance，CV）是由两部分组成的，一是技术差异（Technical CV），另一个是生物学差异（Biological coefficient of variance，BCV）。前者是会随着测序通量的提升而消失的，而后者则是样本间真实存在的差异。所以，对于一个基因g而言，它的BCV在样本间足够大的话，就可以认为基因g是一个差异表达基因。而
edgeR正是通过估计dispersion来估计BCV（其中的数理不在此展开），进而拟合出线性回归模型的参数。
estimateGLMCommonDisp(x,design)：为所有基因都计算同样的dispersion
estimateGLMTrendedDisp(x,design)：根据不同基因的均值--方差之间的关系来计算dispersion，并拟合一个trended model
estimateGLMTagwiseDisp(x,design)：计算每个基因的dispersion，并利用经验性贝叶斯方法（empirical bayes）压缩至Trend Didspersion。个人认为这一项相当于GLM中每个基因的beta值
最后，根据design进行拟合，并利用likehood ratio test（LRT）进行统计检验
fit <- glmFit(degs.norm, designMat)
# LRT=likehood ratio test
# group1-group2
lrt.1vs2 <- glmLRT(fit, contrast = c(1,-1,0))
# group1-group3
lrt.1vs3 <- glmLRT(fit, contrast = c(1,0,-1))
# group2-group3
lrt.2vs3 <- glmLRT(fit, contrast = c(0,1,-1))
下面以MM组和 WW组的比较为例
topTags 提取出差异分析的数据；
decideTestsDGE 可以根据条件筛选差异基因，返回-1,0,1三个数值，分别代表下调，不显著和上调；
plotSmear 画一个简单的表达量与fold change的关系图。
degs.res.1vs3 <- topTags(lrt.1vs3, n = Inf, adjust.method = 'BH', sort.by = 'PValue')
degs.res.1vs3[1:5, ]
## Coefficient: 1*mobGroupsMM -1*mobGroupsWW
## logFC logCPM LR PValue FDR
## 74 -10.364020 9.042115 130.95167 2.537092e-30 7.258620e-27
## 490 -6.043444 8.401692 119.28833 9.056111e-28 1.295477e-24
## 1717 -7.056255 9.921304 114.37749 1.077199e-26 1.027289e-23
## 1963 6.492175 6.868420 102.37925 4.584800e-24 3.279278e-21
## 1111 -9.565662 8.284244 97.26908 6.051792e-23 3.462836e-20
deGenes.1vs3 <- decideTestsDGE(lrt.1vs3, p=0.05, lfc = 1)
summary(deGenes.1vs3)
## 1*mobGroupsMM -1*mobGroupsWW
## Down 430
## NotSig 2094
## Up 337
detag <- rownames(lrt.1vs3)[as.logical(deGenes.1vs3)]
plotSmear(lrt.1vs3, de.tags=detag)
abline(h=c(-1, 1), col='blue')
图中红色的是统计学上的显著差异表达基因
至于edgeR与DESeq2的比较其实已经有很多benchmark的文献做过了，这里就先鸽一下，以后有机会再来填坑。
benchmark ref：https://bioconductor.org/packages/release/bioc/vignettes/SummarizedBenchmark/inst/doc/Feature-Iterative.html
参考文章：
完。
