RNA-seq的標(biāo)準(zhǔn)化方法

對(duì)于RNA-seq而言，由于 技術(shù)誤差， 測(cè)序深度不同， 基因長度不同，為了能夠比較不同的樣本，比較不同的基因的表達(dá)量，以及使表達(dá)水品分布符合統(tǒng)計(jì)方法的基本假設(shè)，就需要對(duì)原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化。

對(duì)于一個(gè)新興的領(lǐng)域，通常會(huì)有50多種算法，但是最后常用的，其實(shí)也就那么幾個(gè)。在RNA-seq標(biāo)準(zhǔn)化這個(gè)領(lǐng)域也是如此，目前用的最多也就是, RPKM/FPKM, TPM,但是注意，有些時(shí)候一個(gè)方法出現(xiàn)的多，單純是因?yàn)楣緵]有修改他們的分析流程。

為了方便理解，假設(shè)目前你在一次測(cè)序中(即剔除批次效應(yīng))檢測(cè)了一個(gè)物種的3個(gè)樣本，A,B,C，這個(gè)物種有三個(gè)基因G1,G2,G3, 基因長度分別為100, 500, 1000. 通過前期數(shù)據(jù)預(yù)處理，你得到了尚未標(biāo)準(zhǔn)化的表達(dá)量矩陣，如下所示。

事先聲明，以下數(shù)字完全是我瞎編，方便用于后續(xù)的計(jì)算，如有雷同，留個(gè)微信，線下交友。

基因表達(dá)量矩陣

先說三個(gè)簡單的策略，也就是最容易想到的方法

• Total Count, TC, 每個(gè)基因計(jì)數(shù)除以總比對(duì)數(shù), 即文庫大小, 然后乘以不同樣本的總比對(duì)數(shù)的均值
• Upper Quartile, UQ, 和TC方法相似, 即用上四分位數(shù)替代總比對(duì)數(shù)
• Median, Med, 和TC方法相似, 用中位數(shù)代替總比對(duì)數(shù)

上面方法都相似，考慮到我的例子只有三個(gè)基因，所以只展示TC方法的結(jié)果. 可以發(fā)現(xiàn)，原本比其他組觀測(cè)值的A-G2，目前反而是最高。

如果省去TC中的 "乘以不同樣本的總比對(duì)數(shù)的均值" 這一步，那么差不多就是CPM (counts per million)的策略，也就是根據(jù)直接根據(jù)深度對(duì)每個(gè)樣本單獨(dú)進(jìn)行標(biāo)準(zhǔn)化. 在edgeR和limma/voom里面都有出現(xiàn)過。

CPM

TMM（trimmed mean of M value)方法出現(xiàn)在2010年，比TC、 UQ、Med, CPM方法高級(jí)一點(diǎn)，基本假設(shè)是絕大數(shù)的基因不是差異表達(dá)基因.計(jì)算方法有點(diǎn)復(fù)雜，簡單的說就是移除一定百分比的數(shù)據(jù)后，計(jì)算平均值作為縮放因子，對(duì)樣本進(jìn)行標(biāo)準(zhǔn)化。這次我們用R/edgeR來算. 和之前不同，A組的G2基因標(biāo)準(zhǔn)化后還是最低，這就是trim所引起。

> library(edgeR)
> expr <- matrix(c(300,400,500,700,750,650,1000,1300,1800),nrow = 3,byrow = TRUE)
> f <- calcNormFactors(expr, method = "TMM")
> mt_norm <- t(t(mt) /f )
> mt_norm
          [,1]      [,2]      [,3]
[1,]  303.0164  402.5300  491.9114
[2,]  707.0382  754.7438  787.0583
[3,] 1010.0545 1308.2226 1770.8811

DESeq2/DESeq有自己專門的計(jì)算縮放因子(scaling factor)的策略，它的基本假設(shè)就是絕大部分的基因表達(dá)在處理前后不會(huì)有顯著性差異，表達(dá)量應(yīng)該相似，據(jù)此計(jì)算每個(gè)基因在所有樣本中的幾何平均值(geometri mean), 每個(gè)樣本的各個(gè)基因和對(duì)應(yīng)的幾何平均數(shù)的比值的中位數(shù)就是縮放因子(scaling factor). 這里僅僅提到思想，不做計(jì)算。

上述方法都是對(duì)樣本整體進(jìn)行標(biāo)準(zhǔn)化，標(biāo)準(zhǔn)化的結(jié)果只能比較不同樣本之間的同一個(gè)基因的表達(dá)水平。如果要同時(shí)比較不同樣本不同基因之間的表達(dá)量差異，就得考慮到每個(gè)基因的轉(zhuǎn)錄本長度未必相同，畢竟轉(zhuǎn)錄本越長，打算成片段后被觀察到的概率會(huì)高一點(diǎn)。最長嘗試解決這個(gè)問題，應(yīng)該是單端測(cè)序時(shí)代的RPKM(雙端則是FPKM), 全稱為Reads Per Kilobase Million, 其中K象征的是轉(zhuǎn)錄本長度, M象征的是測(cè)序量.

RPKM

對(duì)于A-G1而言，他的表達(dá)量就是300*1000 / (2000 x 100) = 1.5, 其中系數(shù)10e6在這里不需要使用。你可以認(rèn)為它是先對(duì)文庫大小進(jìn)行標(biāo)準(zhǔn)化，然后再根據(jù)基因長度標(biāo)準(zhǔn)化

一波操作下來，感覺數(shù)值就有點(diǎn)可比性了。

下面介紹TPM(transcript per million), 計(jì)算公式如下，其中X表示比對(duì)到基因上的read數(shù)，l表示基因的長度。

TPM

如果手動(dòng)計(jì)算的話，那就是先進(jìn)行 長度標(biāo)準(zhǔn)化，

然后進(jìn)行文庫標(biāo)準(zhǔn)化，即每一列各個(gè)數(shù)值除以每一列的和

目前還有一種技術(shù)叫做UMI，可以進(jìn)行轉(zhuǎn)錄本的絕對(duì)定量，一個(gè)轉(zhuǎn)錄本有多少特異的UMI,就一定程度上代表了它的表達(dá)量。

一些看法

RNA-seq數(shù)據(jù)標(biāo)準(zhǔn)化其實(shí)要分為兩種，樣本間標(biāo)準(zhǔn)化和樣本內(nèi)標(biāo)準(zhǔn)化。

對(duì)于差異表達(dá)分析而言，樣本內(nèi)對(duì)不同轉(zhuǎn)錄本長度進(jìn)行標(biāo)準(zhǔn)化 毫無必要，是的，真的是沒啥意義，粗略比較同一個(gè)基因在兩個(gè)樣本間是否有差異，只要處理好測(cè)序深度這個(gè)問題就行，任何要求用RPKM/FPKM而不是raw count作為差異表達(dá)分析的原始數(shù)據(jù)的分析方法都需要被淘汰掉，而沿用這種分析策略的公司都已經(jīng)很久沒有更新自己流程了。

此外，如果你要比較不同樣本內(nèi)的基因表達(dá)情況，那么目前更推薦用TPM。 因?yàn)樵谒挠?jì)算方法中，能更有效的標(biāo)準(zhǔn)化不同轉(zhuǎn)錄本組成上的差異，而不是簡單除以文庫大小。

題外話，我最近看到一篇發(fā)在Nature上ATAC-seq文章，Method部分提到他用RPKM這個(gè)方法對(duì)每個(gè)bin的read count進(jìn)行標(biāo)準(zhǔn)化?？紤]到每個(gè)bin的大小都一樣，我覺得這個(gè)標(biāo)準(zhǔn)化的方法從定義上更接近CPM。

Correcting for gene length is not necessary when comparing changes in gene expression within the same gene across samples, but it is necessary for correctly ranking gene expression levels within the sample to account for the fact that longer genes accumulate more reads.

對(duì)于差異表達(dá)分析而言，標(biāo)準(zhǔn)化不但要考慮測(cè)序深度的問題，還要考慮到某些表達(dá)量超高或者極顯著差異表達(dá)的基因?qū)е耤ount的分布出現(xiàn)偏倚, 推薦用TMM, DESeq方法進(jìn)行標(biāo)準(zhǔn)化。

參考文獻(xiàn)

• TMM: A scaling normalization method for differential expression analysis of RNA-seq data
• TPM: RNA-Seq gene expression estimation with read mapping uncertainty
• RNA-seq必看綜述: A survey of best practices for RNA-seq data analysis
• Estimating number of transcripts from RNA-Seq measurements (and why I believe in paywall)
• What the FPKM? A review of RNA-Seq expression units