如何采用DESeq2對表達量進行PCA和聚類分析

如何采用DESeq2對表達量進行PCA和聚類分析，相信很多沒有經驗的人對此束手無策，為此本文總結了問題出現的原因和解決方法，通過這篇文章希望你能解決這個問題。

得到基因/轉錄本的表達量之后，通常會通過以下三種類型的圖表來檢驗和分析生物學樣本和實驗設計間關系。

1.  樣本的聚類樹

利用所有樣本的表達量數據，對樣本進行聚類。理論上如果樣本和實驗操作都沒有問題，那么屬于同一組的生物學重復樣本會聚到一起。示意圖如下

上圖中，樣本的名稱用組別代替，可以看到，同一條件的樣本聚在了一起。

2. PCA圖

通過主成分分析進行降維，在二維或者三維平面上展示樣本點的分布，根據點的位置，也可以看出屬于同一組的樣本是否在一起，不同組之間的樣本有沒有明顯分開，示意如下

從圖中可以看到，不同條件的樣本區分的很明顯，而生物學重復之間距離較近，表明生物學重復的一致性和不同分組的差異性較好。

3.  熱圖

相比樣本的聚類樹，熱圖包含了更多的信息，比如可以直觀的展示不同分組間表達量的差異，也是常見的可視化手段之一，示意如下

只要有樣本的表達量矩陣，DESeq2可以輕松的畫出以上3種圖表。但是我們應該選擇原始的表達量矩陣，還是歸一化之后的表達量矩陣來畫呢？或者有沒有其他的選擇呢？

輸入的矩陣不同，得出的結論也會不同。由于基因的表達水平在不同樣本間本身就存在一定的差異，所以無論是采用原始的還是歸一化之后的表達量矩陣，效果都不理想。針對這一問題，DESeq2提出了兩種count值的轉換算法，rlog和VST轉換。

1. rlog 轉換

rlog 轉換的用法如下

rld <- rlog(dds)

2. VST 轉換

用法如下

vsd <- vst(dds)

兩種轉換本質上是在降低生物學重復之間的差異，使得樣本聚類和PCA分析的效果更好。轉換之后的表達量數據可以采用assay函數進行提取，代碼如下

> head(assay(rld)[, 1:2])
       sample1   sample2
gene1 2.049029 1.6828707
gene2 8.151262 6.8552583
gene3 0.818971 0.2964686
gene4 5.340361 4.4766682
gene5 6.316175 6.8345783
gene6 2.157821 1.9264385

對于raw count定量表格，建議采用rlog或者VST轉換之后的數據去進行PCA和聚類分析，效果會更好。

利用DESeq2提供的示例數據pasilla，分別用原始的count, 歸一化之后的count, rlog, vst 轉換的count 進行PCA分析，代碼如下

dds <- estimateSizeFactors(dds)
raw <- SummarizedExperiment(counts(dds, normalized=FALSE),
                                colData=colData(dds))
nor <- SummarizedExperiment(counts(dds, normalized=TRUE),
                                colData=colData(dds))
vsd <- vst(dds)
rld <- rlog(dds)
pdf("PCA.pdf")
plotPCA( DESeqTransform(raw), intgroup=c("condition", "type") )
plotPCA( DESeqTransform(nor), intgroup=c("condition", "type") )
plotPCA(vsd, intgroup=c("condition", "type"))
plotPCA(rld, intgroup=c("condition", "type"))
dev.off()

raw count 的結果如下

歸一化之后count結果如下

VST轉換之后的結果如下

rlog轉換之后的結果如下

可以很明顯看出，原始的count和歸一化之后的count, 其PCA圖是雜亂無序的，沒什么明顯規律，而VST和rlog轉換之后，生物學重復之間更佳的接近，不同分組也區分的較為明顯。

看完上述內容，你們掌握如何采用DESeq2對表達量進行PCA和聚類分析的方法了嗎？如果還想學到更多技能或想了解更多相關內容，歡迎關注億速云行業資訊頻道，感謝各位的閱讀！