RNA-seq 详细教程：搞定count归一化（5）
學習目標
1. 歸一化
差異表達分析工作流程的第一步是計數歸一化，這是對樣本之間的基因表達進行準確比較所必需的。
Normalization每個基因的映射讀數計數是 RNA 表達以及許多其他因素的結果。歸一化是調整原始計數值以解決“無關”因素的過程。以這種方式，表達水平在樣本之間或樣本內更具可比性。
在歸一化過程中經?？紤]的“無關”因素：
1.1. 測序深度
考慮測序深度對于比較樣本之間的基因表達是必要的。在下面的示例中，每個基因在樣本 A 中的表達似乎是樣本 B 的兩倍。然而，這是樣本 A 的測序深度加倍的結果。
sequencing depth1.2. 基因長度
計算基因長度對于比較同一樣本中不同基因之間的表達是必要的。在下面的示例中，基因 X 和基因 Y 具有相似的表達水平，但映射到基因 X 的讀數數量將比映射到基因 Y 的讀數多得多，因為基因 X 更長。
Gene length1.3. RNA組成
樣本之間一些高度差異表達的基因、樣本之間表達的基因數量的差異或污染的存在可能會扭曲某些類型的歸一化方法。建議考慮 RNA 組成以準確比較樣本之間的表達，這在進行差異表達分析時尤為重要。
在下面的示例中，假設樣本 A 和樣本 B 之間的測序深度相似，并且除了基因差異表達之外的每個基因在樣本之間呈現相似的表達水平。樣本 B 中的計數會受到 差異表達基因的極大影響，它占據了大部分計數。因此，樣本 B 的其他基因的表達似乎低于樣本 A 中的相同基因。
RNA composition歸一化不僅對于差異表達分析必不可少，對于探索數據分析、數據可視化以及探索或比較樣本之間或樣本內的計數也是必要的。
2. 歸一化方法
幾種常見的歸一化方法：
| CPM (counts per million) | 按照reads總數縮放計數 | 測序深度 | 同一樣本組重復之間的基因計數比較；不適用于樣本內比較或差異表達分析 |
| TPM (transcripts per kilobase million) | 每百萬讀取reads比對的轉錄本長度 (kb) 計數 | 測序深度與基因長度 | 樣本內或同一樣本組樣本之間的基因計數比較；不適用于差異表達分析 |
| RPKM/FPKM (reads/fragments per kilobase of exon per million reads/fragments mapped) | 類似于TPM | 測序深度與基因長度 | 樣本中基因之間的基因計數比較；不適用于樣本比較或差異表達分析 |
| DESeq2’s median of ratios | 計數除以特定于樣本的大小因子，該因子由基因計數相對于每個基因的幾何平均值的中位數比率確定 | 測序深度和RNA組成 | 樣品之間的基因計數比較和差異表達分析；不適用于樣本內比較 |
| EdgeR’s trimmed mean of M values (TMM) | 使用樣本之間對數表達比率的加權修剪平均值 | 測序深度和RNA組成 | 樣品之間的基因計數比較和差異表達分析；不適用于樣本內比較 |
- RPKM/FPKM：不推薦用于樣本間比較
雖然 TPM 和 RPKM/FPKM 歸一化方法都考慮了測序深度和基因長度，但不推薦使用 RPKM/FPKM。原因是RPKM/FPKM方法輸出的歸一化計數值在樣本之間沒有可比性。
使用 RPKM/FPKM 歸一化，每個樣本的 RPKM/FPKM 歸一化計數總數會有所不同。因此，您不能在樣本之間平均比較每個基因的歸一化計數。
RPKM-歸一化計數表：
| XCR1 | 5.5 | 5.5 |
| WASHC1 | 73.4 | 21.8 |
| … | … | … |
| Total RPKM-normalized counts | 1,000,000 | 1,500,000 |
例如，在上表中，樣本 A 的 XCR1 (5.5/1,000,000) 計數比例高于樣本 B (5.5/1,500,000)，即使 RPKM 計數值相同。因此，我們不能直接比較樣本 A 和樣本 B 之間 XCR1（或任何其他基因）的計數，因為樣本之間的歸一化計數總數不同。
- DESeq2-歸一化計數：比率方法的中值(Median of ratios method)
由于用于差異表達分析的工具正在比較樣本組之間相同基因的計數，因此該工具不需要考慮基因長度。然而，確實需要考慮測序深度和 RNA 組成。為了標準化測序深度和 RNA 組成，DESeq2 使用比率中值方法。在用戶端只有一個步驟，但在后端涉及多個步驟，如下所述。
對于每個基因，都會創建一個偽參考樣本，該樣本等于所有樣本的幾何平均值。
| EF2A | 1489 | 906 | sqrt(1489 * 906) = 1161.5 |
| ABCD1 | 22 | 13 | sqrt(22 * 13) = 17.7 |
| … | … | … | … |
對于每個樣本中的每個基因，計算比率（樣本/參考）（如下所示）。由于大多數基因沒有差異表達，因此每個樣本中的大多數基因在樣本中的比例應該相似。
| EF2A | 1489 | 906 | 1161.5 | 1489/1161.5 = 1.28 | 906/1161.5 = 0.78 |
| ABCD1 | 22 | 13 | 16.9 | 22/16.9 = 1.30 | 13/16.9 = 0.77 |
| MEFV | 793 | 410 | 570.2 | 793/570.2 = 1.39 | 410/570.2 = 0.72 |
| BAG1 | 76 | 42 | 56.5 | 76/56.5 = 1.35 | 42/56.5 = 0.74 |
| MOV10 | 521 | 1196 | 883.7 | 521/883.7 = 0.590 | 1196/883.7 = 1.35 |
| … | … | … | … |
給定樣本的所有比率的中值（上表中的列）被視為該樣本的歸一化因子（大小因子），計算如下。請注意，差異表達的基因不應影響中值：
normalization_factor_sampleA <- median(c(1.28, 1.3, 1.39, 1.35, 0.59))
normalization_factor_sampleB <- median(c(0.78, 0.77, 0.72, 0.74, 1.35))
下圖說明了單個樣本的所有基因比率分布的中值（y 軸是頻率）。
figure比率中位數法假設并非所有基因都差異表達；因此，歸一化因子應考慮樣本的測序深度和 RNA 組成（大的離群基因不會影響中值比率值）。該方法對上調/下調和大量差異表達基因的不平衡具有魯棒性。
這是通過將給定樣本中的每個原始計數值除以該樣本的歸一化因子來執行的，生成歸一化計數值。這是針對所有計數值（每個樣本中的每個基因）執行的。例如，如果樣本 A 的中值比率為 1.3，樣本 B 的中值比率為 0.77，則可以按如下方式計算歸一化計數：
- Raw Counts
| EF2A | 1489 | 906 |
| ABCD1 | 22 | 13 |
| … | … | … |
- Normalized Counts
| EF2A | 1489 / 1.3 = 1145.39 | 906 / 0.77 = 1176.62 |
| ABCD1 | 22 / 1.3 = 16.92 | 13 / 0.77 = 16.88 |
| … | … | … |
請注意，歸一化計數值不是整數。
”★
以上步驟僅作為演示，在實際使用DESeq2過程中，只需要一步命令，即可完成計算。
”3. Mov10 歸一化
現在我們知道了計數歸一化的理論，我們將使用 DESeq2 對 Mov10 數據集的計數進行歸一化。這需要幾個步驟：
3.1. 數據匹配
我們應該始終確保樣本名稱在兩個文件之間匹配，并且樣本的順序相同。如果不是這種情況，DESeq2 將輸出錯誤。
#?檢查兩個文件中的樣本名稱是否匹配all(colnames(txi$counts)?%in%?rownames(meta))
all(colnames(txi$counts)?==?rownames(meta))
如果數據不匹配，可以使用 match() 函數重新排列它們。
3.2. 創建對象
讓我們從創建 DESeqDataSet 對象開始，然后可以更多地討論其中存儲的內容。要創建對象，我們需要將計數矩陣和元數據表作為輸入。我們還需要指定一個設計公式。設計公式指定元數據表中的列以及它們在分析中的使用方式。對于我們的數據集，我們只有一列感興趣，即 ~sampletype。此列具有三個因子水平，它告訴 DESeq2 對于每個基因，我們要評估相對于這些不同水平的基因表達變化。
我們的計數矩陣輸入存儲在 txi 列表對象中。所以我們需要指定使用 DESeqDataSetFromTximport() 函數，這將提取計數組件并將值四舍五入到最接近的整數。
#?對象創建dds?<-?DESeqDataSetFromTximport(txi,?colData?=?meta,?design?=?~?sampletype)
3.3. 生成歸一化counts
下一步是對計數數據進行歸一化，以便在樣本之間進行正確的基因比較。
normalize為了執行歸一化比率方法的中位數，DESeq2 有一個 estimateSizeFactors() 函數可以生成大小因子。我們將在下面的示例中演示此功能，但在典型的 RNA-seq 分析中，此步驟由 DESeq() 函數自動執行，我們將在后面討論。
dds?<-?estimateSizeFactors(dds)通過將結果分配回 dds 對象，我們用適當的信息填充了 DESeqDataSet 對象。我們可以使用以下方法查看每個樣本的歸一化因子：
sizeFactors(dds)現在，要從 dds 中檢索歸一化計數矩陣，我們使用 counts() 函數并添加參數 normalized=TRUE。
normalized_counts?<-?counts(dds,?normalized=TRUE)我們可以將這個歸一化的數據矩陣保存到文件中以備后用：
write.table(normalized_counts,?file="data/normalized_counts.txt",?sep="\t",?quote=F,?col.names=NA)★
注意：DESeq2 實際上并不使用歸一化計數，而是使用原始計數并對廣義線性模型 (GLM) 中的歸一化進行建模。這些歸一化計數對于結果的下游可視化很有用，但不能用作 DESeq2 或任何其他使用負二項式模型執行差異表達分析的工具的輸入。
”歡迎Star -> 學習目錄
更多教程 -> 學習目錄
本文由 mdnice 多平臺發布
總結
以上是生活随笔為你收集整理的RNA-seq 详细教程：搞定count归一化（5）的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 网络云存储技术Windows serve
- 下一篇： 个人成长：简单写写《乔布斯传》读后感