學(xué)習(xí)目標(biāo)

1. 了解計(jì)數(shù)數(shù)據(jù)變換方法的重要性
2. 了解 PCA (principal component analysis)
3. 了解如何使用 PCA 和層次聚類(lèi)評(píng)估樣本質(zhì)量

1. 質(zhì)控

DESeq2 工作流程的下一步是 QC，其中包括樣本和基因程度上，以對(duì)計(jì)數(shù)數(shù)據(jù)執(zhí)行 QC 檢查，以幫助我們確保樣本或重復(fù)看起來(lái)良好。

QC

2. 樣本QC

RNA-seq 分析中一個(gè)有用的初始步驟通常是評(píng)估樣本之間的整體相似性：

1. 哪些樣本彼此相似，哪些不同？
2. 這是否符合實(shí)驗(yàn)設(shè)計(jì)的預(yù)期？
3. 數(shù)據(jù)集中的主要變異來(lái)源是什么？

為了探索樣本的相似性，我們將使用主成分分析 (PCA) 和層次聚類(lèi)方法執(zhí)行樣本級(jí) QC。這些方法或工具使我們能夠檢查重復(fù)彼此之間的相似程度（聚類(lèi)），并確保實(shí)驗(yàn)條件是數(shù)據(jù)變化的主要來(lái)源。樣品級(jí) QC 還可以幫助識(shí)別任何表現(xiàn)出異常值的樣品；我們可以進(jìn)一步探索任何潛在的異常值，以確定是否需要在 DE 分析之前將其刪除。

Sample-level QC

這些無(wú)監(jiān)督聚類(lèi)方法使用 log2 變換的歸一化計(jì)數(shù)運(yùn)行。log2 轉(zhuǎn)換改進(jìn)了可視化的距離。我們將不使用普通的 log2 變換，而是使用正則化對(duì)數(shù)變換 (rlog)，以避免因大量低計(jì)數(shù)基因而產(chǎn)生的任何偏差；

transformation

• 為什么需要進(jìn)行數(shù)據(jù)轉(zhuǎn)換？

許多用于多維數(shù)據(jù)探索性分析的常用統(tǒng)計(jì)方法，尤其是聚類(lèi)和排序方法（例如，主成分分析等），最適合（至少近似地）同方差數(shù)據(jù)；這意味著可觀察量的方差（即，這里是基因的表達(dá)值）不依賴(lài)于均值。然而，在 RNA-seq 數(shù)據(jù)中，方差隨平均值增加。例如，如果直接對(duì)歸一化讀取計(jì)數(shù)矩陣執(zhí)行 PCA，則結(jié)果通常僅取決于少數(shù)高表達(dá)的基因，因?yàn)樗鼈冊(cè)跇颖局g顯示出最大的絕對(duì)差異。避免這種情況的一種簡(jiǎn)單且經(jīng)常使用的策略是取歸一化計(jì)數(shù)值的對(duì)數(shù)加上一個(gè)小的偽計(jì)數(shù)；然而，現(xiàn)在具有低計(jì)數(shù)的基因往往主導(dǎo)結(jié)果，因?yàn)橛捎谛∮?jì)數(shù)值固有的強(qiáng)泊松噪聲，它們?cè)跇颖局g顯示出最強(qiáng)的相對(duì)差異。

• 使用 rlog的優(yōu)勢(shì)：

因此，DESeq2 提供了正則化對(duì)數(shù)轉(zhuǎn)換，或簡(jiǎn)稱(chēng)為 rlog。對(duì)于計(jì)數(shù)高的基因，rlog 轉(zhuǎn)換與普通的 log2 轉(zhuǎn)換差別不大。然而，對(duì)于計(jì)數(shù)較低的基因，這些值會(huì)縮小到所有樣本中基因的平均值。這樣做是為了使 rlog 轉(zhuǎn)換后的數(shù)據(jù)近似同方差。

DESeq2 建議大型數(shù)據(jù)集（100 個(gè)樣本）使用方差穩(wěn)定變換 (vst) 而不是 rlog 來(lái)進(jìn)行計(jì)數(shù)變換，因?yàn)?rlog 函數(shù)可能需要運(yùn)行很長(zhǎng)時(shí)間，而 vst() 函數(shù)在類(lèi)似情況下更快。

3. PCA

主成分分析 (PCA) 是一種用于強(qiáng)調(diào)變化并在數(shù)據(jù)集中降維的技術(shù)。這是一種非常重要的技術(shù)，用于質(zhì)量控制和 Bulk RNA-seq 和單細(xì)胞 RNA-seq 數(shù)據(jù)的分析。

3.1. PCA plots

本質(zhì)上，如果兩個(gè)樣本的基因表達(dá)水平相似，這些基因?qū)o定 PC（主成分）表示的變異有顯著貢獻(xiàn)，則它們將在表示該 PC 的軸上靠近繪制。因此，我們期望生物重復(fù)具有相似的分?jǐn)?shù)（因?yàn)槲覀兊钠谕窍嗤幕蛘诎l(fā)生變化）并聚集在一起。通過(guò)可視化一些示例 PCA 圖最容易理解這一點(diǎn)。

我們?cè)谙旅嬗幸粋€(gè)示例數(shù)據(jù)集和一些相關(guān)的 PCA 圖，以了解如何解釋它們。實(shí)驗(yàn)的元數(shù)據(jù)如下所示。感興趣的主要條件是處理。

dataset

在 PC1 和 PC2 上進(jìn)行可視化時(shí)，我們沒(méi)有看到樣本因處理而分開(kāi)，因此我們決定探索數(shù)據(jù)中存在的其他變異來(lái)源。我們希望我們已經(jīng)在我們的元數(shù)據(jù)表中包含了所有可能的已知變異源，并且我們可以使用這些因素來(lái)為 PCA 圖著色。

PCA_1

我們從cage因子開(kāi)始，但cage因子似乎無(wú)法解釋 PC1 或 PC2 上的變化。

cage

然后，我們按 sex 因素著色，這似乎在 PC2 上分離樣本。這是需要注意的好信息，因?yàn)槲覀兛梢栽谙掠问褂盟鼇?lái)解釋模型中由于 sex 引起的變化并將其回歸。

sex

接下來(lái)我們探索 strain 因子，發(fā)現(xiàn)它可以解釋 PC1 上的變化。

strain

很高興我們能夠確定 PC1 和 PC2 的變異來(lái)源。通過(guò)在我們的模型中考慮它，我們應(yīng)該能夠檢測(cè)到更多因處理而差異表達(dá)的基因。

令人擔(dān)憂的是，我們看到兩個(gè)樣本沒(méi)有與正確的 strain 聚類(lèi)。這表明可能存在樣本交換，應(yīng)進(jìn)行調(diào)查以確定這些樣本是否確實(shí)是標(biāo)記的 strain。如果我們發(fā)現(xiàn)有一個(gè)交換，我們可以交換元數(shù)據(jù)中的樣本。但是，如果我們認(rèn)為它們被正確標(biāo)記或不確定，我們可以從數(shù)據(jù)集中刪除樣本。

我們?nèi)匀粵](méi)有發(fā)現(xiàn)處理是否是 strain 和 sex 后變異的主要來(lái)源。因此，我們探索 PC3 和 PC4 以查看處理是否正在推動(dòng)由這兩個(gè) PC 中的任何一個(gè)代表的變異。

treatment

我們發(fā)現(xiàn)樣本在 PC3 上通過(guò)處理分離，并且對(duì)我們的 DE 分析持樂(lè)觀態(tài)度，因?yàn)槲覀兏信d趣的條件，處理，在 PC3 上分離，我們可以回歸驅(qū)動(dòng) PC1 和 PC2 的變化。

根據(jù)前幾個(gè)主要成分解釋了多少變化，您可能想要探索更多（即考慮更多成分并繪制成對(duì)組合）。即使您的樣品沒(méi)有通過(guò)實(shí)驗(yàn)變量清楚地分離，您仍然可以從 DE 分析中獲得生物學(xué)相關(guān)的結(jié)果。如果您期望效果大小非常小，那么信號(hào)可能會(huì)被無(wú)關(guān)的變異源淹沒(méi)。在您可以識(shí)別這些來(lái)源的情況下，在您的模型中考慮這些來(lái)源很重要，因?yàn)樗鼮闄z測(cè) DE 基因的工具提供了更多功能。

4. 層次聚類(lèi)

與 PCA 類(lèi)似，層次聚類(lèi)是另一種互補(bǔ)的方法，用于識(shí)別數(shù)據(jù)集中的模式和潛在異常值。熱圖顯示數(shù)據(jù)集中所有成對(duì)樣本組合的基因表達(dá)相關(guān)性。由于大多數(shù)基因沒(méi)有差異表達(dá)，樣本之間通常具有很高的相關(guān)性（值高于 0.80）。低于 0.80 的樣本可能表示您的數(shù)據(jù)和/或樣本污染中存在異常值。

沿軸的分層樹(shù)指示哪些樣本彼此更相似，即聚集在一起。頂部的色塊表示數(shù)據(jù)中的子結(jié)構(gòu)，您會(huì)希望看到您的重復(fù)一起作為每個(gè)樣本組的一個(gè)塊。我們的期望是樣本聚集在一起類(lèi)似于我們?cè)?PCA 圖中觀察到的分組。

在下圖中， Wt_3 和 KO_3 樣本沒(méi)有與其他重復(fù)聚類(lèi)在一起。我們想要探索 PCA 以查看我們是否看到相同的樣本聚類(lèi)。

Hierarchical Clustering Heatmap

5. Mov10 QC

現(xiàn)在我們已經(jīng)很好地理解了通常用于 RNA-seq 的 QC 步驟，讓我們?yōu)?Mov10數(shù)據(jù)集進(jìn)行 QC。

5.1. 數(shù)據(jù)轉(zhuǎn)換

• 轉(zhuǎn)換 MOV10 數(shù)據(jù)集的歸一化計(jì)數(shù)

為了促進(jìn) PCA 和層次聚類(lèi)可視化方法的距離或聚類(lèi)，我們需要通過(guò)對(duì)歸一化計(jì)數(shù)應(yīng)用 rlog 變換來(lái)調(diào)節(jié)均值的方差。

歸一化計(jì)數(shù)的 rlog 轉(zhuǎn)換僅在該質(zhì)量評(píng)估期間對(duì)于這些可視化方法是必需的。我們不會(huì)使用這些轉(zhuǎn)換后的計(jì)數(shù)來(lái)確定差異表達(dá)。

# rlog 轉(zhuǎn)換
rld <- rlog(dds, blind=TRUE)

blind=TRUE 參數(shù)是為了確保 rlog() 函數(shù)不考慮我們的樣本組——即以無(wú)偏見(jiàn)的方式進(jìn)行轉(zhuǎn)換。在執(zhí)行質(zhì)量評(píng)估時(shí)，包含此選項(xiàng)很重要。

rlog() 函數(shù)返回一個(gè) DESeqTransform 對(duì)象，這是另一種特定于 DESeq 的對(duì)象。您不只是獲得轉(zhuǎn)換值矩陣的原因是因?yàn)橛糜谟?jì)算 rlog 轉(zhuǎn)換的所有參數(shù)（即大小因子）都存儲(chǔ)在該對(duì)象中。我們使用此對(duì)象繪制 PCA 和層次聚類(lèi)圖以進(jìn)行質(zhì)量評(píng)估。

5.2. PCA

我們現(xiàn)在已準(zhǔn)備好進(jìn)行 QC 步驟，讓我們從 PCA 開(kāi)始吧！

DESeq2 有一個(gè)內(nèi)置函數(shù)，可以在后臺(tái)使用 ggplot2生成 PCA 圖。這很棒，因?yàn)樗刮覀儾槐剌斎氪a行，也不必?cái)[弄不同的 ggplot2 層。此外，它直接將 rlog 對(duì)象作為輸入，從而省去了我們從中提取相關(guān)信息的麻煩。

plotPCA() 需要兩個(gè)參數(shù)作為輸入：DESeqTransform 對(duì)象和 intgroup，即元數(shù)據(jù)中包含有關(guān)實(shí)驗(yàn)樣本組信息列的名稱(chēng)。

# Plot PCA 
plotPCA(rld, intgroup="sampletype")

PCA

默認(rèn)情況下，plotPCA()使用前 500 個(gè)最易變的基因。您可以通過(guò)添加 ntop= 參數(shù)并指定您希望函數(shù)考慮的基因數(shù)量來(lái)更改此設(shè)置。

plotPCA() 函數(shù)將只返回 PC1 和 PC2 的值。如果您想探索數(shù)據(jù)中的其他 PC，或者如果您想識(shí)別對(duì) PC 貢獻(xiàn)最大的基因，您可以使用 prcomp() 函數(shù)。例如，要繪制任何 PC，我們可以運(yùn)行以下代碼：

 # Input is a matrix of log transformed values
 rld <- rlog(dds, blind=T)
 rld_mat <- assay(rld)
 pca <- prcomp(t(rld_mat))

 # Create data frame with metadata and PC3 and PC4 values for input to ggplot
 df <- cbind(meta, pca$x)
 ggplot(df) + geom_point(aes(x=PC3, y=PC4, color = sampletype))

5.3. Hierarchical Clustering

• MOV10 數(shù)據(jù)集層次聚類(lèi)

DESeq2中沒(méi)有內(nèi)置函數(shù)來(lái)繪制熱圖來(lái)顯示所有樣本之間的成對(duì)相關(guān)性和層次聚類(lèi)信息；我們將使用 pheatmap 包中的 pheatmap() 函數(shù)。此函數(shù)不能使用 DESeqTransform 對(duì)象作為輸入，但需要矩陣或數(shù)據(jù)框。因此，要做的第一件事是使用名為 assay() 的函數(shù)，從 rld 對(duì)象檢索該信息，該函數(shù)將 DESeqTransform 對(duì)象中的數(shù)據(jù)轉(zhuǎn)換為簡(jiǎn)單的二維數(shù)據(jù)結(jié)構(gòu)。

# Extract the rlog matrix from the object
rld_mat <- assay(rld)    

接下來(lái)，我們需要計(jì)算所有樣本的成對(duì)相關(guān)值。我們可以使用 cor() 函數(shù)來(lái)做到這一點(diǎn)：

# Compute pairwise correlation values
rld_cor <- cor(rld_mat) 

讓我們看一下相關(guān)矩陣的列名和行名。

head(rld_cor)  

head(meta)  

您會(huì)注意到它們與我們?cè)陂_(kāi)始時(shí)使用的元數(shù)據(jù)數(shù)據(jù)框中為樣本提供的名稱(chēng)相匹配。這很重要，因此我們可以使用下面的注釋參數(shù)在頂部繪制一個(gè)色塊。此塊可輕松實(shí)現(xiàn)層次聚類(lèi)的可視化。

# Load pheatmap package
library(pheatmap)

# Plot heatmap using the correlation matrix and the metadata object
pheatmap(rld_cor, annotation = meta)

當(dāng)您使用 pheatmap() 進(jìn)行繪圖時(shí)，層次聚類(lèi)信息用于將相似的樣本放在一起，并且該信息由沿軸的樹(shù)結(jié)構(gòu)表示。注釋參數(shù)接受一個(gè)數(shù)據(jù)框作為輸入，在我們的例子中它是元數(shù)據(jù)框。

pheatmap

總體而言，我們觀察到高相關(guān)性 (> 0.999)，表明沒(méi)有異常樣本。此外，與 PCA 圖類(lèi)似，您會(huì)看到樣本按樣本組聚集在一起?？傊@些圖向我們表明數(shù)據(jù)質(zhì)量很好，我們有信心可以進(jìn)行差異表達(dá)分析。

歡迎Star -> 學(xué)習(xí)目錄

更多教程 -> 學(xué)習(xí)目錄

本文由mdnice多平臺(tái)發(fā)布