學習目標

1. 瞭解如何在歸一化過程中列出不同的 uninteresting factors(無關因素)
2. 瞭解常用的歸一化方法，已經如何使用
3. 瞭解如何建立 DESeqDataSet 物件及其結構
4. 瞭解如何使用 DESeq2 進行歸一化

1. 歸一化

差異表達分析工作流程的第一步是計數歸一化，這是對樣本之間的基因表達進行準確比較所必需的。

每個基因的對映讀數計數是 RNA 表達以及許多其他因素的結果。歸一化是調整原始計數值以解決“無關”因素的過程。以這種方式，表達水平在樣本之間或樣本內更具可比性。

在歸一化過程中經常考慮的“無關”因素：

1.1. 測序深度

考慮測序深度對於比較樣本之間的基因表達是必要的。在下面的示例中，每個基因在樣本 A 中的表達似乎是樣本 B 的兩倍。然而，這是樣本 A 的測序深度加倍的結果。

1.2. 基因長度

計算基因長度對於比較同一樣本中不同基因之間的表達是必要的。在下面的示例中，基因 X 和基因 Y 具有相似的表達水平，但對映到基因 X 的讀數數量將比對映到基因 Y 的讀數多得多，因為基因 X 更長。

1.3. RNA組成

樣本之間一些高度差異表達的基因、樣本之間表達的基因數量的差異或汙染的存在可能會扭曲某些型別的歸一化方法。建議考慮 RNA 組成以準確比較樣本之間的表達，這在進行差異表達分析時尤為重要。

在下面的示例中，假設樣本 A 和樣本 B 之間的測序深度相似，並且除了基因差異表達之外的每個基因在樣本之間呈現相似的表達水平。樣本 B 中的計數會受到 差異表達基因的極大影響，它佔據了大部分計數。因此，樣本 B 的其他基因的表達似乎低於樣本 A 中的相同基因。

歸一化不僅對於差異表達分析必不可少，對於探索資料分析、資料視覺化以及探索或比較樣本之間或樣本內的計數也是必要的。

2. 歸一化方法

幾種常見的歸一化方法：

• RPKM/FPKM：不推薦用於樣本間比較

雖然 TPM 和 RPKM/FPKM 歸一化方法都考慮了測序深度和基因長度，但不推薦使用 RPKM/FPKM。原因是RPKM/FPKM方法輸出的歸一化計數值在樣本之間沒有可比性。

使用 RPKM/FPKM 歸一化，每個樣本的 RPKM/FPKM 歸一化計數總數會有所不同。因此，您不能在樣本之間平均比較每個基因的歸一化計數。

RPKM-歸一化計數表：

例如，在上表中，樣本 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 使用比率中值方法。在使用者端只有一個步驟，但在後端涉及多個步驟，如下所述。

1. 建立一個偽參考樣本（逐行幾何平均值）

對於每個基因，都會建立一個偽參考樣本，該樣本等於所有樣本的幾何平均值。

2. 計算每個樣本與參考的比率

對於每個樣本中的每個基因，計算比率（樣本/參考）（如下所示）。由於大多數基因沒有差異表達，因此每個樣本中的大多數基因在樣本中的比例應該相似。

3. 計算每個樣本的歸一化因子（大小因子）

給定樣本的所有比率的中值（上表中的列）被視為該樣本的歸一化因子（大小因子），計算如下。請注意，差異表達的基因不應影響中值：

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 軸是頻率）。

比率中位數法假設並非所有基因都差異表達；因此，歸一化因子應考慮樣本的測序深度和 RNA 組成（大的離群基因不會影響中值比率值）。該方法對上調/下調和大量差異表達基因的不平衡具有魯棒性。

4. 使用歸一化因子計算歸一化計數值

這是通過將給定樣本中的每個原始計數值除以該樣本的歸一化因子來執行的，生成歸一化計數值。這是針對所有計數值（每個樣本中的每個基因）執行的。例如，如果樣本 A 的中值比率為 1.3，樣本 B 的中值比率為 0.77，則可以按如下方式計算歸一化計數：

• Raw Counts

• Normalized Counts

請注意，歸一化計數值不是整數。

以上步驟僅作為演示，在實際使用DESeq2過程中，只需要一步命令，即可完成計算。

3. Mov10 歸一化

現在我們知道了計數歸一化的理論，我們將使用 DESeq2 對 Mov10 資料集的計數進行歸一化。這需要幾個步驟：

1. 確保 metadata 資料框的行名存在，並且與 counts 資料框的列名順序相同。
2. 建立一個 DESeqDataSet 物件
3. 生成歸一化 counts

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

下一步是對計數資料進行歸一化，以便在樣本之間進行正確的基因比較。

為了執行歸一化比率方法的中位數，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多平臺釋出