兩周前，我接觸了一個RNA-seq的項目，做完之後，我重新思考了FPKM和RPKM的計算，覺得它們很可能是不對的，後來查閱了一些文獻終於驗證了我的想法。現在我重新將這個過程記錄下來：
1. FPKM和RPKM分別是什麼
RPKM是Reads Per Kilobase per Million的縮寫，它的計算方程非常簡單：
RPKM = (10 ^ 6 × n_r) / (L × N)
其中，n_r 是比對至目標基因的read數量；L 是目標基因的外顯子長度之和除以1000，因此，要注意這裡的L單位是kb，不是bp；N 是總有效比對至基因組的read數量。
FPKM是Fregments Per Kilobase per Million的縮寫，它的計算與RPKM極為類似，如下：
FPKM = (10 ^ 6 × n_f) / (L × N)
其中，n_f是比對至目標基因的fregment數量。唯一的區別為：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被轉錄了mRNA_g次，同時樣本X中所有基因的轉錄總次數假定是mRNA_total, 那么正確描述基因g轉錄豐度的值應該是：
r_g = mRNA_g / mRNA_total
同時，樣本X中其他基因的轉錄豐度的計算也和以上式子類似，除了要把分子mRNA_g換為其他基因對應的轉錄次數之外，分母mRNA_total都一樣。於是有趣的事情就是，所有基因轉錄本豐度的均值mean(r_g)將是一個恆定不變的數，由以上定義這個數就是：
r_mean = 1/g_total * sum(r_g) = 1/g_total * [sum(mRNA_g)/mRNA_total]
而
sum(mRNA_g) = mRNA_total
所以
r_mean = 1 / g_total
這個值是由基因的總數所決定的，也就是說，對於同一個物種，不管它的樣本是哪種組織（正常的或病變的等），也不管有多少個不同的樣本，只要它們都擁有相同數量的基因，那麼它們的r_mean都將是一致的。這是一個在進行比較分析的時候，非常有意義的恆等關係。
但在實際的操作中，我們是難以直接計算這些r值的。好在只要能夠保證模型的自洽性，我們是能通過自建一些統計量來對r值進行間接描述的，比如FPKM和RPKM，本質上它們的目的就是為了描述r。雖然如此，但我們也要注意，所有這些要用來描述轉錄本豐度的統計量，都應該能等價描述這一恆等關係。也就是說，不管我們使用了什麼統計量，它所描述出來的轉錄本豐度應該得是其真實豐度(r_g)的m倍（m必須是一個根據模型定出的不變值），它的均值也將是r_mean的m倍，至少這樣才是得到有意義的結果的前提！
（那麼）現在，我們回過頭來看看FPKM和RPKM的計算式，就會發現它們根本做不到。先舉個例子來說明（以FPKM的計算為例），我們假定有兩個來自同一個個體不同組織的樣本X和Y，這個個體只有5個基因，分別為A、B、C、D和E它們的長度分別如下：
由此，我們可以得到，樣本X和Y的轉錄本的不變量：r_mean值都是r_mean = 1 / 5 = 0.2。如果FPKM或RPKM是一個合適的統計量的話，那麼至少，樣本X和Y的平均FPKM（或RPKM）應該相等。
我們以FPKM的計算的為例子，以下這個表格列出的分別是樣本X和Y在這5個基因分別比對上的fregment數和各自總的fregment數量：
於是，按照以上公式我們可以得到樣本X和Y在這5個基因上的FPKM值分別為：
接下來我們就可以計算FPKM的均值了。我們得到，樣本X在這5個基因上的FPKM均值FPKM_mean = 5,680；而樣本Y的FPKM均值卻是FPKM_mean = 161,840!!它們根本不同，而且差距相當大，那麼究竟為什麼會有如此之大的差異？難道是因為我故意構造出來的例子所造成的嗎？當然不是，這是由其數學計算所決定的。
首先，我們可以把FPKM的計算式拆分成兩部分：等價（其實嚴格來講也沒那麼等價）描述某個基因轉錄本數量的統計量（n_f / L） 和 測序獲得的總有效Fregment數量的百萬分之一（N / 10 ^ 6）；看，FPKM便是這兩部分的商！分開來看它們貌似都有點道理，但是合起來的時候其實很沒邏輯，尤其是第二部分（N/10^6），本來既然式子的第一部分 #Fregments mapped to a gene/L)是為了描述某個基因的轉錄本數量，那么正常來講，第二部分就應該是要描述樣本總體的轉錄本數量（或至少是其等價描述）才能說得通，而且可以看得出FPKM(RPMK)是有此意的，因為這本身就是這一統計量的目的。然，它失敗了！N/10^6的大小其實是由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，它的計算是：
TPM = (n_r * read_l / g_l) * 10^6 / T = n_r * read_l * 10 ^ 6 / (g_l * T)
T = sum(n_r * read_l / g_l)
其中，read_l是比對至基因g的平均read長度，g_l 是基因g的外顯子長度之和（這裡無需將其除以1000了）。在不考慮比對的剪切的情況下，read_l這個值往往都是一個固定值（如100bp或者150bp等），因此我們也可以將read_l統一約掉，那麼分子就會蛻變成RPKM計算式的第一部分，但留著會更合理。這樣子，整個統計量就很好理解了，分子是某個基因g的轉錄本數（的等價描述），分母則為樣本中總轉錄本的數量，兩者的比值TPM——便是基因g的轉錄本豐度！並且，簡單計算我們就可以知道TPM的均值是一個獨立於樣本之外的恆定值：
TPM_mean = 10^6 / N
這個值也剛好是r_mean的10^6倍，滿足上述等價描述的關係。我們仍然通過上面的例子來進作說明，為簡單起見我們只把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的計算為例）RPKM = [T * 10^3 / (N * read_l) ] * TPM ，而且在同一個樣本內部由於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）的統計量在描述轉錄本豐度的時候都應該被摒棄。
Tips：由於微信不支持數學公式的顯示，大家可以選擇登陸泛基因的網站（或點擊閱讀原文）查看展示得更漂亮的數學公式。
參考文獻：
http://bib.oxfordjournals.org/content/early/2012/09/15/bib.bbs046.full
http://www.ncbi.nlm.nih.gov/pubmed/22872506
更多基因組學乾貨敬請關注「泛基因」或訪問網站 www.fungenomics.com 並歡迎留言。