本稿では、マッピング・カウント数の算出を終えたRNA-Seq (シングルセルではなくバルク) データから、Pythonを用いて発現量の正規化を行う方法を解説する。

発現量の正規化方法 FPKMとTPM

参考: https://bi.biopapyrus.jp/rnaseq/analysis/normalizaiton/tpm.html (2022年6月8日アクセス)
遺伝子発現量の正規化手法として、下記の代表的なふたつの手法が知られている。

FPKM/RPKM (Fragments/Reads Per Kilobase of exon per Million mapped fragments/reads, 以降はFPKMと表記する)
TPM (Transcripts Per Million)

最近の論文では TPMを用いた算出を目にすることが多い気がする。実際にFPKMでは転写産物の発現量を正確に表していないことが報告されている（Li et al. Bioinformatics 2010; Wagner et al. Theory Biosci. 2012）。
詳細は上記の論文を参照していただきたいが簡単に説明すると、
FPKMの場合は、
(1) サンプルの総リード数を100万に揃える
(2) 各遺伝子の長さによる補正 (長い遺伝子の方が、多くのリードがマッピングされるから)
の順番で計算が進む。
一方TPMの場合は、
(1) 各遺伝子の長さによる補正
(2) サンプルの総リード数を100万に揃える
の順番で計算が進む。
つまりTPMの場合は、最終的に全ての遺伝子のTPMを合算すると、どのサンプルにおいてもその値が一定になるが、FPKMの場合は若干のずれが生じる。

FPKM/TPMの計算

FPKM

\displaylines{
C_i: 遺伝子 (転写産物) i のリードカウント数 \\
L_i: 遺伝子 (転写産物) i の長さ \\
FPKM_i = \biggl(\frac{C_i}{\sum{C_i}}\biggl)10^6\biggl(\frac{10^3}{L_i}\biggl)
}

左の2つの項で総リード数を100万 (10⁶)に揃え、一番右の項で遺伝子 (転写産物) の長さにより補正しており、先述の順番に計算が行われている (1000 bpあたりの長さなので、10³をかけている)。

TPM

\displaylines{
C_i: 遺伝子 (転写産物) i のリードカウント数 \\
L_i: 遺伝子 (転写産物) i の長さ \\
TPM_i = \biggl(\frac{C_i}{L_i}\biggl)\frac{10^6}{\sum{\frac{C_i}{L_i}}}
}

TPMの場合は、FPKMの場合とは逆で、左の項で、まずそれぞれの遺伝子 (転写産物) のリードカウント数をその長さで補正した後に、補正したリード数の合計が100万 (10⁶)になるように補正している。

Pythonでの実装

今回は、featureCountsで算出したリードカウントデータを、dataframeとして読み込んで計算を行なっている。
データ例

Gene_ID	Chr	Gene_Length (Li)	Gene_Count(Ci)
Gene_A	chr1	12567	3352
Gene_B	chr2	567	32
Gene_C	chr3	123	2

Gene_Length: 遺伝子長データがあるカラムのラベル
Gene_Count: カウントデータがあるカラムのラベル

TPM_FPKM.py

import pandas as pd

# カウント数÷遺伝子長の計算
def count_per_length(df):
    return float(df['Gene_Count'])/float(df['Gene_Length'])

# FPKMの計算
def fpkm(df):
    return (float(df['Counts_per_Length'])/float(count_sum)) * 1000000000

#TPMの計算
def tpm(df):
    return (float(df['Counts_per_Length'])/float(cpt_sum)) * 1000000

# メインスクリプト
data_path = 'data_path' #featureCountの出力ファイルのパス
data = pd.read_table(data_path, header=1) #featureCountファイルの読み込み
count_sum = data['Gene_Count'].sum() #リード数の合計を算出
print('Total count reads', count_sum)

data['Counts_per_Length'] = data.apply(count_per_length, axis=1) #各遺伝子・転写産物における、カウント数÷長さ を算出
counts_per_length_sum = data['Counts_per_Length'].sum() #カウント÷長さ、の合計を算出
print('Sum of counts per length >', counts_per_length_sum)

data['FPKM'] = data.apply(fpkm, axis=1)
data['TPM'] = data.apply(tpm, axis=1)
data.to_excel('Result_File.xlsx')

おわりに

上記のスクリプトにおいて、'Gene_Length', 'Gene_Count', 'data_path'の部分は、使用するデータによって変える必要がある点に注意してほしい。また今回はPythonを用いて算出を行ったが、すべてのサンプルにおいて手作業で行うことを厭わなければ、表計算ソフト (エクセルなど) を用いて計算を行うことも可能である。