1. FPKM 和 RPKM 分別是什麼
RPKM 是 Reads Per Kilobase per Million 的縮寫,它的計算方程非常簡單:
其中,nr是比對至目標基因的 read 數量;L 是目標基因的外顯子長度之和除以 1000,因此,要注意這裡的 L 單位是 kb,不是 bp;N 是總有效比對至基因組的 read 數量。
FPKM 是 Fregments Per Kilobase per Million 的縮寫,它的計算與 RPKM 極為類似,如下:
其中,nf是比對至目標基因的 fregment 數量。FPKM 與 RPKM 的唯一的區別為:F 是 fragments,R 是 reads,如果是 pair-end 測序,每個 fragments 會有兩個 reads,FPKM 只計算兩個 reads 能比對到同一個轉錄本的 fragments 數量,而 RPKM 計算的是可以比對到轉錄本的 reads 數量(不管 pair-end 的兩個 reads 是否能比對到同一個轉錄本上)。如果是 single-end 測序,那麼 FPKM 和 RPKM 計算的結果將是一致的。
以上是這兩個量的計算方式,它們這樣計算的目都是為了解決在計算 RNA-seq 轉錄本豐度的兩個 bias:
(1)相同表達豐度的轉錄本,往往會由於其基因長度上的差異,導致測序獲得的 Read(Fregment)數不同。總的來說,越長的轉錄本,測得的 Read(Fregment)數越多。
(2)由測序文庫的不同大小而引來的差異。即同一個轉錄本,其測序深度越深,通過測序獲得的 Read(Fregment)數就越多。
FPKM 和 RPKM 通過同時除以 L(轉錄本長度)和 N(有效比對的 Read(Fregment)總數)的辦法,最終將不同樣本(或者同個樣本在不同條件下)的轉錄本豐度歸一化到一個能夠進行量化比較的標準上。上面的式子看起來似乎合情合理,但是它們卻都做錯了。
2. 為什麼 FPKM/RPKM 是錯的
要回答這個問題,我們需要先撇開所有形式上的計算,重新思考到底什麼是 RNA 轉錄本的表達豐度這個問題。事實上,對於任何一個製備好的樣本,它上面任何一個基因的表達量(或者說豐度),都將已是一個客觀存在的值,這個值是不管你改變了多少測序環境都不會變的。而且總共有多少個不同的基因在表達,實際上也已經是客觀定好了的。但我們開始以這樣一種「先知」的形式來理解的時候,有趣的事情就開始出現了。
此刻,我們可以假定,對於樣本 X,其有一個基因 g 被轉錄了 mRNAg次,同時樣本 X 中所有基因的轉錄總次數假定是 mRNAtotal, 那么正確描述基因 g 轉錄豐度的值應該是:
同時,樣本 X 中其他基因的轉錄豐度的計算也和以上式子類似,除了要把分子mRNAg換為其他基因對應的轉錄次數之外,分母mRNAtotal都一樣。於是有趣的事情就是,所有基因轉錄本豐度的均值 rmean將是一個恆定不變的數,由以上定義這個數就是:
而
所以
這個值是由基因的總數所決定的,也就是說,對於同一個物種,不管它的樣本是哪種組織(正常的或病變的等),也不管有多少個不同的樣本,只要它們都擁有相同數量的基因,那麼它們的rmean都將是一致的。(這是一個)這是一個在進行比較分析的時候,非常有意義的恆等關係。
但在實際的操作中,我們是難以直接計算這些 r 值的。好在只要能夠保證模型的自洽性,我們是能通過自建一些統計量來對 r 值進行間接描述的,比如 FPKM 和 RPKM,本質上它們的目的就是為了描述 r。雖然如此,但我們也要注意,所有這些要用來描述轉錄本豐度的統計量,都應該能等價描述這一恆等關係。也就是說,不管我們使用了什麼統計量,它所描述出來的轉錄本豐度應該得是其真實豐度 rg的 m 倍(m 必須是一個根據模型定出的不變值),它的均值也將是rmean的 m 倍,至少這樣才是得到有意義的結果的前提!
(那麼)現在,我們回過頭來看看 FPKM 和 RPKM 的計算式,就會發現它們根本做不到。先舉個例子來說明(以 FPKM 的計算為例),我們假定有兩個來自同一個個體不同組織的樣本 X 和 Y,這個個體只有 5 個基因,分別為 A、B、C、D 和 E 它們的長度分別如下:
由此,我們可以得到,樣本 X 和 Y 的轉錄本的不變量:rmean值都是 rmean=1/5=0.2。如果 FPKM 或 RPKM 是一個合適的統計量的話,那麼至少,樣本 X 和 Y 的平均 FPKM(或 RPKM)應該相等。
我們以 FPKM 的計算的為例子,以下這個表格列出的分別是樣本 X 和 Y 在這 5 個基因分別比對上的 fregment 數和各自總的 fregment 數量:
於是,按照以上公式我們可以得到樣本 X 和 Y 在這 5 個基因上的 FPKM 值分別為:
接下來我們就可以計算 FPKM 的均值了。我們得到,樣本 X 在這 5 個基因上的 FPKM 均值 FPKMmean=5,680;而樣本 Y 的 FPKM 均值卻是 FPKMmean=161,840!! 它們根本不同,而且差距相當大,那麼究竟為什麼會有如此之大的差異?難道是因為我故意構造出來的例子所造成的嗎?當然不是,這是由其數學計算所決定的。
首先,我們可以把 FPKM 的計算式拆分成兩部分:等價(其實嚴格來講也沒那麼等價)描述某個基因轉錄本數量的統計量 nf/L 和測序獲得的總有效 Fregment 數量的百萬分之一 N/106;看,FPKM 便是這兩部分的商!分開來看它們貌似都有點道理,但是合起來的時候其實很沒邏輯,尤其是第二部分 N/106,本來既然式子的第一部分(nf/L)是為了描述某個基因的轉錄本數量,那么正常來講,第二部分就應該是要描述樣本總體的轉錄本數量(或至少是其等價描述)才能說得通,而且可以看得出 FPKM(RPMK) 是有此意的,因為這本身就是這一統計量的目的。然,它失敗了!N/106的大小其實是由 RNA-seq 測序深度所決定,並且是一個和總轉錄本數量無直接線性關係的統計量——N 與總轉錄本數量之間的關係還受轉錄本的長度分布所決定,而這個分布往往在不同樣本中是有差異的!比如,雖然有些基因,有效比對到它們的 Fregment 數是相等的,但總的來講長度越長的基因,其被轉錄的次數就越少。也就是說,N 必須將各個被轉錄的基因的長度考慮進去才能正確描述總體的轉錄本數!而 FPKM(RPKM)顯然沒有做到這一點,這便是 FPKM(RPKM)出錯的內在原因。
3. 那麼應該是用什麼樣統計量才合適
其實,通過以上分析,我們已經可以確定一個更加合理的統計量來描述 RNA 轉錄本的豐度了。這個統計量在 12 年所發表的一篇關於討論 RPKM 的文章(RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ.Theory Biosci. 2012.)中已被提到過了,它稱之為 TPM —— Transcripts Per Million,它的計算是:
其中,readl 是比對至基因 g 的平均 read 長度,gl 是基因 g 的外顯子長度之和(這裡無需將其除以 1000 了)。在不考慮比對的剪切的情況下,readl 這個值往往都是一個固定值(如 100bp 或者 150bp 等),因此我們也可以將 readl 統一約掉,那麼分子就會蛻變成 RPKM 計算式的第一部分,但留著會更合理。這樣子,整個統計量就很好理解了,分子是某個基因 g 的轉錄本數(的等價描述),分母則為樣本中總轉錄本的數量,兩者的比值 TPM——便是基因 g 的轉錄本豐度!並且,簡單計算我們就可以知道 TPM 的均值是一個獨立於樣本之外的恆定值:
這個值也剛好是 rmean的 106倍,滿足上述等價描述的關係。我們仍然通過上面的例子來進作說明,為簡單起見我們只把 fregment 換為 read,其他數字都一樣,並且統一假設 read_l 都是一樣的:
接著,我們可以分別計算樣本 X 和 Y 的 TPM_mean, 並且很明顯它們都是 200000 = 10^6 / 5. 而且,經過這樣的標準化之後,X 和 Y 就處於同樣的一個標準上了,此刻,彼此之間的比較分析才是真正有意義的。
4. 既然 FPKM/RPKM 是錯的,那為什麼大家還在用,而且還真找到了(能被實驗所驗證)有價值的結果呢?
關於對於這個問題,我也思考過。而且我們都知道 2008 那篇關於 RPKM 的文章用實驗結果證明了,RPKM 是一個合適的統計量,符合 qPCR 的驗證結果。但歸根到底,眼見未必為實,實驗是表象的,我們更應該從其本質意義和原理上去考慮。FPKM/RPKM 之所以看起來會是一個合適的值,我想主要原因有二:
其一,它們和 TPM 之間存在一定的正比關係。這在通過它們各自的數學計算式就可以看出來(以 RPKM 的計算為例):
而且在同一個樣本內部由於 T,N 和 read_l 實際上都是定值,因此同個樣本內的 RPKM 和 TPM 是可以恆等轉換的,只是在樣本與樣本之間就不行,因為它們之間的轉換因子大小不一了!如以上例子,對於樣本 X,TPM 轉換到 RPKM 的轉換因子為:0.0284,但樣本 Y 的轉換因子為:0.8092。而由於這個標準的改變,會導致其原本所要描述的「轉錄本豐度」變得不可比較。然,這其實不是最根本的原因,更本質的原因是,這個轉換會對本來已經正確標準化了的結果,再次做了一次無意義的不等變換,最終導致了結果不可解釋。如何理解呢,後文會有補充,這裡先簡單說一下:這個數學轉換式子僅是告訴了我們這樣子來計算是可行的,但是在 RNA-seq 的實際應用場景中,它其實是無生物(或物理)意義的;
其二,實驗驗證的精度是有限的,常用的 qPCR 也只能給出定性的比較結果,而且實驗驗證也未必總能成功。
5. 總結
現在回過頭來總結一下。事實上,FPKM/RPKM 最大的問題就在於其無意義性。我們所要表達出來的任何統計量,它的變化都應該要能對應到其物理或生物過程中的變化,如果做不到這一點,那麼這個統計量往往都是無意義的,用它得到的結果就算看起來符合預期也只不過是數值上的巧合,本質上是不可解釋的。FPKM/RPKM 的分母 (N/10^6) 並不具有任何形式的生物意義,它所能表達出來的這個量,只能代表測序深度的變化,而無法作為表達生物過程的量,比如無法代表(等價代表)樣本中轉錄本的總量。
一個統計量該如何計算,說到底都只是一個「術」的問題,而我們應該儘可能在接近其本質意義的地方去確定。
FPKM/RPKM 和 TPM 存在一定的正比關係,因此我們在使用 FPKM/RPK 時,有些時候確實也能獲得可以被實驗所驗證的「好」結果,但其實它是一個橡皮筋,它的單位刻度是會隨著樣本的不同而改變的。到頭來,樣本之間的差異比較實際上也只是在不同的標準下進行的,這樣的比較就算得到了所謂的「好」結果,那又有什麼意義,根本就是個錯誤的東東。想想就是由於這種統計量,我們一定已經獲得了許多的假陽性結果,同時也肯定錯過了許多本來真正有意義的差異,真是彎路走盡也不知,而且還浪費了大堆的心情和時間。
這篇文章:A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Briefings in Bioinformatics.10.1093/bib/bbs046. 對 7 種主要的 RNA-seq 標準化方法(但不包含本文提到的 TPM)做了一個詳細的比較,它用實際結果進行比較(不同於本文所用的數學方式)也得出了 RPKM 這個統計量是應該被摒棄的結論,因為它所描述出來的結果是最不合理的,其實所有類似於 RPKM(包括 FPKM)的統計量在描述轉錄本豐度的時候都應該被摒棄。
http://bib.oxfordjournals.org/content/early/2012/09/15/bib.bbs046.full
http://www.ncbi.nlm.nih.gov/pubmed/22872506
來源:泛基因