全文5,314字，閱讀16分鐘。
這本是三年多之前我發在公眾號上的一篇舊文，一些偶然的機會，發現不少朋友也在討論這個問題，因此我重新做了梳理並發出來。這樣的討論是有益的，放在今天也不過時，我也願意和更多的朋友一起進一步討論這個問題。
曾經（2015年），我接觸了一個RNA-seq的項目，做完之後，我重新思考了FPKM和RPKM的計算，覺得它們很可能是不對的（當時是第一次接觸RNA-seq數據，還不知道TPM的存在），後來查閱了一些文獻終於驗證了我的想法。現在我重新將這個過程記錄下來。
FPKM和RPKM分別是什麼
RPKM是Reads Per Kilobase per Million的縮寫，它的計算方程非常簡單：
其中，
n_r 是比對至目標基因的read數量；
L 是目標基因的外顯子長度之和除以1000。
因此，要注意這裡的L單位是kb，不是bp；N 是總有效比對至基因組的read數量。
FPKM是Fregments Per Kilobase per Million的縮寫，它的計算與RPKM極為類似，如下：
其中，n_f是比對至目標基因的fregment數量。
與RPKM唯一的區別為：F是fragments，R是reads，如果是PE(Pair-end)測序，每個fragments會有兩個reads，FPKM只計算兩個reads能比對到同一個轉錄本的fragments數量，而RPKM計算的是可以比對到轉錄本的reads數量而不管PE的兩個reads是否能比對到同一個轉錄本上。如果是SE(single-end)測序，那麼FPKM和RPKM計算的結果將是一致的。
以上是這兩個量的計算方式，它們這樣計算的目是為了解決計算RNA-seq轉錄本豐度時的兩個bias：
FPKM和RPKM通過同時除以L（轉錄本長度）和N（有效比對的Read（Fregment）總數）的方式，最終將不同樣本（或者同個樣本在不同條件下）的轉錄本豐度歸一化到一個能夠進行量化比較的標準上。
上面的式子看起來似乎合情合理，但是它們卻都做錯了！
為什麼FPKM/RPKM是錯的
要回答這個問題，我們需要先撇開所有形式上的計算，重新思考到底什麼是RNA轉錄本的表達豐度這個問題。
事實上，對於任何一個製備好測序文庫的待測樣本，它上面任何一個基因的表達量（或者說豐度），都將已是一個客觀存在的值，這個值是不管你改變了多少測序環境都不會變的。而且總共有多少個不同的基因在這這一刻進行了表達，實際上也已經是客觀定好了——難道不是嗎？！
當我們開始以這樣一種「先知」的形式來理解問題時，有趣的事情就開始出現了。
此刻，我們可以假定，對於樣本X，其有一個基因g被轉錄了mRNA_g次，同時樣本X中所有基因的轉錄總次數假定是mRNA_total, 那么正確描述基因g轉錄豐度的值應該是：
同時，樣本X中其他基因的轉錄豐度的計算也應該和上述公式類似。除了要把分子mRNA_g換為其他基因對應的轉錄次數之外，分母mRNAtotal都一樣。於是有趣的事情就是，所有基因轉錄本豐度的均值r_mean將是一個恆定不變的數，由以上定義這個數就是：
而
所以
這個期望值竟然和測序狀態無關！僅僅由樣本中基因的總數所決定的。也就是說，對於同一個物種，不管它的樣本是哪種組織（正常的或病變的），也不管有多少個不同的樣本，只要它們都擁有相同數量的基因，那麼它們的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值是:
如果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在這5個基因上的FPKM均值FPKM_mean = 161,840
看到了嗎？它們根本不同，而且差距相當大，那麼究竟為什麼會有如此之大的差異？難道是因為我故意構造出來的例子所造成的嗎？當然不是，我們分析一下就會發現，這是完全是其數學公式所決定的。
首先，我們可以把FPKM的計算式拆分成如下兩部分。
第一部分的統計量是對一個基因轉錄本數量的一個等價描述（雖然嚴格來講也沒那麼等價）:
第二部分的統計量是測序獲得的總有效Fregment數量的百萬分之一:
看到了嗎？FPKM其實就只是這兩部分的商！這有道理嗎？分開來看它們似乎都有點道理，但是合起來的時候其實很沒邏輯。
尤其是第二部分（N/10^6），本來式子的第一部分是為了描述一個基因的轉錄本數量，那么正常來講，第二部分就應該是樣本的轉錄本總數量（或至少是其總數量的等價描述）才能形成合理的比例關係，而且可以看出來FPKM/RPMK是有此意的，這本來就是這個統計量的目的。
可是，它卻失敗了！
N/10^6的大小其實是由RNA-seq測序深度所決定的，並且是一個和總轉錄本數量無直接線性關係的統計量——N與總轉錄本數量之間的關係還受轉錄本的長度分布所決定，而這個分布往往在不同樣本中是有差異的！
比如，在有些基因中，雖然有效比對到基因的Fregment數是相等的，但一般來說長度越長的基因，其被轉錄的次數就越少。那也就是說，N必須將各個被轉錄的基因的長度考慮進去才能正確描述總體的轉錄本數！而FPKM/RPKM顯然沒有做到這一點，這便是FPKM和RPKM出錯的內在原因。
那麼應該是用什麼樣的統計量才合適？
其實，通過以上分析，我們已經可以確定一個更加合理的統計量來描述RNA轉錄本的豐度了。
這個統計量在2012年所發表的一篇討論RPKM的文章（RPKM measure is inconsistent among samples. Wagner GP, Kin K, Lynch VJ. Theory Biosci. 2012.）中就被提出來了，稱之為TPM —— Transcripts Per Million，它的計算是：
其中，read_l是比對至基因g的平均read長度，g_l 是基因g的外顯子長度之和（這裡無需將其除以1000，這是沒必要的）。在不考慮read剪切的情況下，read_l這個值是測序長度，它往往都是一個固定值（如100bp或者150bp等），所以我們也可以將read_l統一消掉，這時分子就會蛻變成RPKM計算式的第一部分，只是留著會更好。
這樣子，整個統計量就很好理解了，分子是某個基因g的轉錄本數（或者說是轉錄本數量的等價描述），分母為樣本中總轉錄本的數量，兩者的比值TPM——便是基因g的轉錄本豐度！並且，簡單計算之後我們就可以發現TPM的均值是一個獨立於樣本之外的恆定值，它等於：
這個值剛好是r_mean的一百萬倍，滿足等價描述的關係。
我們仍然通過上面的例子來作進一步的說明，為簡單起見我只把fregment換為read，其他數字都一樣，並且統一假設read_l都是一樣的：
接著，我們可以分別計算樣本X和Y的TPM_mean,並且很明顯它們都是
而且，經過這樣的標準化之後，X和Y就處於同樣的一個標準上了，此刻，在同一個坐標系下的差異比較才是真正有意義的。
既然FPKM/RPKM是錯的，那為什麼大家還在用，而且還真找到了（能被實驗所驗證）的有價值結果呢？
關於對於這個問題，我也思考過。而且我們都知道2008年那篇提出RPKM的文章也用實驗結果來證明RPKM是一個合適的統計量，符合qPCR的驗證結果。
但歸根到底，眼見未必為實，實驗也是表象的一個結果，我們應該回到其本質意義和原理上去思考。
FPKM/RPKM之所以看起來會是一個合適的值，我想主要原因有三：
其一，它們和TPM之間存在一定的正比關係，如下。通過比較它們各自的數學計算式就可以看出來（以RPKM的計算為例）：
可以看到，在RPKM的公式中，原來TPM是其中的一個因子，而且在同一個樣本內部由於T、N和read_l實際上都是定值——也就是說式子的第一個因子此時變成常數了！這就使得在同個樣本內的RPKM和TPM是可以恆等轉換的，可是這個情況在不同樣本之間就不成立了！
因為不同樣本之間的轉換因子大小是不一樣的！還是上面的例子，對於樣本X，TPM轉換到RPKM的轉換因子為：0.0284，但樣本Y的轉換因子卻是：0.8092。而由於這個標準的改變，會導致其原本所要描述的「轉錄本豐度」變得不可比較。
然，這其實還不是最根本的原因，更本質的原因是，這個轉換會對本來已經正確標準化了的結果，再次做了一次無意義的不等變換，最終導致了結果不可解釋。如何理解呢？後文會有補充，這裡先簡單說一下：這個數學轉換式子僅是告訴了我們這樣子來計算是可行的，但是在RNA-seq的實際應用場景中，它其實是無生物（或物理）意義的；
其二，實驗驗證的精度是有限的，常用的qPCR也只能給出定性的比較結果，而且實驗驗證也未必總能成功；
其三，習慣了，形成了路徑依賴，大家都這樣做，所以也這樣跟隨，這應該是最不可取的一個情況。
總結
現在回過頭來總結一下。事實上，FPKM/RPKM最大的問題就在於其無意義性。我們所要表達出來的任何統計量，它的變化都應該能對應到物理或生物過程中的變化，如果做不到這一點，那麼這個統計量往往都是無意義的，用它得到的結果就算看起來符合預期也只不過是數值上的巧合，本質上是不可解釋的。FPKM/RPKM的分母(N/10^6)並不具有任何形式的生物意義，它所能表達出來的這個量，只能代表測序深度的變化，而無法作為對生物基因表達量的描述，即無法代表（等價代表）樣本中轉錄本的總量。
一個統計量該如何計算，說到底都只是一個「術」的問題，而我們應該儘可能在接近其本質意義的地方去確定。
FPKM/RPKM和TPM存在一定的正比關係，因此我們在使用FPKM/RPKM時，有些時候確實也能獲得可以被實驗所驗證的「好」結果，但其實它是一個橡皮筋，它的單位刻度是會隨著樣本的不同而改變的。到頭來，樣本之間的差異分析實際上也只是在不同的標準下進行的，這樣的比較就算得到了所謂的「好」結果，那又有什麼意義？根本就是個錯誤的東西。
想想就是由於這種統計量，我們一定已經獲得了許多的假陽性結果，同時也肯定錯過了許多本來真正有意義的差異，真是彎路走盡卻不自知，而且過程中還浪費了多少的心情和時間？
這篇文章：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）的統計量在描述轉錄本豐度的時候都應該被摒棄。
時至今日，這個反思的經歷依然可以不斷給我帶來有益的啟發，其中最重要的一點就是絕不盲從、絕不輕信，如果我不能理解其原理和來龍去脈，那麼我寧可不用，這也是為什麼後來我總是更願意自己去把方法實現的一個重要原因。想起一句話：
What I cannot create, I do not understand.
- Richard P.Feynman(理察.菲利普斯.費曼)
費曼
願與你共勉。
參考文獻：
http://bib.oxfordjournals.org/content/early/2012/09/15/bib.bbs046.fullhttp://www.ncbi.nlm.nih.gov/pubmed/22872506
----/ END /----
DESeq2差異基因分析和批次效應移除
高通量數據中批次效應的鑑定和處理 - 系列總結和更新
Nature重磅綜述 |關於RNA-seq，你想知道的都在這
後臺回復「生信寶典福利第一波」或點擊閱讀原文獲取教程合集