【生物信息】RPKM, FPKM和TPM

注：這幾個名詞是RNA-Seq數據分析中的基礎，在此小結一下。

 

在RNA-Seq的分析中，對基因或轉錄本的read counts數目進行標准化（normalization）是一個極其重要的步驟，因為落在一個基因區域內的read counts數目取決於基因長度和測序深度。很容易理解：一個基因越長、測序深度越深，落在（比對到）該基因所在區域的read counts數目就會相對越多。當我們進行基因差異表達的分析時，往往是在多個樣本中比較不同基因的表達量，如果不進行數據標准化，比較結果是沒有意義的（因為差異可能是由於測序深度導致的，而不是本身表達上的差異）。因此，我們需要標准化的兩個關鍵因素就是基因長度和測序深度，常常用RPKM (Reads Per Kilobase Million), FPKM (Fragments Per Kilobase Million) 和 TPM (Transcripts Per Million)作為標准化數值。

 

RPKM和FPKM

---

計算RPKM主要包括以下三步：

1. 計算與測序深度有關的系數：計算每個樣本中reads的總數並除以$10^6$——此時就可以得到"per million"的縮放系數；
2. 去除測序深度的影響：每個reads數除以上面得到的"per million"縮放系數，得到對應基因在每百萬reads中所占比例，即reads per million (RPM)；
3. 去除基因長度的影響：用上一步的結果再除以對應基因的長度（通常是所有外顯子長度的總和，以kb為單位），得到每百萬reads每一千鹼基對中包含的reads數，即RPKM。

計算RPKM的公式可以表示如下（如圖1中的示例）：

$$RPKM(x)=\frac{Reads \ per \ transcript}{\frac{total \ reads}{10^6} \times \frac{transcript \ length}{1000}} = \frac{Reads \ per \ transcript}{million \ reads \times transcript \ length (kb)} = \frac{10^9 \times C_x}{R \times L_x} \quad \cdots (1)$$

其中：

• x可以表示一個基因或轉錄本，或基因組上一段特定的區域（下面統一當做基因）；
• $C_x$表示比對到基因x外顯子區域的reads數；
• R表示當前樣本中包含的全部reads數，$R = \sum_{i=1}^{N} C_i$，其中N表示該樣本中的基因總數；
• $L_x$表示基因x外顯子區域包含的鹼基數（bp）.

圖1：RPKM的計算

FPKM與RPKM的計算過程相同，它們的區別是：RPKM用於單端測序結果，FPKM用於雙端測序結果（如圖2所示）。因為雙端測序中，每一個fragment會包含兩個reads，使用FPKM來計算基因的表達豐度時，可以避免把同一個fragment的兩個reads計算兩次（實際上只需要計算一次）。

圖2：單端測序read與雙端測序read

單端read與雙端read比對到基因組的示意圖如圖3所示：

圖3：單端測序read與雙端測序read比對到參考基因組

 

TPM

---

TPM與RPKM最大的區別在於消除兩種影響的次序：在TPM中先消除基因長度的影響，再消除測序深度的影響。計算TPM的過程也可以分為三個步驟：

1. 將每個read counts除以對應基因的長度（外顯子區域的長度，單位為kb），此時得到每千個鹼基包含的reads數，即（reads per kilobase, RPK）；
2. 將一個樣本中的RPK加起來的總數除以$10^6$，得到"per million"縮放系數（這是兩種方法計算結果不同的主要來源，因為這里的總數是消除了基因長度的影響之后得到的RPK，而不是原始read counts之和）；
3. 用RPK除以"per million"縮放系數，得到TPM。

計算公式表示如下：

$$TPM(x) = \frac{C_x \times r \times 10^6}{L_x \times T} = \frac{C_x / L_x \times 10^6}{\sum_{i=1}^{N} C_i / L_i} \quad \cdots (2)$$

其中：

• x可以表示一個基因或轉錄本，或基因組上一段特定的區域（下面統一當做基因）；
• r表示reads的長度（例如所有reads包含的平均鹼基數），一般為常數；
• $C_x$表示比對到基因x外顯子區域的reads數；
• T表示當前樣本中所有基因除以自身長度（kb）后求和的結果，$T = \sum_{i=1}^{N} r \times C_i / L_i$，其中N表示該樣本中的基因總數；
• $L_x$表示基因x外顯子區域包含的鹼基數（kp）.

因為交換了兩次計算的次序，TPM最終得到的結果中，每個樣本總的TPM值是相同的，這樣的結果便於樣本間差異的比較。

 

RPKM與TPM的比較

---

 

根據提出TPM的論文（Günter P. Wagner et al., 2012）中的介紹：從理論上來說，任意基因的表達量（$mRNA_x$）占樣本中所有基因總的表達量（$mRNA_T$）比值（文中提到的相對摩爾濃度，$rmc$）的均值應該是一個固定的值。用公式可表示如下：

$$rmc_x = \frac{mRNA_x}{mRNA_T}$$

$$E[rmc] = \frac{1}{N} \sum_{i=1}^{N} \frac{mRNA_i}{mRNA_T} = \frac{1}{N} \frac{\sum_{i=1}^{N} mRNA_i}{mRNA_T} = \frac{1}{N}$$

其中：

• x可以表示一個基因或轉錄本，或基因組上一段特定的區域（下面統一當做基因）；
• $mRNA_x$表示基因x的表達量（豐度）；
• $mRNA_T$表示所有基因在該樣本中總的表達量；
• N表示樣本中包含的基因總數；
• $E[rmc]$表示相對摩爾濃度的均值。

這里的均值與下面例子中的均值的含義相同：假設現在有100萬元，20個人分，無論這些人怎么分，最終人均都是5萬元，即$100/20$；人均占比都是$\frac{5}{100} = \frac{1}{20}$。上述相對摩爾濃度的均值就相當於這里的人均占比，只與人數有關。

由於以上的約束條件，我們在構造統計量對任意基因在某個樣本中的相對表達豐度進行估計時，該估計量的均值也應該是一個只與基因總數相關的值。於是就有了對RPKM的反對以及TPM的提出。因為任意一個樣本中所有基因的TPM之和都等於$10^6$，因此其均值都等於$\frac{10^6}{N}$。

 

計算RPKM和TPM的代碼示例

---

 下面展示Python代碼：

import pandas as pd
import numpy as np


def read_counts2tpm(df, sample_name):
    """
    convert read counts to TPM (transcripts per million)
    :param df: a dataFrame contains the result coming from featureCounts
    :param sample_name: a list, all sample names, same as the result of featureCounts
    :return: TPM
    """
    result = df
    sample_reads = result.loc[:, sample_name].copy()
    gene_len = result.loc[:, ['Length']]
    rate = sample_reads.values / gene_len.values
    tpm = rate / np.sum(rate, axis=0).reshape(1, -1) * 1e6
    return pd.DataFrame(data=tpm, columns=sample_name, index=df['Gene'])

def read_counts2rpkm(df, sample_name):
    result = df
    sample_reads = result.loc[:, sample_name].copy()
    gene_len = result.loc[:, ['Length']]
    total_reads = np.sum(sample_reads.values, axis=0).reshape(1, -1)
    rate = sample_reads.values / gene_len.values
    tpm = rate / total_reads * 1e6
    return pd.DataFrame(data=tpm, columns=sample_name, index=df['Gene'])

從代碼上可以看到，計算RPMK與計算TPM的差別有兩點，在RPKM中計算total_reads用的是原始reads數（第23行），且用該值來消除測序深度對結果的影響（第25行）；在TPM中，使用的是消除基因長度的影響之后的值來求和（第16行）。

下面通過示例數據來計算結果：

# raw data
a = pd.DataFrame(data = {
    'Gene': ("A","B","C","D","E"),
    'Length': (100, 50, 25, 5, 1),
     'S1': (80, 10,  6,  3,   1),
     'S2': (20, 20, 10, 50, 400)
})

tpm = read_counts2tpm(a, ['S1', 'S2'])

rpkm = read_counts2rpkm(df=a, sample_name=['S1', 'S2'])

構造的示例數據包括基因名稱、長度、樣本1和樣本2的read counts，該結果可以通過featureCounts或Rsubread等上游工具得到。

注：featureCounts（2014）和Rsubread（2019）由同一個實驗室先后開發。

 

結果如下：

# TPM
Gene    S1              S2	
A	281690.140845	486.618005
B	70422.535211	973.236010
C	84507.042254	973.236010
D	211267.605634	24330.900243
E	352112.676056	973236.009732

# RPKM
Gene	S1	S2
A	8000.0	400.0
B	2000.0	800.0
C	2400.0	800.0
D	6000.0	20000.0
E	10000.0	800000.0

可以驗證得到，在兩個樣本中TPM之和都等於$10^6$，且均值都等於200000。

 

R相關的代碼可參考：

https://gist.github.com/slowkow/6e34ccb4d1311b8fe62e

 

References

---

https://www.plob.org/article/16013.html

https://www.rna-seqblog.com/rpkm-fpkm-and-tpm-clearly-explained/

https://izabelcavassim.wordpress.com/2015/03/09/rpkm-and-fpkm-normalization-units-of-expression/

http://www.mi.fu-berlin.de/wiki/pub/ABI/GenomicsLecture13Materials/rnaseq2.pdf

https://yourgene.pixnet.net/blog/post/68952443-paired-end-v.s.-single-end

https://slideplayer.com/slide/14812922/

https://zhuanlan.zhihu.com/p/55988984

https://link.springer.com/article/10.1007%2Fs12064-012-0162-3 

http://bioinf.wehi.edu.au/featureCounts/

https://gist.github.com/slowkow/6e34ccb4d1311b8fe62e