RNA-seq中的那些統(tǒng)計學問題（二）FPKM/RPKM之外的那些標準化方法

微笑如酒 2021-11-14

展開全文

1. 標準化

1.1. House-keeping gene(s)

1.2. spike-in

1.3. CPM

1.4. TCS

1.5. Quantile

1.6. Median of Ration

1.7. TMM

2. 為什么說FPKM和RPKM都錯了？

2.1. FPKM和RPKM分別是什么

2.2. 什么樣才算好的統(tǒng)計量

2.3. FPKM和RPKM犯的錯

2.4. TPM是一個合適的選擇

1. 標準化

由于不同文庫測序深度不同，比較前當然要進行均一化！用總reads進行均一化可能最簡單，其基于以下兩個基本假設(shè)：

絕大多數(shù)的gene表達量不變；
高表達量的gene表達量不發(fā)生改變；

但在轉(zhuǎn)錄組中，通常一小部分極高豐度基因往往會貢獻很多reads，如果這些“位高權(quán)重”的基因還是差異表達的，則會影響所有其它基因分配到的reads數(shù)，而且，兩個樣本總mRNA量完全相同的前提假設(shè)也過于理想了。那如何比較呢，各個方家使出渾身解數(shù)，有用中位數(shù)的，有用75分位數(shù)的，有用幾何平均數(shù)的，有用TMM(trimmed mean of Mvalues)的等等，總之要找一個更穩(wěn)定的參考值。

1.1. House-keeping gene(s)

矯正的思路很簡單，就是在變化的樣本中尋找不變的量

那么在不同RNA-seq樣本中，那些是不變的量呢？一個很容易想到的就是管家基因 (House-keeping gene(s))

那么 Human 常用的 House-keeping gene 怎么確定？

目前大家用的比較多的一個human housekeeping gene list 來源于下面這篇文章，是2013年發(fā)表在 Cell系列的 Trends in Genetics 部分的一篇文章

1.2. spike-in

使用Housekeeping gene的辦法來進行相對定量，這種辦法在一定程度上能夠解決我們遇到的問題。但其實這種辦法有一個非常強的先驗假設(shè)：housekeeping gene的表達量不怎么發(fā)生變化。其實housekeeping gene list有幾千個，這幾千個基因有一定程度上的變化是有可能的

spike-in方法：在RNA-Seq建庫的過程中摻入一些預(yù)先知道序列信息以及序列絕對數(shù)量的內(nèi)參。這樣在進行RNA-Seq測序的時候就可以通過不同樣本之間內(nèi)參（spike-in）的量來做一條標準曲線，就可以非常準確地對不同樣本之間的表達量進行矯正

比較常用的spike-in類型：ERCC Control RNA

ERCC = External RNA Controls Consortium

ERCC就是一個專門為了定制一套spike-in RNA而成立的組織，這個組織早在2003年的時候就已經(jīng)宣告成立。主要的工作就是設(shè)計了一套非常好用的spike-in RNA，方便microarray，以及RNA-Seq進行內(nèi)參定量

1.3. CPM

CPM(count-per-million)

$\text{cpm}_i=\frac{\text{read count of Gene}_i}{\text{total reads}/10^6}$

1.4. TCS (Total Count Scaling)

簡單來說，就是找出多個樣本中l(wèi)ibrary size為中位數(shù)的樣本，作為參考樣本，將所有的樣本的library size按比例縮放到參考樣本的水平

選擇一個library size為中位數(shù)的sample，以它為baseline，計算出其它sample對于baseline的normalization factor，即一個縮放因子：

$d_i=\frac{S_{baseline}}{S_i}$

然后基于該縮放因子對特定的sample中的每個基因的read count進行標準化（縮放）：

$x_i'=d_ix_i$

1.5. Quantile

簡單來說，就是排序后求平均，然后再回序

在R里面，推薦用preprocessCore 包來做quantile normalization，不需要自己造輪子啦！
但是需要明白什么時候該用quantile normalization，什么時候不應(yīng)該用，就復(fù)雜很多了

1.6. Median of Ratio (DESeq2)

該方法基于的假設(shè)是，即使處在不同條件下的不同個樣本，大多數(shù)基因的表達是不存在差異的，實際存在差異的基因只占很小的部分那么我們只需要將這些穩(wěn)定的部分找出來，作為標準化的內(nèi)參，依據(jù)內(nèi)參算出各個樣本的標準化因子

（1）對每個基因計算幾何平均數(shù)，得到一個假設(shè)的參考樣本(pseudo-reference sample)

gene	sampleA	sampleB	pseudo-reference sample
EF2A	1489	906	sqrt(1489 * 906) = 1161.5
ABCD1	22	13	sqrt(22 * 13) = 17.7
…	…	…	…

（2）對每個樣本的每個基因?qū)τ趨⒖紭颖居嬎鉌old Change

gene	sampleA	sampleB	pseudo-reference sample	ratio of sampleA/ref	ratio of sampleB/ref
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
…	…	…	…	…	…

（3）獲取每個樣本中Fold Change的中位數(shù)，我們就得到了非DE基因代表的Fold Change，該基因就是我們選擇的該樣本的內(nèi)參基因，它的Fold Change就是該樣本的標準化因子

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))

1.7. TMM (Trimmed Mean of M value, edgeR)

該方法的思想與DESeq2的Median of Ratio相同，假設(shè)前提都是：大多數(shù)基因的表達是不存在差異的

它與DESeq2的不同之處在于對內(nèi)參的選擇上：

DESeq2選擇一個內(nèi)參基因，它的Ratio/Fold-Change就是標準化因子

edgeR選擇一組內(nèi)參基因集合，它們對標準化因子的計算均有貢獻：加權(quán)平均

（1）移除所有未表達基因

（2）從眾多樣本中找出一個數(shù)據(jù)趨勢較為平均的樣本作為參考樣本

對所有樣本求總Read數(shù)；

各樣本除以各自的總Read數(shù)，得到修正Read數(shù)；

求出各自樣本修正Read數(shù)的Q3值（第3個四分位數(shù)）；

所有的Q3值求平均，與平均Q3相差最小的樣本即是參考樣本。

（3）找出每個樣本中的代表基因集，參考這些代表基因集的fold change，計算出該樣本的標準化因子

尋找樣本的代表基因集：依據(jù)基因的偏倚程度和Reads數(shù)大小選出——偏倚程度小、reads數(shù)居中的基因

衡量偏倚度的量：LFC (log fold change)

$LFC=\log_2\frac{\text{sample}_i}{ref}$

LFC過大或過小都表示具有偏倚性，LFC越大表示reads數(shù)在sample_i中越高，即偏向sample_i；LFC越小表示reads數(shù)在ref中越高，即偏向ref

衡量reads數(shù)的量：read的幾何平均數(shù) (read geometric mean, RGM)

RGM越大表示基因reads越多，RGM越小表示基因reads越少

結(jié)合偏倚度、read數(shù)畫出散點圖：

偏倚度小、表達量居中的那些基因落在圖中的紅線附近

由參考代表基因集計算樣本的標準化因子：

對這些代表基因集計算加權(quán)平均數(shù)：

$\frac{\sum_i^n LFC_i\times \text{ReadCount}_i}{\sum_i^n \text{ReadCount}_i}$

該加權(quán)平均數(shù)就能代表該樣本的標準化因子，只是經(jīng)過了log變換，所以需要恢復(fù)為它的正值：

$\text{ScalingFactor}=2^{\text{weight-average}}$

2. 為什么說FPKM和RPKM都錯了？

2.1. FPKM和RPKM分別是什么

FPKM和RPKM分別是什么

RPKM是Reads Per Kilobase per Million的縮寫，它的計算方程非常簡單：

$RPKM=\frac{10^6 \times n_r}{L \times N}$

FPKM是Fregments Per Kilobase per Million的縮寫，它的計算與RPKM極為類似，如下：

$FPKM=\frac{10^6 \times n_f}{L \times N}$

與RPKM唯一的區(qū)別為：F是fragments，R是reads，如果是PE(Pair-end)測序，每個fragments會有兩個reads，F(xiàn)PKM只計算兩個reads能比對到同一個轉(zhuǎn)錄本的fragments數(shù)量，而RPKM計算的是可以比對到轉(zhuǎn)錄本的reads數(shù)量而不管PE的兩個reads是否能比對到同一個轉(zhuǎn)錄本上。如果是SE(single-end)測序，那么FPKM和RPKM計算的結(jié)果將是一致的。

這兩個量的計算方式的目的是為了解決計算RNA-seq轉(zhuǎn)錄本豐度時的兩個bias：

相同表達豐度的轉(zhuǎn)錄本，往往會由于其基因長度上的差異，導(dǎo)致測序獲得的Read（Fregment）數(shù)不同?？偟膩碚f，越長的轉(zhuǎn)錄本，測得的Read（Fregment）數(shù)越多；
由測序文庫的不同大小而引來的差異。即同一個轉(zhuǎn)錄本，其測序深度越深，通過測序獲得的Read（Fregment）數(shù)就越多。

2.2. 什么樣才算好的統(tǒng)計量

首先，到底什么是RNA轉(zhuǎn)錄本的表達豐度這個問題

對于樣本X，其有一個基因g被轉(zhuǎn)錄了mRNA_g次，同時樣本X中所有基因的轉(zhuǎn)錄總次數(shù)假定是mRNA_total, 那么正確描述基因g轉(zhuǎn)錄豐度的值應(yīng)該是：

$r_g=\frac{\text{mRNA}_g}{\text{mRNA}_{total}}$

則一個樣本中基因表達豐度的均值為

$r_{mean}=\frac{1}{g_{total}} \sum_g^G r_g = \frac{1}{g_{total}} \frac{\sum_g^G \text{mRNA}_g}{\text{mRNA}_{total}}$

而

$\sum_g^G \text{mRNA}_g=\text{mRNA}_{total}$

所以

$r_{mean}=\frac{1}{g_{total}}$

這個期望值竟然和測序狀態(tài)無關(guān)！僅僅由樣本中基因的總數(shù)所決定的

也就是說，對于同一個物種，不管它的樣本是哪種組織（正常的或病變的），也不管有多少個不同的樣本，只要它們都擁有相同數(shù)量的基因，那么它們的r_mean都將是一致的

由于上面的結(jié)果是在理論情況下推導(dǎo)出來的，實際上我們無法直接計算這個r，那么我們可以嘗試通過其他方法來近似估計r，只要這些近似統(tǒng)計量可以隱式地包含這一恒等關(guān)系即可

2.3. FPKM和RPKM犯的錯

實際數(shù)據(jù)來證明

假定有兩個來自同一個個體不同組織的樣本X和Y，這個個體只有5個基因，分別為A、B、C、D和E它們的長度分別如下：

$r_{mean}$ 值是:

$r_{mean}=\frac 15=0.2$

如果FPKM或RPKM是一個合適的統(tǒng)計量的話，那么，樣本X和Y的平均FPKM（或RPKM）應(yīng)該相等。

以下這個表格列出的分別是樣本X和Y在這5個基因分別比對上的fregment數(shù)和各自總的fregment數(shù)量：

可以算出：樣本X在這5個基因上的FPKM均值FPKM_mean = 5,680;樣本Y在這5個基因上的FPKM均值FPKM_mean = 161,840

很明顯，它們根本不同，而且差距相當大

究竟為什么會有如此之大的差異？

可以從其公式上找到答案

首先，我們可以把FPKM的計算式拆分成如下兩部分。

第一部分的統(tǒng)計量是對一個基因轉(zhuǎn)錄本數(shù)量的一個等價描述（雖然嚴格來講也沒那么等價）：

$\frac{n_f}{L}$

第二部分的統(tǒng)計量是測序獲得的總有效Fregment數(shù)量的百萬分之一：

$\frac{N}{10^6}$

這么一拆，就可以看出這個公式的問題了：邏輯上根本說不通嘛！

尤其是第二部分（N/10⁶），本來式子的第一部分是為了描述一個基因的轉(zhuǎn)錄本數(shù)量，那么正常來講，第二部分就應(yīng)該是樣本的轉(zhuǎn)錄本總數(shù)量（或至少是其總數(shù)量的等價描述）才能形成合理的比例關(guān)系，而且可以看出來FPKM/RPMK是有此意的，這本來就是這個統(tǒng)計量的目的。

但是N/10⁶并不能描述樣本的轉(zhuǎn)錄本總數(shù)量。N/10⁶的大小其實是由RNA-seq測序深度所決定的，并且是一個和總轉(zhuǎn)錄本數(shù)量無直接線性關(guān)系的統(tǒng)計量——N與總轉(zhuǎn)錄本數(shù)量之間的關(guān)系還受轉(zhuǎn)錄本的長度分布所決定，而這個分布往往在不同樣本中是有差異的！

2.4. TPM是一個合適的選擇

這個統(tǒng)計量在2012年所發(fā)表的一篇討論RPKM的文章（RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ. Theory Biosci. 2012.）中就被提出來了，稱之為TPM —— Transcripts Per Million，它的計算是：

$TPM=\frac{\frac{n_r \times read_l}{g_l}}{T}=\frac{n_r \times read_l \times 10^6}{g_l \times T}$

$T=\sum_{g=i}^G (\frac{n_r \times read_l}{g_l})_i$

簡單計算之后我們就可以發(fā)現(xiàn)TPM的均值是一個獨立于樣本之外的恒定值，它等于：

$TPM_{mean}=\frac{10^6}{N}$

這個值剛好是r_mean的一百萬倍，滿足等價描述的關(guān)系。

參考資料：

(1) 孟浩巍《生物信息學100個基礎(chǔ)問題 —— 第38題當轉(zhuǎn)錄組普遍變化時RNA-Seq怎么進行分析(1)？》

(2) 孟浩巍《生物信息學100個基礎(chǔ)問題 —— 第38題當轉(zhuǎn)錄組普遍變化時RNA-Seq怎么進行分析(2)？》

(3) 【生信菜鳥團】quantile normalization到底對數(shù)據(jù)做了什么？

(4) Introduction to DGE

(5) 生信菜鳥團：StatQuest生物統(tǒng)計學專題 - library normalization進階之edgeR的標準化方法

(6) 【簡書】為什么說FPKM和RPKM都錯了？