聊聊轉錄組測序——2.數據分析與解讀（下）

前言

nn

上次我們聊到獲取read比對信息（BAM文件），下面就是喜聞樂見的，通過BAM文件獲取基因表達量和差異表達的計算。

表達量計算

nn

基因表達量最直接的分析手段就是計算比對到每個基因的reads有多少條，在轉錄組測序中，我們通常稱這個數字為count。前文說過，使用基因組比對的方法，我們獲得了每條read比對到基因組上的位置，因此，我們需要把基因組的位置對應到基因。這個提供基因位置的注釋文件通常以GTF或GFF3格式呈現，GTF格式示例如下：

第一列是染色體編號，第三列是本行的特徵（feature），如gene、transcript、exon、CDS等（實際上大多數情況下，計算表達量只要帶exon的行就夠了），第四列和第五列是基因組起始和終止位置，第七列是正負鏈，第九列是注釋信息（可以包括類似基因ID、轉錄本ID、基因名等信息）。詳細的說明見GTF2.2: A Gene Annotation Format。

有GTF文件後，就可以利用注釋信息計算每個基因/轉錄本/外顯子比對了多少reads，從而獲取counts值。值得一提的是，為提高比對的準確性，HISAT2和STAR等軟體在比對的過程中就已經結合GTF文件中提供的轉錄本剪接信息進行了優化。

nn

然而，由於每個測序樣品的起始RNA量不同，文庫量不同，測序數據量不同… 原始的counts值好比沒做內參的qPCR，不適合直接作為表達量用於樣品間比較。因此，我們需要對counts值進行校正，或均一化。

Counts值的校正有幾種常見的演算法，最早流行的是RPKM（Reads Per Kilobase per Million mapped reads）或在此基礎上衍生的FPKM（FragmentsnPer Kilobase per Million mapped reads）。這種校正方法是用基因的Counts值除以基因長度，再除以測序數據量（實際使用的是比對成功的reads總數），從而得到每個基因相對的表達水平，見http://www.arrayserver.com/wiki/index.php?title=RPKM（或者直接百度百科，難得比WIKI寫得好）。公式如下（公式來源What the FPKM? A review of RNA-Seq expression units）：

Xi為Counts值，li為轉錄本長度，N為比對成功的reads總數，FKPM和RPKM本質上差別不大

此外，近幾年TPM（Transcripts PernMillion）也越來越常用，這種演算法計算的是某個轉錄本在一個樣品測到的所有轉錄本中所佔的比例。因此使用基因長度校正Counts值後，還要除以所有基因的 Counts值/基因長度n的和，從而獲取這個比例（雖然這個描述有點繞，可以慢慢理解，理解不了也沒關係）公式如下（公式來源What the FPKM? A review of RNA-Seq expression units）：

Xi為Counts值，li為轉錄本長度

第三類常見的校正方法（其實是很多種不同的校正方法）是基於大多數表達量中等的基因在樣品間不存在表達差異的假設，使用TMM（trimmed mean of M values）或類似方法針對每個樣本計算一個校正係數，再用校正係數直接生成Normalized counts值，具體的演算法就不展開講了（例如下面的公式是DESeq的校正係數的演算法，kij是基因i在樣品j中的counts值，m是樣本數，Sj是樣品j的校正係數）。由於這類方法通常是差異表達計算軟體生成的（edgeR、DESeq等），並且符合負二相分布/泊松分布，因此非常適合直接用於後續分析。

nn

此外還有中位值（使用每個樣品中表達量位於中位數的基因校正）、四分位校正（選取表達量位於25~75%的基因的總counts值校正）等等，由於用的人相對較少，這裡就不一一介紹了。

我個人比較推崇DESeq等軟體校正出來的NormalizednCounts，一方面是評價比較高，一方面用來算差異表達的p值也更方便。但Normalized Counts最大的問題在於，每次加入新樣本時都需要重新校正，對於多批次完成的大樣本量研究、或資料庫的維護非常麻煩。因此，作為TCGA數據首選的校正方式，TPM同樣是一個不錯的選擇。此外，人工合成的Spike-in同樣可以作為外參校正表達量，不過這就是另一個話題了。

獲取表達量之後，就可以做測序文章的第一張圖——Heatmap（熱圖）了。

Heatmap通過顏色的深淺變化，直觀體現每個基因在不同樣品中的表達情況（如上圖，紅色越深表達量越高，藍色越深表達量越低），並且看到不同樣品/分組之間表達譜的差異。當然，為了圖片好看，通常會對表達量取對數（由於部分基因Count等於0，所以取對數前可以給所有基因的表達量加1，或加一個比較小的數值，如0.001），並使用每個基因在所有樣品中表達量的均值對基因表達量進行均一化。當然，取對數和均一化都是可選步驟，作圖時需根據不同情況調整。

另外，Heatmap的上方和左側進化樹一樣的結構，是聚類（Clustering）的結果，表達趨勢越接近的兩個基因/樣品，在聚類結果中也更為接近。通常認為能聚到一起的基因簇/樣品簇有著更為接近的生物學特徵。

#常用軟體#

Counts值計算常用HTSeq和featureCounts，此外部分軟體自帶counts值計算，如RSEM、Salmon等。

TPM和RPKM用RSEM都能算，或者其實直接寫個代碼手算都可以。

TMM之類的校正有不少R包可以用，我一般用DESeq（DESeq1和DESeq2沒區別）來計算，edgeR也可以。

nn

Heatmap同樣可以用R包畫，pheatmap應該是裡面最簡單的，通常三到四行代碼就能畫一張最簡單的帶有聚類結果的heatmap。如果不擅長R的話，可以用Cluster和Treeview的組合。

差異表達分析

使用基因表達量對基因的差異表達水平進行計算。如果樣品多，非要對每個基因用T檢驗也不是不行，但是這種方法不是轉錄組測序的主流做法（雖然某些研究認為測序獲取的基因表達量符合正態分布，但不是主流觀點，而樣本量太小時秩和檢驗的檢驗功效偏低，並且從應用的角度講，每組不足三個重複樣品，也沒法做直接T檢驗或秩和檢驗）。

通常認為基因的counts值服從某種統計學的分布模型（負二相分布、泊松分布等），然後利用該分布模型分析每個基因的差異表達情況，並計算組間差異。計算差異後，軟體會輸出每個基因對應的Fold-change以及p值。

在大多數轉錄組測序文章中，通常使用q值或FDR來衡量表達差異的顯著水平。這是使用多重檢驗校正控制錯誤率後的結果。舉個例子，當我們使用p<0.05為基準時，那麼每20個基因，就可能有一個單純因為概率被認為p<0.05，而當我們對上萬個基因分析表達差異時，必須排除這種因為概率導致的假陽性，因此需要使用Bonferroni校正等方法控制false discovery rate（FDR），並得到更加嚴格的q值。（當然，無論是p值還是q值，都只是基於統計學的計算，雖然p<0.05是一個約定俗成的閾值，但q值完全沒有必要限制使用0.05的閾值。另外，我們樣品的差異大小與樣本本身性質相關，比如病人組織間的表達差異通常較大，而用比較溫和的藥物刺激細胞系的前後表達差異通常較小，因此在做後續篩選時建議根據實際情況設定閾值。）

得到Fold-change以及p值/q值後，就可以繪製火山圖（VolcanonPlot）和其它散點圖，火山圖中通常橫軸代表Fold change，縱軸代表p值，類似的散點圖只要認清坐標軸就可以知道圖片中描述的是什麼了。

最後，在轉錄組測序中，唯p值/q值論是不可取的，一方面p/q值的可信區間的範圍是隨著|Fold-change|的增大而縮小的，因此只有|Fold-change|足夠大時，才能確保差異分析的假陽性率足夠小。況且在|Fold-change|太小的情況下，即使p/q值顯著，這種差異也未必具有生物學意義。此外在篩選差異表達基因時，每個基因的表達水平也需要考慮，因為丰度越低的基因，受測序誤差的影響也越顯著（同樣是1個count的誤差，對於表達量為1和100的基因，影響顯然是不同的）；雖然有些軟體會低丰度的基因進行校正，但是低丰度基因的p值計算在大多數軟體中都存在較大誤差，況且一個本底表達量比較低的基因，生物學功能可能也不會太強（用qPCR也很難檢測）。因此，建議大家手動去除丰度過低的基因，我個人習慣counts>50（12G數據量），也可以使用RPKM>1等條件。

#常用軟體#

差異表達分析（都是R包）：

nn

DESeq/DESeq2，大概是使用率最高的R包之一，但不支持非整數的counts值；

edgeR：這個有點老了，但還是有人在用，不太推薦；

EBSeq：RSEM使用的分析軟體，可以使用非整數counts值進行分析；

此外還有sSeq、BaySeq等等使用率較高的R包，就不一一介紹了。

Cufflinks自帶計算差異表達的Cuffdiff，可以直接使用FKPM值進行計算，但個人建議如果不是萬不得已，還是不要使用Cuffdiff進行分析。

火山圖和其它散點圖可以用R包ggplot做（其實用Excel都能做），因為我自己不常做這個圖，用R做也不是很麻煩，所以也沒有專門找過可以直接出結果的軟體，大家有推薦的請踴躍留言吧。

推薦閱讀：