對(duì)樣本的無(wú)監(jiān)督聚類

在我們看來(lái)，用于檢查基因表達(dá)分析的最重要的探索性圖表之一便是MDS圖或其余類似的圖。這種圖表使用無(wú)監(jiān)督聚類方法展示出了樣品間的相似性和不相似性，能讓我們?cè)谶M(jìn)行正式的檢驗(yàn)之前對(duì)于能檢測(cè)到多少差異表達(dá)基因有個(gè)大致概念。理想情況下，樣本會(huì)在不同的實(shí)驗(yàn)組內(nèi)很好的聚類，且可以鑒別出遠(yuǎn)離所屬組的樣本，并追蹤誤差或額外方差的來(lái)源。如果存在技術(shù)重復(fù)，它們應(yīng)當(dāng)互相非常接近。

這樣的圖可以用limma中的plotMDS函數(shù)繪制。第一個(gè)維度表示能夠最好地分離樣品且解釋最大比例的方差的引導(dǎo)性的倍數(shù)變化（leading-fold-change），而后續(xù)的維度的影響更小，并與之前的維度正交。當(dāng)實(shí)驗(yàn)設(shè)計(jì)涉及到多個(gè)因子時(shí)，建議在多個(gè)維度上檢查每個(gè)因子。如果在其中一些維度上樣本可按照某因子聚類，這說(shuō)明該因子對(duì)于表達(dá)差異有影響，在線性模型中應(yīng)當(dāng)將其包括進(jìn)去。反之，沒(méi)有或者僅有微小影響的因子在下游分析時(shí)應(yīng)當(dāng)被剔除。

suppressMessages(library(RColorBrewer))
col.group <- pheno
levels(col.group) <- brewer.pal(nlevels(col.group), "Set1") 
col.group <- as.character(col.group)
plotMDS(lcpm, labels=pheno,col = col.group)
title(main="Sample groups")

MDS聚類.png

差異表達(dá)分析

主要說(shuō)明一下edgeR中的glmQLFTest,exactTest以及limma中的voom這幾種獲取差異基因的不同方式

1基于limma的差異分析

Limma包基于線性模型建模。 它最初設(shè)計(jì)用于分析微陣列數(shù)據(jù)，但最近已擴(kuò)展到RNA-seq數(shù)據(jù)。 根據(jù)limma用戶指南的當(dāng)前建議是使用edgeR包的TMM標(biāo)準(zhǔn)化和“voom”轉(zhuǎn)換，其本質(zhì)上將標(biāo)準(zhǔn)化數(shù)據(jù)取對(duì)數(shù)并估計(jì)它們的均值 - 方差關(guān)系以確定在線性建模之前每次觀察的權(quán)重。 默認(rèn)情況下，Benjamini-Hochberg程序用于估計(jì)FDR 。

首先先建立design矩陣,在此研究中，我們想知道哪些基因在我們研究的兩組之間以不同水平表達(dá)。在我們的分析中，假設(shè)基礎(chǔ)數(shù)據(jù)是正態(tài)分布的，為其擬合一個(gè)線性模型。在此之前，需要?jiǎng)?chuàng)建一個(gè)包含細(xì)胞類型design矩陣。

design <- model.matrix(~0+pheno)
colnames(design) <- gsub("pheno", "", colnames(design))
rownames(design) <- colnames(y)
#注意使用colnames()，而不能使用names(),后者只會(huì)得到：[1] "counts"  "samples"
design

據(jù)顯示對(duì)于RNA-seq計(jì)數(shù)數(shù)據(jù)而言，當(dāng)使用原始計(jì)數(shù)或當(dāng)其被轉(zhuǎn)換為log-CPM值時(shí)，方差并不獨(dú)立于均值。使用負(fù)二項(xiàng)分布來(lái)模擬計(jì)數(shù)的方法假設(shè)均值與方差間具有二次的關(guān)系。在limma中，假設(shè)log-CPM值符合正態(tài)分布，并使用由voom函數(shù)計(jì)算得到的精確權(quán)重來(lái)調(diào)整均值與方差的關(guān)系，從而對(duì)log-CPM值進(jìn)行線性建模。

當(dāng)操作DGEList對(duì)象時(shí)，voom從x中自動(dòng)提取文庫(kù)大小和歸一化因子，以此將原始計(jì)數(shù)轉(zhuǎn)換為log-CPM值。在voom中，對(duì)于log-CPM值額外的歸一化可以通過(guò)設(shè)定normalize.method參數(shù)來(lái)進(jìn)行。

下圖左側(cè)展示了這個(gè)數(shù)據(jù)集log-CPM值的均值-方差關(guān)系。通常而言，方差是測(cè)序?qū)嶒?yàn)中的技術(shù)差異和不同細(xì)胞類型的重復(fù)樣本之間的生物學(xué)差異的結(jié)合，而voom圖會(huì)顯示出一個(gè)在均值與方差之間遞減的趨勢(shì)。 生物學(xué)差異高的實(shí)驗(yàn)通常會(huì)有更平坦的趨勢(shì)，其方差值在高表達(dá)處穩(wěn)定。 生物學(xué)差異低的實(shí)驗(yàn)更傾向于急劇下降的趨勢(shì)。

不僅如此，voom圖也提供了對(duì)于上游所進(jìn)行的過(guò)濾水平的可視化檢測(cè)。如果對(duì)于低表達(dá)基因的過(guò)濾不夠充分，在圖上表達(dá)低的一端，受到非常低的表達(dá)計(jì)數(shù)的影響，可以觀察到方差水平的下降。如果觀察到了這種情況，應(yīng)當(dāng)回到最初的過(guò)濾步驟并提高用于該數(shù)據(jù)集的表達(dá)閾值。

par(mfrow=c(1,2))
v <- voom(y, design, plot=TRUE)

fit <- lmFit(v, design)

head(coef(fit),3)

contr <- makeContrasts(trt - untrt, levels = colnames(coef(fit)))
contr
#       Contrasts
#Levels  trt - untrt
# trt             1
# untrt          -1
diff <- contrasts.fit(fit, contr)
diff <- eBayes(diff)
plotSA(diff, main="Final model: Mean-variance trend")
top.table <- topTable(diff, sort.by = "P", n = Inf)
DEG <- na.omit(top.table)
summary(decideTests(diff))

前后數(shù)據(jù)對(duì)比

2基于edgeR的差異基因分析

edgeR使用經(jīng)驗(yàn)貝葉斯估計(jì)和基于負(fù)二項(xiàng)模型的精確檢驗(yàn)來(lái)確定差異基因。 特別地，經(jīng)驗(yàn)貝葉斯用于通過(guò)在基因之間來(lái)調(diào)節(jié)跨基因的過(guò)度離散程度。 使用類似于Fisher精確檢驗(yàn)但適應(yīng)過(guò)度分散數(shù)據(jù)的精確檢驗(yàn)用于評(píng)估每個(gè)基因的差異表達(dá)。edgeR 在默認(rèn)情況下，執(zhí)行TMM標(biāo)準(zhǔn)化程序以考慮樣本之間的不同測(cè)序深度，通過(guò)Benjamini-Hochberg用于控制FDR 。

利用exactTest估計(jì)差異基因

精確檢驗(yàn)適用于多組實(shí)驗(yàn)的精確統(tǒng)計(jì)法，使用函數(shù)exactTest估計(jì)差異基因，人們將其稱為classic edgeR。

estimateDisp函數(shù)在一組離散網(wǎng)格點(diǎn)上為每個(gè)標(biāo)簽(基因)計(jì)算一個(gè)似然矩陣，然后應(yīng)用加權(quán)似然經(jīng)驗(yàn)貝葉斯方法獲得后驗(yàn)離散度估計(jì)。如果沒(méi)有設(shè)計(jì)矩陣，它計(jì)算每個(gè)標(biāo)簽的分位數(shù)的條件似然，然后將其最大化。在這種情況下，它類似于函數(shù)estimateCommonDisp和estimateTagwiseDisp。

同樣首先利用calcNormFactors對(duì)因子進(jìn)行矯正

#對(duì)因子矯正
y <- calcNormFactors(y)
#對(duì)每個(gè)基因估測(cè)了一個(gè)經(jīng)驗(yàn)貝葉斯文件離散值，一個(gè)公共離散值，一個(gè)趨勢(shì)離散值
y <- estimateDisp(y,design = design)


### 等同于下面兩步
# 估計(jì)變異系數(shù)，即估計(jì)方差；估計(jì)內(nèi)部差異程度，看組間差異是否比內(nèi)部差異大，如果大，可選為差異基因
y <- estimateCommonDisp(y)
# 估計(jì)逐段分布
y <- estimateTagwiseDisp(y)

# 繪制基因的趨勢(shì)圖
plotBCV(y)

BCV

# 比較獲得差異基因
et <- exactTest(y,pair = c("normal","tumor"))

利用glmQLFTest估計(jì)差異基因

似然比檢驗(yàn)是用廣義線性模型（glms）的統(tǒng)計(jì)學(xué)方法，適用于不同復(fù)雜程度的多因素實(shí)驗(yàn)，而QLFTest則是首選，因?yàn)樗从沉斯烙?jì)每個(gè)基因離散度的不確定性。當(dāng)重復(fù)次數(shù)較少時(shí)，它可以提供更可靠的錯(cuò)誤率控制。

y <- calcNormFactors(y)
y <- estimateDisp(y, design)  
plotBCV(y)
fit <- glmQLFit(y, design, robust=TRUE)
plotQLDisp(fit)
con <- makeContrasts(trt- untrt, levels = design)
res <- glmQLFTest(fit, contrast=con)
summary(decideTestsDGE(res))

plotMD(res)
DEmarkers <- topTags(res, n=Inf)$table
head(DEmarkers)
is.sig <- DEmarkers$FDR <= 0.05
plotSmear(res, de.tags=rownames(DEmarkers)[is.sig], cex=0.1)

差異基因分布