リードカウントの分布 | 負の二項分布に従うRNA-seqのカウントデータ

RNA-seq の実験で、発現変動遺伝子を同定するとき、最終的に RNA-seq のデータをリードカウントデータあるいは FPKM の形で出力する。ここでは、リードカウントデータの方について述べる。

1、2、3、・・・と数え上げるデータはカウントデータ（計数データ）と呼ばれる。カウントデータをモデル化するときは、離散値の確率分布であるポアソン分布、負の二項分布などが用いられる。RNA-seq の実験より得られたリードカウントデータもリードカウントの 1 種であるため、これらのモデルでモデル化されます。しかし、biological replicate の RNA-seq データの場合、ポアソン分布を用いると過分散とよばれる現象が生じる。そこで現在はこのようなデータを解析するとき、ポアソン分布よりも自由度の高い負の二項分布を用いるのが一般的である。もしろん、ポアソン分布を拡張した分布なども用いられるデータ解析方法もある。

ポアソン分布と負の二項分布

リードカウントをポアソン分布を用いて解析するとき、そのモデルは以下のように立てることができる。ただし、リードカウントのリード数を Y とする。

\[ P(Y=y) = \frac{\lambda ^{y}e^{-\lambda}}{y!} \]

ポアソン分布を利用したモデルのパラメーターでは λ 一つだけであり、その特徴から「λ = 平均 = 分散」の関係が成り立つ。そのため、「分散 > 平均」のような状況に対応できない。そこで、ガンマ分布 λ が Gamma(φ^-1, μφ) に従うものとして、Y の確率分布関数は次のように拡張できる。

\[ \begin{eqnarray} P(Y=y|\mu, \phi) &=& \int_{0}^{\infty}P(Y|\lambda)P(\lambda |\mu,\phi)d\lambda \\ &=& \int_{0}^{\infty} \frac{\lambda^{y}e^{-\lambda}}{y!} \frac{\lambda^{\phi^{-1}-1}e^{-\lambda\mu^{-1}\phi^{-1}}}{(\mu\phi)^{\phi^{-1}}\Gamma (\phi^{-1})}d\lambda \\ &=& \frac{\Gamma(y + \phi^{-1})}{\Gamma(\phi^{-1})\Gamma(y+1)}\left(\frac{1}{1+\mu\phi}\right)^{\phi^{-1}}\left(\frac{\mu}{\phi^{-1} + \mu}\right)^{y} \end{eqnarray} \]

このように導出された分布関数が負の二項分布の確率分布関数となる。

負の二項分布

負の二項分布は、確率 p で成功するベルヌーイ試行を行ったとき、r 回の成功するまでの、失敗した回数の分布である。したがって、負の二項分布は、試行の成功確率 p と成功回数 r の二つのパラメーターを必要とする。そこで、失敗した回数の確率分布 P(Y=y) は、成功確率 p かつ成功回数 r の負の二項分布を用いて次のように書ける。

\[ P(Y=y) = \frac{\Gamma (y+r)}{\Gamma (r) \Gamma (y+1)} p^{r}(1-p)^{y} \]

この時、負の二項分布の期待値および分散は r と p を用いて、次のように表すことができる。

\[ E(Y) = \frac{r(r+1)(1-p)^{2}}{p^{2}} \] \[ Var(Y) = \frac{r(1-p)}{p^{2}}\]

このように、負の二項分布の分散は期待値よりも大きいことが確認できる。

リードカウントデータと負の二項分布

上述のように、負の二項分布は成功確率 p と成功回数 r によって決定される。リードカウントデータの場合は、成功確率と成功回数に関して興味はない。生物学的には興味があるのは期待値と分散である。そこで、まず、期待値を μ とおく。（もちろん、真の μ は不明。）

\[ E(Y) = \mu \]

次に分散を期待値よりも大きくなるようにモデル化する。例えば、φ > 0 となるパラメーター φ を導入し、次のように決めることができる。

\[ Var(Y) = \mu + \phi\mu^{2}\]

分散の決め方は期待値よりも大きければよく、上記のほかにも様々な決め方が提唱されている。例えば、以下のような決め方がある。

分散	論文や採用しているRのパッケージ
\[Var(Y) = \mu + \phi\mu^{2}\]	Robinson らのモデル。edgeR パッケージに採用されている。
\[Var(Y) = \mu + s^{2}\nu_{\rho}\]	Anders らのモデル。DESeq パッケージに採用されている。第二項は raw variance と呼ばれている、平均に依存するなめらかな関数である。

こうして、期待値と分散が μ と φ によって表されたとする。続いて、期待値と分散を利用して、r と p を μ と φ で表すと、次のようになる。（r と p を次のように変換するのは上述したポアソン分布から導いた二項分布の確率分布関数にフィットさせるためである。）

\[ r = \phi^{-1} \] \[ p = \frac{1}{1+\mu\phi} \]

この時、負の二項分布の分布関数は、μ と φ を用いて次のように書ける。この式が Robinson らの文献によく見られるリードカウントの分布関数になる。

\[ P(Y=y) = \frac{\Gamma(y + \phi^{-1})}{\Gamma(\phi^{-1})\Gamma(y+1)}\left(\frac{1}{1+\mu\phi}\right)^{\phi^{-1}}\left(\frac{\mu}{\phi^{-1} + \mu}\right)^{y} \]

References

Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010, 11(10):R106. PubMed Abstract
Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010, 11(3):R25. PubMed Abstract