RNA-Seqの数理―生成モデルによる発現量推定（2）Count Based ModelとRPKM

Count Based Model

今回はCount Based Modelと呼ばれる最も基本的な生成モデルを作って発現量推定を定式化し，それがRPKMの計算式と同じ意味を持つことを示してみよう．なお，以降の内容はPachter 2011をベースにしている．

Models for transcript quantification from RNA-Seq. q-bio.GN, (2011).

http://arxiv.org/abs/1104.3889v2

Count Based Modelにおける仮定

Count based Modelは非常にシンプルな生成モデルだ．単純化のために3つの前提条件を置いて，RNA-Seqをモデル化している．

  • NGSのシーケンスはシングルエンド
  • ショートリードはただ1箇所だけにマッピングされる
  • すべての転写物はただ1つのアイソフォームを持つ

NGSのシーケンスに関しても，トランスクリプトームに関しても，かなり限定して考えていることがわかる．遺伝子の構造が単純な生物において初期のNGS解析をしていると思えばいいだろう．現実的にはこのような状態というのはあり得ないが，転写物からショートリードが生成されるという核の部分は残しつつ単純化をしている．

Count Based Modelの生成過程

生成モデルというものは，まず初めに観測データであるショートリードの生成過程を立てる．今回のCount Based Modelでは

トランスクリプトーム全体からある転写物が選ばれ，ある開始点からショートリードがシーケンスされた

と考える．前回の記事で概要を示したように，RNA-Seqの実験プロトコルは幾つかの複雑な段階を踏んでショートリードをシーケンスしていた．その多くの段階は偶然の連続によって，たまたま転写物がフラグメントとしてある部分から断片化され，たまたまシーケンスされてショートリードとして出力されたとみなすことができる．これらを確率的な事象として扱うことによって，確率的な視点からシーケンスデータが生じた背景にあるRNA-Seqの発現量というものを推定することができる．


さて，長ったらしい文章の前置きはここまでにして，次からは実際に生成過程を数式に落としこんでいく．まず最初にこれまで言葉で説明してきたものを変数や集合として定義をし，事象の起こりやすさを尤度として表現する．概念と数式の対応関係が少し分かりにくいかもしれないが，もし要領が掴めないようならこの後に書いてある具体例の図も参照しつつ，どの変数が何を表しているかを丁寧に追ってほしい．

尤度の定式化

転写物の集合T に含まれる転写物t は，\tilde{l_t} の有効配列長を持ち，\sum_{t\in T} \rho_t = 1 となるような\rho_t という相対的な発現量を持つ．RNA-Seqによって得られたショートリードをf ，その全体をF とし，転写物t にマッピングされたショートリードの集合をF_t とする．

このとき，ショートリードf が転写物t から生じる確率は，以下の式で表すことができる．

p(f\in t)= \frac{\rho_t \tilde{l_t}}{\sum_{r\in T}\rho_r \tilde{l_t}}

次に，選ばれた転写物のどの位置からショートリードがシーケンスされるかはランダムに選ばれるため，この確率p(f\in t) に有効配列長の逆数\frac{1}{\tilde{l_t}} を掛けあわせた値が，ショートリードf が転写物t のある開始点から生成される確率となる．

このとき，観測されたショートリード全体F の尤度は，\rho_t の関数として以下のように表される．

\mathcal{L}(\rho)= \prod_{t\in T}\prod_{f\in F_t} \bigg( \frac{\rho_t \tilde{l_t}}{\sum_{r\in T}\rho_r \tilde{l_r}}\cdot\frac{1}{\tilde{l_t}} \bigg)

そして，転写物t にマッピングされたショートリードの本数X_t を使って，以下のように書き換えられる．

\mathcal{L}(\rho)= \prod_{t\in T} \bigg( \frac{\rho_t \tilde{l_t}}{\sum_{r\in T}\rho_r \tilde{l_r}}\cdot\frac{1}{\tilde{l_t}} \bigg)^{X_t}

これがCount Based Modelにおける尤度である．

具体例を用いて尤度を導出する

確率の考え方

といっても，これだけでは何のことか分からないので，具体例を使って詳しく解説していこう．単純化のために，以下の図のような3つの転写物t_1, t_2, t_3 と6本のショートリードf_1,f_2\dots,f_6 が次のようにマッピングされたという状態を想定する．

3つの転写物にはそれぞれ発現量を持っている．例えば，t_1 は20，t_2 は5，t_3 は15という具合だ．今回は発現量を確率値として考えるため，\frac{20}{40},\frac{5}{40},\frac{15}{40} のように相対値として表す．これが\rho_t であり，総和が1になる理由だ．この\rho_t が今回の生成モデルで最終的に求めたい値である．

さて，ショートリードf_1 が転写物t_1 から生成される確率は，先ほどの定式化でp(f\in t)= \frac{\rho_t \tilde{l_t}}{\sum_{r\in T}\rho_r \tilde{l_t}} としていた．これはある転写物t の長さと相対的な発現量を掛けあわせたものと，その全体との比率を求めており，ただの発現量の比率だけでなく転写物の長さも関係しているということを表している．ランダムに一つ選んでシーケンスするということを考えれば，多く転写されていてなおかつ長い転写物のほうが選ばれやすいのは明らかだろう．この式は，そういった意味を数式として厳密に表現しているだけなのだ．

次に，転写物t が選ばれた後に，今度はその転写物のどの位置の配列が読まれたかということを考える．これは端っこの方の配列が読まれたのか真ん中の方の配列が読まれたのかといった情報だ．このときに出てくるのが有効配列長Effective lengthと呼ばれる値で，ショートリードの長さをm としたときの\tilde{l_t}=l_t - m + 1 で表される．これは単純に長さl_t から長さm の配列を取ってくる時の場合の数である．今回のモデルでは，シーケンスされる場所は場所に依存せず一様なランダムで決まるとしているので，有効配列長で割り算をしている．

そして，これをすべてのショートリードであるf_1 からf_6 まですべて掛けあわせることで，今回得られたショートリードが生成される確率というものを求めることができる．

\mathcal{L}(\rho)=\frac{p(f_1\in t_1)}{\tilde{l_1}}\times\frac{p(f_2\in t_1)}{\tilde{l_1}}\times\frac{p(f_3\in t_1)}{\tilde{l_1}} \times\frac{p(f_4\in t_2)}{\tilde{l_2}}\times\frac{p(f_5\in t_2)}{\tilde{l_2}}\times\frac{p(f_6\in t_3)}{\tilde{l_3}}

この確率\mathcal{L}(\rho) が尤度と呼ばれる．尤度とは「もっともらしさ」という意味だ．上で仮定したようなf_1,f_2\dots,f_6 のショートリードのデータが実際に手元に得られたのだから，それが起こる確率が高かったに違いない，ということを尤度では考えている．

尤度式の展開

では次に，この式をもう少し綺麗にしていこう．まずは2行目で転写物ごとに式をまとめ，3行目ですべての転写物を一つの式で表している．そして最後にp(f\in t) を置き換えると，先ほどの\mathcal{L}(\rho)= \prod_{t\in T}\prod_{f\in F_t} \bigg( \frac{\rho_t \tilde{l_t}}{\sum_{r\in T}\rho_r \tilde{l_r}}\cdot\frac{1}{\tilde{l_t}} \bigg) まで導出することができる．

そして最後のまとめとして，転写物ごとにマッピングされたショートリードの生成確率というのは同じであることを利用して，X_t でまとめる． つまり，ショートリードの配列は生成確率に依存せずp(f_1\in t_1) = p(f_2\in t_1) = p(f_3\in t_1)なので，p(f_1\in t_1) \times p(f_2\in t_1) \times p(f_3\in t_1) = (p(f\in t_1))^{X_1} と変形している．こうすることにより，\mathcal{L}(\rho)= \prod_{t\in T} \bigg( \frac{\rho_t \tilde{l_t}}{\sum_{r\in T}\rho_r \tilde{l_r}}\cdot\frac{1}{\tilde{l_t}} \bigg)^{X_t} となる．

尤度を最大にするパラメータを解析的に求める

では，数式として表された尤度を最大化する値を求めよう．まずは上の式を少し整理するところから始める．\alpha_t という新たな変数を用意して\rho_t を\alpha_t に置き換える．

\alpha_t = p(f\in t)= \frac{\rho_t \tilde{l_t}}{\sum_{r\in T}\rho_r \tilde{l_t}}

この\alpha_t を用いて，先ほどの尤度式を書き換えると

\mathcal{L}(\alpha_t)=\prod_{t\in T}\bigg( \frac{\alpha_t}{\tilde{l_t}}\bigg)^{X_t} \propto \prod_{t\in T} {\alpha_t}^{X_t}

となる．ここで，尤度\mathcal{L} の変数が\alpha_t になっていることに注意．最尤推定の目的としては尤度\mathcal{L} を最大化する\alpha_t を求めたいだけなので，関係ない\tilde{l_t} は消している．

そして，この尤度を最大にする\alpha_t をラグランジュの未定乗数法を用いて求めると，\hat{\alpha_t} = \frac{X_t}{N} となる．以下の図は，ラグランジュの未定乗数法の導出部分の詳細．

このときの\rho_t は\rho_t = \frac{\alpha_t / \tilde{l_t}}{\sum_{r\in T} \alpha_t / \tilde{l_t}} で表されるので，先ほど求めた\hat{\alpha_t} を代入して，

\hat{\rho_t} = \frac{X_t}{N} \cdot \frac{1}{\tilde{l_t}} \cdot \bigg( \frac{1}{\sum_{r\in T}\frac{X_r}{N l_r}} \bigg)

と表すことができる．これが，Count Based Modelによって求められた発現量\rho_t の推定値である．

Count Based Modelの解析解とRPKMとの関係

さて，ようやくRPKMとの関連の話に入ることができる．上で求めた\hat{\rho_t} を，次のように式変形する．

\hat{\rho_t} = \frac{X_t}{N} \cdot \frac{1}{\tilde{l_t}} \cdot \bigg( \frac{1}{\sum_{r\in T}\frac{X_r}{N l_r}} \bigg)

この式をよく見ると，後半の\bigg( \frac{1}{\sum_{r\in T}\frac{X_r}{N l_r}} \bigg) は全ての転写物に関して計算しているだけなので，ただの定数である．ということで，特定の転写物の発現量\hat{\rho_t} に関係してくる変数だけを抜き出すと，

\hat{\rho_t} \propto \frac{X_t}{\tilde{l_t} N}

ここに10^9 を掛ければ，生成モデルにおける発現量の推定式はRPKMの式RPKM \ for \ transcript \ t = 10^9 \times\frac{X_t}{l_t N}と全く同じ式になることがわかる！！！

ということで，Count Based Modelにおける最尤推定によって推定される発現量と，マッピングされたショートリードの数だけで計算されたRPKMの発現量は，基本的には同じ計算式によって求めることができることがわかった．

まとめ

後半はだいぶ駆け足になってしまったが，Count Based Modelの導出とRPKMとの関係を解説してみた．生成モデルというややこしい数式を立てておいて結局RPKMと同じかよと思うかもしれないが，生成モデルはここから様々な仮定をモデルに組み込むことによって，より複雑かつ現実に則した一般的なモデルへと改良していくことができる．今回のCount Based Modelはその枠組みとして，非常に基本的な考え方である．しかしながら，実際にはこの方法だけを使ってソフトウェアが作られたり論文として出ているわけでもなく，あまり詳しい解説が無い部分だと思われる．今回はLior Pachter氏のレビューをもとに具体例をはさみながらモデルを読み解いていったが，次回以降はもう少しCufflinksで使われているような技術も導入しつつ，より複雑なモデルへと進んでいきたい．

参考