RNA-seq 詳細教程：詳解DESeq2流程（9）
學習目標
- 瞭解
DESeq2涉及的不同步驟 - 瞭解變異的來源並檢查
size factors - 檢查基因水平的離散估計
- 瞭解差異表達分析過程中離散的重要性
DESeq2流程
前面，我們使用設計公式建立了 DESeq2 物件，並使用下面兩行程式碼執行DESeq2：
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)
dds <- DESeq(dds)
我們用 DESeq2 完成了差異基因表達分析的整個工作流程。分析中的步驟如下：
我們將詳細瞭解這些步驟中的每一個，以更好地瞭解 DESeq2
1. size factors
差異表達分析的第一步是估計大小因子，這正是我們已經對原始計數進行歸一化所做的。
DESeq2 在執行差異表達分析時會自動估計大小因子。但是，如果您已經像我們之前所做的那樣使用 estimateSizeFactors() 生成了大小因子，那麼 DESeq2 將使用這些值。
為了歸一化計數資料，DESeq2 使用前面教程中討論的比率中值方法計算每個樣本的大小因子。
-
MOV10 DE分析：檢查size factors
讓我們瀏覽一下每個樣本的大小因子值：
# Check the size factors sizeFactors(dds) Irrel_kd_1 Irrel_kd_2 Irrel_kd_3 Mov10_kd_2 Mov10_kd_3 Mov10_oe_1 Mov10_oe_2 1.1149694 0.9606733 0.7492240 1.5633640 0.9359695 1.2262649 1.1405026 Mov10_oe_3 0.6542030
這些數字應該與我們最初執行函式 estimateSizeFactors(dds) 時生成的數字相同。檢視每個樣本的總讀取數：
# Total number of raw counts per sample
colSums(counts(dds))
- 這些數字如何與尺寸因子相關聯？
我們看到較大的大小因子對應於具有較高測序深度的樣本，這是有道理的，因為要生成我們的歸一化計數，我們需要將計數除以大小因子。這解釋了樣本之間測序深度的差異。
現在看看歸一化後的總深度：
# Total number of normalized counts per sample colSums(counts(dds, normalized=T))
- 樣本間的值與每個樣本的總計數相比如何？
您可能期望歸一化後樣本中的計數完全相同。然而，DESeq2 還在歸一化過程中考慮了 RNA 組成。通過使用大小因子的中值比值，DESeq2 不應偏向於被少數 DE 基因吸收的大量計數；然而，這可能導致大小因素與僅基於測序深度的預期大不相同。
2. gene-wise dispersion
差異表達分析的下一步是估計基因差異。在我們進入細節之前，我們應該對 DESeq2 中的離散指的是什麼有一個很好的瞭解。
在 RNA-seq 計數資料中，我們知道：
- 為了確定差異表達的基因，我們評估組間表達的變化並將其與組內（重複之間）的變化進行比較。
- 對於每個單獨的基因，均值不等於方差。 高表達的基因將具有更一致的變異水平，但會高於平均值。 低表達的基因將表現出徘徊在平均值附近的變異（但具有更高的變異性）。
這種複雜的關係意味著我們不能只使用觀察到的方差來解釋組內變異。相反，DESeq2 使用離散。
- 什麼是離散（
dispersion）？
離散引數通過描述方差偏離均值的程度來模擬組內變異性。離散度為 1 表示沒有偏離均值（即均值 == 方差）。一個典型的 RNA-seq 資料集，將在重複中表現出一定數量的生物變異性，因此我們將始終具有小於 1 的離散值。
- 離散值是如何計算的？
在 DESeq 中，我們知道給定基因的計數方差由均值和離散度建模：
現在讓我們重新排列公式，以便我們可以看到離散引數等同於什麼，以便我們可以更好地理解它與均值和方差的關係：
這也與以下內容相同：
| Effect on dispersion | |
|---|---|
| Variance increases | Dispersion increases |
| Mean expression increases | Dispersion decreases |
-
DESeq2中的離散
DESeq2 根據基因的表達水平（組內重複的平均計數）和重複觀察到的方差來估計每個基因的離散度，正如我們用上面的公式所證明的那樣。這樣，具有相同均值的基因的離散估計將僅基於它們的方差而不同。因此，離散估計反映了給定平均值的基因表達的方差。
下面，有一個離散圖，其中每個黑點都是一個基因，離散是針對每個基因的平均表達繪製的。如上所述，您可以看到均值和離散之間的反比關係。黑點是根據我們擁有的資料進行的離散估計。每組只有少數 (3-6) 次重複，每個基因的變異估計通常不可靠。
為了解決這個問題，DESeq2 使用一種稱為 shrinkage 的方法在基因之間共享資訊，以根據基因的平均表達水平生成更準確的變異估計。 DESeq2 假定具有相似表達水平的基因應該具有相似的離散度。藍點代表縮小的離散值。
3. 擬合曲線
流程的下一步是將曲線擬合到基因方面的離散估計。將曲線擬合到資料背後的想法是，不同的基因將具有不同規模的生物變異性，但是，在所有基因中，將存在合理的離散估計分佈。
這條曲線在下圖中顯示為一條紅線，它繪製了給定表達強度的基因的預期離散值的估計值。每個黑點都是一個基因，具有相關的平均表達水平和離散的最大似然估計 (MLE)（步驟 1）。
4. Shrink
流程的下一步是將基因方面的離散估計縮小到預期的離散值。
當樣本量較小時，該曲線可以更準確地識別差異表達的基因，並且每個基因的收縮強度取決於：
- 基因離散離曲線有多近
- 樣本量（更多樣本 = 更少收縮）
這種收縮方法對於減少差異表達分析中的誤報尤為重要。具有低離散估計的基因向曲線收縮，並且輸出更準確、更高收縮值用於模型擬合和差異表達測試。這些縮小的估計值代表了確定跨組基因表達是否顯著不同所需的組內變異。
略高於曲線的離散估計也會向曲線收縮，以獲得更好的離散估計；然而，具有極高離散值的基因則不然。這是由於該基因可能不遵循建模假設，並且由於生物學或技術原因比其他基因具有更高的變異性。向曲線收縮值可能會導致誤報，因此這些值不會收縮。這些基因顯示在下面的藍色圓圈中。
這是一個很好的檢查方式，以確保您的資料非常適合 DESeq2 模型。您希望您的資料通常散佈在曲線周圍，散佈隨著平均表達水平的增加而降低。如果您看到雲或不同的形狀，那麼您可能想要更多地探索您的資料以檢視您是否有汙染（線粒體等）或異常樣本。請注意，對於任何低自由度的實驗，在 plotDispEsts() 圖中的整個均值範圍內收縮了多少。
令人擔憂的離散圖示例如下所示：
下圖顯示了離散值雲，這些值通常不遵循曲線。這會令人擔憂，並表明資料與模型的擬合不佳。
下圖顯示離散值最初下降，然後隨著較大的表達值而增加。根據我們的預期，較大的平均表達值不應該有較大的離散——我們期望離散隨著均值的增加而減小。這表明比預期的更高度表達的基因的變異更少。這也表明我們的分析中可能存在異常樣本或汙染。
5. MOV10
讓我們看一下 MOV10 資料的離散估計：
# Plot dispersion estimates
plotDispEsts(dds)
由於我們的樣本量很小，因此對於許多基因，我們看到了相當多的收縮。那麼我們的資料適合該模型嗎？
隨著平均表達的增加，我們看到離散度有很好的降低，這是好的。我們還看到離散估計通常圍繞曲線，這也是預期的。總的來說，看起來不錯。我們確實看到了強烈的收縮，這可能是因為我們的一個樣本組只有兩個重複。我們擁有的重複越多，對離散估計應用的收縮越少，能夠識別的 DE 基因就越多。我們通常建議每個條件至少進行 4 次生物重複，以便更好地估計變異。
歡迎Star -> 學習目錄
更多教程 -> 學習目錄
本文由mdnice多平臺釋出