遺伝子発現量カウントの補正方法(RPM, RPKM, TPM)【Python】

2020-11-29

RNA-Seqにおける遺伝子発現量カウントは遺伝子長やリード数などで補正をする必要があります。ここではその補正方法として、RPM/FPM、RPKM/FPKM、TPMについて説明し、Pythonを用いたプログラムで実装してみます。

OS

Windows10 + WSL (WindowsユーザーのためのPythonを用いたゲノム解析環境)

Python

Python 3.7.9

モジュール

BCBio-GFF 0.6.6
pandas 1.1.3

[toc]

遺伝子発現量カウントの補正の必要性

RNA-SeqではmRNAをランダムに短い断片に切断することで効率よくその断片の塩基配列を同定していきます。その短い断片の塩基配列をリードと呼びます。そのリードをリファレンス配列にマッピングし、目的の遺伝子にマッピングされたリード数で発現の程度を評価するのがRNA-Seqにおける発現量解析の基本的な考え方になります。

ここで重要なのは、1つ1つのmRNAに等しい確率でリードがマッピングされるわけではないということです。長いmRNAでは多くの断片に切断されるので、1つのmRNAを表す配列に多数のリードがマッピングされることになりますが、短いmRNAの場合は切断されてできる断片の数も少ないので、1つのmRNAを表す配列にマッピングされるリードの数も少なくなります。つまり、転写産物を表すmRNAが長ければ長いほどリード数も多くなるという関係があるので、リード数の比較から転写産物の発現量を比べるためには、そのmRNAの長さも補正してそろえてやる必要があります。

また、RNA-Seqのシークエンス深度が深い(=合計のリード数が多い)ほど、1つの転写産物を表すmRNAに多くのリードがマッピングされることになります。合計のリード数が2倍であれば、単純に考えてそれぞれのmRNAにマッピングされるリードの数も2倍になってしまいます。そのため異なるサンプル同士でmRNAにマッピングされたリード数を比較してその発現量を比べるためには、まずそれぞれのリード数の合計をそろえてやる必要があります。

遺伝子発現量カウントを補正するための準備

遺伝子発現量カウントの補正には総リード数とそれぞれの転写産物を表すmRNAの長さの情報が必要になります。総リード数は遺伝子発現量カウントのデータそのものから取得できますが、mRNAの長さの情報はアノテーションファイル(*.gff, *.gtfなど)から取得しましょう。

ここでは、次の記事のようにhtseq-countを用いて遺伝子発現量カウントを行った場合の補正の方法を説明します。

この記事で説明している通り、アノテーションファイルにおけるfeature領域を表すIDが「gene_id」として保存されているものとして説明していきます。mRNAの長さはその遺伝子のエクソンの長さの合計として取得します。次のようにして、アノテーションファイルからgene_idごとのエクソンの長さの合計を「gene_length.tsv」として取得しましょう。

from BCBio import GFF
import pandas as pd
 
# 読み込むアノテーションファイルを指定する
file_name = 'sample.gff'

gene_id_list = []
length_list = []
 
# 指定されたGFFファイルを開く
with open(file_name) as handle:

    last_gene_id = ''
 
    # 親子関係を無視して、1つのアノテーションを1つのrecordとして取得して順番に処理する
    for record in GFF.parse(handle, target_lines=1):
        if record.features[0].type == 'exon':
            gene_id = record.features[0].qualifiers['gene_id'][0]
            length = record.features[0].location.end - record.features[0].location.start
            
            # 今取得したgene_idが1つ前に取得したgene_idと同じ場合は、前のlengthに加える
            if gene_id == last_gene_id:
                length_list[-1] += length
            # 新たなgene_idであれば、新たにgene_idとlengthをリストに加える
            else:
                gene_id_list.append(gene_id)
                length_list.append(length)
                last_gene_id = gene_id

# タブ区切りのテキストファイルで保存する
df = pd.DataFrame({'gene_id':gene_id_list, 'length':length_list})
df.to_csv('gene_length.tsv', sep='\t', index=False)

アノテーションファイルの扱い方については以下の記事もご覧ください。

htseq-count以外の発現量カウントプログラムを用いた場合は出力結果に転写産物の長さも含まれている場合もあるので、その場合はその値を用いてください。

なお、htseq-countを用いると、出力ファイルの一番最後に以下のような情報が付加されますが、この部分は削除しておいてgene_xxxxの部分だけにしておいてください。

主な遺伝子発現量カウントの補正方法

RPM / FPM

サンプル間の総リード数をそろえることで、サンプル同士のリード数の比較を可能にする方法です。具体的には、ゲノムにマッピングされた総リード数を100万としたときのリード数に変換して正規化します。

それでは、htseq-countの出力結果「count_raw.tsv」をRPM(FPM)に補正してみましょう。

import pandas as pd

# サンプルの総リード数を100万とした場合のリード数に変換する
def normalize_totalreads(df):
    return 10**6 * df / df.sum()

df_count_raw = pd.read_table('count_raw.tsv', index_col=0, header=None)  # データの読み込み
df_rpm = normalize_totalreads(df_count_raw)  # RPMに補正
df_rpm.to_csv('rpm.tsv', sep='\t', header=False)  # ファイルに出力

しかし、mRANの長さによってマッピングされるリード数が異なってしまうので、この方法ではその補正はできていません。つまり、同じ発現量の転写産物でもそれを表すmRNAにマッピングされるリード数は異なってしまい、転写産物同士の発現量の比較はできません。

そこで考え出された方法が次のRPKM(FPKM)です。

RPKM / FPKM

RPKM(FPKM)は、RPM(FPM)と同様にまずサンプルの総リード数を100万に正規化して、サンプル間での比較を可能にしたうえで、転写産物のmRNAの長さを1000塩基としたときのリード数に変換した値です。これによって、サンプル間の発現量の比較とともに、転写産物同士の転写産物の比較も可能にすることを目指しています。

それでは、htseq-countの出力結果「count_raw.tsv」をRPKM(FPKM)に補正してみましょう。転写産物の長さは「gene_length.tsv」に用意しておきます。

import pandas as pd

# サンプルの総リード数を100万とした場合のリード数に変換する
def normalize_totalreads(df):
    return 10**6 * df / df.sum()

# 転写産物の長さを1000とした場合のリード数に変換する
def normalize_genelength(df, gene_length):
    return (df.T * 10**3 / gene_length).T

# データの読み込み
df_count_raw = pd.read_table('count_raw.tsv', index_col=0, header=None)  
df_gene_length = pd.read_table('gene_length.tsv', index_col=0)
gene_length = df_gene_length['length']

# RPKMに補正
df_temp = normalize_totalreads(df_count_raw)
df_rpkm = normalize_genelength(df_temp, gene_length)
df_rpkm.to_csv('rpkm.tsv', sep='\t', header=False)

しかし、このRPKM(FPKM)には1つの大きな落とし穴が潜んでいます。

サンプル間の発現量の比較が可能になるように、まずはじめに総リード数が100万になるようにそろえているのですが、その後でmRNAの長さの違いによるマッピングされるリード数の違いを補正するために、それぞれの転写産物の発現量を表すリード数をmRNAが1000塩基長である場合の値に補正してしまっています。つまり、せっかく総リード数を100万にそろえた後でさらにそのリード数を変化させてしまっているのです。このため最終的な総リード数は100万ではなくなってしまっており、サンプル間の発現量の比較ができなくなってしまっています。

RPKM(FPKM)の問題点は論文(Wagner et al, 2012)でも指摘されており、現在ではその問題点を解消したTPMが主に用いられています。

TPM

RPKM(FPKM)では、総リード数を補正した後でさらに転写産物ごとのリード数を補正してしまったことが問題だったので、転写産物ごとのリードを補正した後で総リード数を補正したのがTPMです。TPMでもRPKM(FPKM)と同様に転写産物を表すmRNAの長さを1000塩基とし、総リード数を100万になるように正規化しています。

それでは、htseq-countの出力結果「count_raw.tsv」をTPMに補正してみましょう。転写産物の長さは「gene_length.tsv」に用意しておきます。

import pandas as pd

# サンプルの総リード数を100万とした場合のリード数に変換する
def normalize_totalreads(df):
    return 10**6 * df / df.sum()

# 転写産物の長さを1000とした場合のリード数に変換する
def normalize_genelength(df, gene_length):
    return (df.T * 10**3 / gene_length).T

# データの読み込み
df_count_raw = pd.read_table('count_raw.tsv', index_col=0, header=None)  
df_gene_length = pd.read_table('gene_length.tsv', index_col=0)
gene_length = df_gene_length['length']

# TPMに補正
df_temp = normalize_genelength(df_count_raw, gene_length)
df_tpm = normalize_totalreads(df_temp)
df_tpm.to_csv('tpm.tsv', sep='\t', header=False)

プログラムを見てわかる通り、TPMは総リード数の正規化と転写産物の長さの正規化の順番をRPKM(FPKM)から逆にしただけです。