RNA-Seq データから得られたリードカウントデータは、そのまま転写物(遺伝子)発現量を表すわけではない。RNA-Seq によりシーケンシングされたリードは、mRNA などの転写産物の断片である。その断片の元となる転写産物の長さがながければ長いほど、リード数が多くなる。つまり、短い転写産物からシーケンシングされたリードの数は少なく、長い 転写産物からシーケンシングされたリードの数は多いという状況が起こる。遺伝子間の発現量比較を行うのが目的であれば、このような転写産物の長さに由来する影響を取り除く必要がある。その解決方法の一つとして、転写産物 1,000 bp あたりのリード数を求めればよい。このような操作を正規化という。RNA-Seq データから得られたリードカウントデータに対して、総リード数による補正の後に、転写産物長による補正を行なって得られるデータを FPKM(single-end リードの場合は RPKM)という。FPKM は fragments per kilobase of exon per million reads mapped の略である。
ただし、FPKM/RPKM は転写産物の発現量を正しく表せないことが報告されており(Wagner et al, 2012)、最近では、FPKM/RPKM の代わりに TPM が用いられるようになった。
FPKM/RPKM の計算
リファレンスにマッピングできた全リード数を N、そのうち転写産物 i の領域にマッピングされたリード数を Yi、転写産物 i の長さを Li とすると、転写産物 i の FPKMi は以下のように計算できる。
この式をわかりやすく整理すると、次のように書き換えることができる。
最初の項は、ある転写物の長さ 1,000 bp あたりのリード数を表している。つまり、長い転写物も、短い転写物も同じ長さになるように補正してからリード数を計算している。この補正後のリード数は、相対発現量とみなすことができ、同じライブラリー内での遺伝子の発現量の比較を可能にしている。次の項は、そのライブラリーの総リード数を揃える操作を行なっている。高速シーケンサーの出力結果では、ライブラリーによって総リード数が大きく異なる。そのため、ライブラリー同士(サンプル同士)の比較を可能にするためには、総リード数が同じとなるように補正する必要がある。PFKM はこのような操作で計算される。
解析対象が遺伝子の場合でも FPKM を計算することができる。この場合、遺伝子の長さを解析目的に応じて適当に定義する必要がある。例えば、遺伝子領域内の全 exon の長さの和(重複除く)を遺伝子長と定義したりすることができる。
FPKM/RPKM の不正確さ
上述のように、FPKM の計算式を長さによる補正部分と総リード数による補正部分を分けて、それぞれ考えたときに、正しいように見える。しかし、両者を合わせて考えたときに、間違いが生じるようになる。特に、分母の方の N/106 に着目すると、この値はサンプル中の転写物の総量を表すべきである。しかし、実際には、N/106 は RNA-Seq シーケンシングの深さによって決定され、サンプル中の転写物総量と線形関係を示さない。その理由として、細胞の中で、長い遺伝子と短い遺伝子が転写される頻度が異なっているから、N(あるいは N/106)と転写物総量の間に、遺伝子長に由来する要因も絡んでいるためである。
次に、具体的な計算を使って FPKM をみていく。もし発現量を正確に定量できるのであれば、定量に用いた RNA サンプルが同じであれば、定量手法の違いに関わらず、同じ遺伝子であれば同じ発現量が得られる。ここで、あるサンプル中に存在する全転写産物の個数を N 個、このうち転写物 i に由来する転写産物を ni 個と仮定する。このとき、このサンプル中における転写産物 i の存在率 ri は、ri = ni/N と計算される。同様に、転写産物 j の割合も rj = nj/N と計算される。次に、サンプル中に含まれる計 T 個の転写産物の存在率の平均値を計算する。
\[ \bar{r} = \frac{1}{T}\sum_{i}^{T}r_{i} = \frac{1}{T}\frac{\sum_{i}^{T}n_{i}}{N} = \frac{1}{T}\frac{N}{N} = \frac{1}{T} \]この式からみられるように、全転写産物の存在率の平均値は、総遺伝子数の影響のみを受ける。つまり、同一生物種ならば、この値は変化しない。たとえ、正常細胞とがん細胞の比較、あるいは組織間の比較において、種が同じならば、この全転写産物の存在率の平均値は変化しないはずである。
転写産物の発現量定量において、本来ならば各転写産物に対してこの存在率を求めていくことになる。しかし、現実問題として、それが難しい。そこで、カウントデータに対して正規化を行い、間接的に各転写産物の存在率を推定していく。FPKMは、存在率を間接的に表す値の 1 種類である。この PFKM が、もし正確に転写産物の存在率(あるいはその定数倍)を表しているならば、全転写産物の FPKM の平均値は、サンプルの種類に関わらず一定でなければならない。
そこで、転写産物が 5 種類のみからなるサンプル A と B のデータに対して、FPKM を計算し、全転写産物の FPKM の平均値を計算してみることにする。
t1 | t2 | t2 | t4 | t5 | |
転写物長 | 1000 bp | 1200 bp | 800 bp | 500 bp | 2000 bp |
サンプル A | 20 | 30 | 20 | 50 | 150 |
サンプル B | 100 | 60 | 140 | 120 | 180 |
このサンプルデータに対して PFKM を計算すると次のようになる。サンプル A とサンプル B の FPKM の平均値について計算すると、それぞれ 181481.5、218333.3 となり、まったく違う値が得られた。このことは、FPKM が各転写産物の存在率を正確に表していないことを示す。つまり、FPKM が転写産物の発現量を正確に表していないことを示している。
t1 | t2 | t2 | t4 | t5 | |
サンプル A | 74074.07 | 92592.59 | 92592.59 | 370370.37 | 277777.78 |
サンプル B | 166666.67 | 83333.33 | 291666.67 | 400000.00 | 150000.00 |
一方で、TPM と呼ばれる正規化法を用いると、サンプル A およびサンプル B の FPKM の平均量はともに 106/5 と計算される。つまり、転写産物の総量の定数倍(103)となっている。このとき、サンプル A およびサンプル B は同じ基準となり、両者の比較が可能になる。
R を利用して FPKM を計算する方法
FPKM を計算するためには、マッピングソフトから出力されるカウントデータと遺伝子長の情報が必要である。ここでは ReCount データベースに公開されているデータを例に FPKM を計算する。まず、カウントデータを読み込んで、整形する。
data <- read.table("http://bowtie-bio.sourceforge.net/recount/countTables/katzmouse_count_table.txt", header = TRUE)
## gene SRX026633 SRX026632 SRX026631 SRX026630
## 1 ENSMUSG00000000001 842 765 437 221
## 2 ENSMUSG00000000003 0 0 0 0
## 3 ENSMUSG00000000028 42 60 10 17
## 4 ENSMUSG00000000031 0 0 0 0
## 5 ENSMUSG00000000037 0 0 0 0
## 6 ENSMUSG00000000049 1 0 0 1
y <- data[, -1]
rownames(y) <- data[, 1]
head(y)
## SRX026633 SRX026632 SRX026631 SRX026630
## ENSMUSG00000000001 842 765 437 221
## ENSMUSG00000000003 0 0 0 0
## ENSMUSG00000000028 42 60 10 17
## ENSMUSG00000000031 0 0 0 0
## ENSMUSG00000000037 0 0 0 0
## ENSMUSG00000000049 1 0 0 1
次に長さの情報を読み込む。ただし、ここで読み込む遺伝子の長さは、イントロンなどの転写されない領域の長さも含まれている。そのため、この遺伝子長を用いて計算した FPKM は間違っている。実データの解析に当たる場合は、必ず正確な遺伝子長または転写物長を用いること(遺伝子長の定義と計算方法)。
data.len <- read.table("http://bowtie-bio.sourceforge.net/recount/ivals/mouse_genes.txt", header = TRUE)
gene.len <- data.len[, 5] - data.len[, 4] + 1
names(gene.len) <- data.len[, 1]
head(gene.len)
## ENSMUSG00000087647 ENSMUSG00000081805 ENSMUSG00000081824 ENSMUSG00000086074
## 8642 810 319 5457
## ENSMUSG00000081968 ENSMUSG00000074467
## 474 425
カウントデータ中の遺伝子名の順番と遺伝子の長さのベクトルの遺伝子名の順番が異なるので、ここで遺伝子長ベクトルをソートしておく。
gene.list.order <- rownames(y)
gene.len.sorted <- gene.len[gene.list.order]
head(gene.len.sorted)
## ENSMUSG00000000001 ENSMUSG00000000003 ENSMUSG00000000028 ENSMUSG00000000031
## 38867 15723 31541 2615
## ENSMUSG00000000037 ENSMUSG00000000049
## 141021 71043
念のため、遺伝子発現量行列中の遺伝子数と遺伝子長ベクトル中の遺伝子数が一致しているかどうかをチェックする。
nrow(y)
## [1] 36536
length(gene.len.sorted)
## [1] 36536
最後に FPKM を計算する。
cpm <- sweep(y, 2, 1e6 / colSums(y), "*")
fpkm <- sweep(cpm, 1, 1e3 / gene.len.sorted, "*")
head(fpkm)
## SRX026633 SRX026632 SRX026631 SRX026630
## ENSMUSG00000000001 14.579167320 10.0588002 13.8070378 8.1219532
## ENSMUSG00000000003 0.000000000 0.0000000 0.0000000 0.0000000
## ENSMUSG00000000028 0.896139212 0.9721685 0.3893361 0.7698794
## ENSMUSG00000000031 0.000000000 0.0000000 0.0000000 0.0000000
## ENSMUSG00000000037 0.000000000 0.0000000 0.0000000 0.0000000
## ENSMUSG00000000049 0.009472843 0.0000000 0.0000000 0.0201061
References
- Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. 2012, 131(4):281-5. DOI: 10.1007/s12064-012-0162-3
- 为什么说 FPKM 和 RPKM 都错了? 简书 2019, website