FPKM / RPKM

RNA-Seq より得られたリードカウントデータを解析する前に正規化を行うのが一般的である。発現変動遺伝子の検出を目的とする正規化の場合は、TMM 正規化DESeq 正規化などの正規化手法が提唱されている。この他によく利用されている正規化手法には FPKM または RPKM とよばれるものがある。

RNA-Seq によりシーケンシングされたリードは mRNA などのトランスクリプトームの断片である。その断片の元となるトランスクリプトームの長さがながければ長いほど、リード数が多くなる。つまり、短い mRNA からシーケンシングされたリードの数は少なく、長い mRNA からシーケンシングされたリードの数は多いという状況が起こる。場合によって、このような mRNA の長さに起因するリード数の変化を除去する必要がある。そこで、解決法としては、リード数を mRNA 1,000 bp 長さあたりのリード数に変換すれば、mRNA の長さを揃えることができるようになる。このような操作が正規化で、この種の正規化によって得られるデータが FPKM (シングルエンドリードの場合は RPKM) と呼ばれている。FPKM は fragments per kilobase of exon per million reads mapped の略である。

シーケンシングされた全遺伝子のリード数の合計(ライブラリーサイズ)を N とし、そのうち遺伝子 i のリード数を Yi とし、遺伝子 i のトランスクリプトームの長さを Li とすると、遺伝子 i の FPKMi は以下のように計算できる。

\[ FPKM_{i} = Y_{i}\frac{1000}{L_{i}} \frac{1000000}{N} = \frac{Y_{i}}{L_{i}N}10^{9}\]

R を利用して FPKM を求める方法

FPKM を計算するためには、まず遺伝子のカウントデータとその遺伝子の長さの 2 つの情報が必要である。

ここでは 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

次に長さの情報を読み込む。ここで読み込む遺伝子の長さは、実際の転写物の長さではない。そのため、計算される RPKM/FPKM も間違いを含んでいる。実際に、RPKM/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

対数化された FPKM の分布は一般に正規分布に似た形になる(厳密に言えば正規分布にも似てません)。例えば、SRX026633 のライブラリーのデータの分布をヒストグラムにすると以下のようになる。ヒストグラムを描く際に MASS パッケージ中の truehist 関数を利用する。

library(MASS)
truehist(log10(fpkm[, 1]), h = 0.1)
FPKMの分布

edgeR パッケージ中の関数を利用した方法

edgeR パッケージ中に rpkm 関数があり、この関数を利用しても FPKM を計算できる。

rpkm.edger <- rpkm(y, gene.length = gene.len.sorted)
head(rpkm.edger)
##                       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