經(jīng)典的轉(zhuǎn)錄組差異分析通常會(huì)使用到三個(gè)工具limma/voom, edgeR和DESeq2， 今天我們同樣使用一個(gè)小規(guī)模的轉(zhuǎn)錄組測(cè)序數(shù)據(jù)來(lái)演示edgeR的簡(jiǎn)單流程。
由于，edgeR和DESeq2都是使用基于負(fù)二項(xiàng)分布的廣義線性回歸模型（GLM）來(lái)對(duì)RNA-seq數(shù)據(jù)進(jìn)行擬合和差異分析，所以我們都用同一個(gè)數(shù)據(jù)來(lái)分析。
文中使用的數(shù)據(jù)來(lái)自Standford 大學(xué)的一個(gè)擬南芥的small RNA-seq數(shù)據(jù)（https://bios221.stanford.edu/data/mobData.RData）
該數(shù)據(jù)包含6個(gè)樣本：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"
簡(jiǎn)而言之，
WW組可以認(rèn)為是實(shí)驗(yàn)的對(duì)照組，而MM和WM則是兩個(gè)處理組。
P.S. 本文的分析是基于有生物學(xué)重復(fù)的單因子差異分析，關(guān)于無(wú)生物學(xué)重復(fù)或者多因子的情況，以后有機(jī)會(huì)再做展開(kāi)。
對(duì)于edgeR的分析流程而言，我們需要輸入的數(shù)據(jù)包括：
- 表達(dá)矩陣（
counts） - 分組信息（
group） - 擬合信息（
design）：指明如何根據(jù)樣本的分組進(jìn)行建模
下面就以mobData 中的數(shù)據(jù)為例簡(jiǎn)單介紹edgeR 的分析流程
載入數(shù)據(jù)及生成DGEList
由于mobData 中的行名沒(méi)有提供基因的ID，我們也不是為了探究生物學(xué)問(wèn)題，就以mobData 的行數(shù)作為其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將數(shù)據(jù)存儲(chǔ)在列表形式的DGEList對(duì)象中，需要指定的參數(shù)包括counts 和group . DGEList 將表達(dá)矩陣存儲(chǔ)在$counts中，將樣本的信息，例如分組情況和文庫(kù)大小等存儲(chǔ)在$samples中。
預(yù)處理
在進(jìn)行差異分析之前，需要對(duì)樣本數(shù)據(jù)的表達(dá)矩陣進(jìn)行預(yù)處理，包括：
- 去除低表達(dá)量基因
- 探索樣本分組信息 -- 有助于挖掘潛在的差異樣本
這里我們根據(jù)CPM normalization后的基因表達(dá)量作為過(guò)濾低表達(dá)基因的指標(biāo)
# cpm normalization
countsPerMillion <- cpm(degs)
countCheck <- countsPerMillion > 1
keep <- which(rowSums(countCheck) >= 2)
degs.keep <- degs[keep,]
dim(degs.keep)
## [1] 2861 6
edgeR默認(rèn)使用 trimmed mean of M-values (TMM) 計(jì)算文庫(kù)的scale factor進(jìn)行normalization，以最大程度地縮小樣本間基因表達(dá)量的log-fold change。這是因?yàn)門(mén)MM 法認(rèn)為樣本間大部分的基因都沒(méi)有發(fā)生差異表達(dá)，而那些真正差異表達(dá)的基因并不會(huì)受到normalization的嚴(yán)重影響。如此一來(lái)，便將那些由于測(cè)序引起的差異表達(dá)基因的表達(dá)量給校正了，消除了一部分的假陽(yáng)性。
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 查看樣本的分組情況（通過(guò)logFC），各組都分得很開(kāi)。plotMDS在多因子的情況下可以更好地觀察各個(gè)樣本組是否有良好的分組。
關(guān)于Normalization
在差異分析中，我們常常更關(guān)注的是相對(duì)表達(dá)量的變化，例如處理組的A基因表達(dá)量相對(duì)于對(duì)照組的而言是上調(diào)還是下調(diào)了。而基因表達(dá)量的定量準(zhǔn)確性則在差異分析中不太重要，因此，在進(jìn)行差異分析時(shí)，像RPKM/FPKM這種對(duì)轉(zhuǎn)錄本長(zhǎng)度進(jìn)行normalization方法是并不常用，也是沒(méi)有必要的。
在常規(guī)的RNA-seq中，影響基因表達(dá)量更大的技術(shù)因素往往是測(cè)序深度以及有效文庫(kù)大?。╡ffective libraries size）。這也是一般的差異分析軟件會(huì)進(jìn)行normalize的部分。
差異分析
首先，我們構(gòu)建出design矩陣，指明差異分析所要比較的關(guān)系
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"
然后，進(jìn)行dispersion的估計(jì)
degs.norm <- estimateGLMCommonDisp(degs.norm,design=designMat)
degs.norm <- estimateGLMTrendedDisp(degs.norm, design=designMat)
degs.norm <- estimateGLMTagwiseDisp(degs.norm, design=designMat)
plotBCV(degs.norm)
plotBCV 反映不同表達(dá)量的基因與模型擬合的情況，如果模型擬合得好則"Tagwise"點(diǎn)的分布會(huì)擬合到“Trend”這條曲線上，如上圖所展示的情況。但也可以看到低表達(dá)量的基因有點(diǎn)離散。
關(guān)于dispersion estimation
一般而言，樣本間的變異系數(shù)（coefficient of variance，CV）是由兩部分組成的，一是技術(shù)差異（Technical CV），另一個(gè)是生物學(xué)差異（Biological coefficient of variance，BCV）。前者是會(huì)隨著測(cè)序通量的提升而消失的，而后者則是樣本間真實(shí)存在的差異。所以，對(duì)于一個(gè)基因g而言，它的BCV在樣本間足夠大的話，就可以認(rèn)為基因g是一個(gè)差異表達(dá)基因。而
edgeR正是通過(guò)估計(jì)dispersion來(lái)估計(jì)BCV（其中的數(shù)理不在此展開(kāi)），進(jìn)而擬合出線性回歸模型的參數(shù)。
estimateGLMCommonDisp(x,design)：為所有基因都計(jì)算同樣的dispersion
estimateGLMTrendedDisp(x,design)：根據(jù)不同基因的均值--方差之間的關(guān)系來(lái)計(jì)算dispersion，并擬合一個(gè)trended model
estimateGLMTagwiseDisp(x,design)：計(jì)算每個(gè)基因的dispersion，并利用經(jīng)驗(yàn)性貝葉斯方法（empirical bayes）壓縮至Trend Didspersion。個(gè)人認(rèn)為這一項(xiàng)相當(dāng)于GLM中每個(gè)基因的beta值
最后，根據(jù)design進(jìn)行擬合，并利用likelihood ratio test（LRT）進(jìn)行統(tǒng)計(jì)檢驗(yàn)
fit <- glmFit(degs.norm, designMat)
# LRT=likelihood 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 提取出差異分析的數(shù)據(jù)；
decideTestsDGE 可以根據(jù)條件篩選差異基因，返回-1,0,1三個(gè)數(shù)值，分別代表下調(diào)，不顯著和上調(diào)；
plotSmear 畫(huà)一個(gè)簡(jiǎn)單的表達(dá)量與fold change的關(guān)系圖。
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')
圖中紅色的是統(tǒng)計(jì)學(xué)上的顯著差異表達(dá)基因
至于edgeR與DESeq2的比較其實(shí)已經(jīng)有很多benchmark的文獻(xiàn)做過(guò)了，這里就先鴿一下，以后有機(jī)會(huì)再來(lái)填坑。
benchmark ref：https://bioconductor.org/packages/release/bioc/vignettes/SummarizedBenchmark/inst/doc/Feature-Iterative.html
參考文章：
完。
