Wolfeyes Bioinformatics beta

研究者をネタにした音楽や動画などはいろいろとあるけれども，Life Technologies社のPh.Diddyシリーズに関してあまり日本語で紹介されているウェブサイトがなかったので，ちょっとまとめてみた．このシリーズは珍しく日本語に翻訳されており，歌詞まで字幕で付けられているという充実っぷりだ．今までに3作が公開されており，Ph.Dを目指すPh.DiddyやPh.Divaがウェットのラボで奮闘する姿が描かれている．

Ph.Diddy登場！ (Ph.Diddy is on the scene)

アーノルド・ヤングは、生物工学Ph.D候補です（別名Ph.Diddy）。

アーノルドはA+学生で、現在、サイエンス界の将来を背負って立つ準備中！彼はすぐに発見します。意欲的、献身的な研究室生活で味わう、失敗、達成感、激務、夜遅くまで度重なる実験、査読、そして、尊敬と信頼を勝ち取るための奮闘の日々を！果たして、研究室生活での浮き沈みは、我々のPh.Diddyをどのように形成していくのか？

そして、彼の初めての論文が発表されるまでの旅路は？

http://www.youtube.com/watch?v=9i5M3nLWW78

Ph.Diva と謎のバンド (Ph.Diva and the Mystery Band)

エイプリル・ブライトは生物工学のポスドクです（別名 Ph.Diva）。抄録が口頭発表に選ばれて興奮気味。さっそくプレゼンテーションの準備に取り掛かります。PCR が予想通りの結果にならず、失望と期待、興奮、決断を繰り返す日々が始まりました。果たして、Ph.Divaは晴れて会議の日を迎えることができるのでしょうか？

http://www.youtube.com/watch?v=5H6_-V5iFmI

Ph.Diddy 学会へ行く (Ph.Diddy at the Conference)

アーノルド・ヤングは、生物工学 Ph.D 候補です（別名 Ph.Diddy）。変化の激しいサイエンス界の将来を背負って立つべく、日々奮闘しています。ポスターも準備万端、いよいよ初めての学会に向けて出発です。ラボから飛び出て出張! 飛行機で! 学会に! ドキドキワクワク、テンション MAX! アーノルドは科学に魅せられた仲間に出会い、交流できることを楽しみにしてきました。いろんなブースで貰えるノベルティにも期待しているようです。今回の学会には PI とエイプリル・ブライト（別名 Ph.Diva）も同行しています。アーノルドは隣のブースの女性が気になっています。彼女はアーノルドに気が付いてくれるでしょうか。

http://www.youtube.com/watch?v=W-foZdGdyIM

曲はiTunesで購入できる

ちなみに，これらの曲はiTunesで購入することができる．

The Ph.Diddy Is On the Scene - Single
Ph.Diddy
価格： 150円

Ph. Diva and the Mystery Band - Single
Ph. Diva
価格： 200円

posted with sticky on 2013.7.26

Conference Yeah - Single
Ph.Diddy
価格： 200円

posted with sticky on 2013.7.26

追記

有難いことにライフテクノロジーズジャパンの方に反応を頂いた．どうやら4部作だそうなので，次の作品にも期待！

@yag_ays まとめていただきありがとうございます！！実はPh.Diddyシリーズは4部作らしいのですが、まだ4作目が完成しておりません…。その他にも動画を多数作成しておりますのでYoutubeをチェックしてみてください！ http://t.co/CDp9oZLjvQ
— ライフテクノロジーズジャパン株式会社 (@LifetechJPN) July 26, 2013

特段留学する予定などはなく，単にアメリカの研究事情についての知識が欲しかったので図書館で借りて軽く読み通してみたのだが，かなり内容が濃くて実際に留学する上で必要な知識が詰まった1冊だった．もしアメリカ留学する機会があれば，この本は絶対に買って読み込むだろうと思う．また，ただの留学に必要な手続き紹介だけではなく，ポスドクとしての研究費獲得方法やアメリカと日本のラボの運営形態の違いなども紹介されており，通常の学生や研究者にとっても興味深い内容となっている．

本書はアメリカ合衆国への留学に必要なありとあらゆる情報が網羅された留学ガイド本である．元は「研究留学ネット」という，著者個人がアメリカ留学を期に立ち上げたウェブサイトを書籍化したもので，2002年に第一版が出版され，2012年にビザなどの手続きの変更や特別鼎談などが加わった第二版が出版されている．内容としては著者の個人的な留学経験を元に一般的なアメリカ留学に関して，受け入れ研究室探しの方法から始まり，ビザの発行や保険の加入，銀行口座の開設などの渡航前後に必要となる手続きであったり，日常生活レベルで必要な知識が細かく紹介されている．著者自体が医者のポスドクとして渡航しているからか，主な対象として日本でPhDを取得しアメリカでポスドクをするという人向けの色合いが強い．そのため，配偶者や子供を連れて行く場合のビザの情報であったり，車の購入に関する情報もカバーされている．といっても，当然ながらアメリカ留学に関する基本的な枠組みにそれほど差があるわけではないので，非医学の生命科学系研究者にとっても有用な情報が詰め込まれている．また，本書後半ではアメリカ留学を経験した9人の体験談が収録されており，大学や専門がそれぞれ違う人たちのリアルな留学経験が語られる．そして最後に特別鼎談として，本書著者の門川氏，先ほどの体験談を起稿された方の一人である広田氏，そして「ハーバード大学医学部留学・独立日記」で有名な島岡氏による特別対談が収録されている．

こういった本を読んでいると，将来の選択肢として留学があること，そして国や文化を超えて研究を行うということに関して色々と考えてしまう．かなり個人的な話になってしまうが，今ちょうど同じ研究室にいる中国からの留学生の世話をしているのだけれども，やはり外国で生活をするというのは並々ならぬ労力が必要のようだ．たとえ言語の壁はなんとか突破できたとしても，公共サービスや役所への提出書類の作成であったり，住居などの暮らしに関する細かな違い，そして何より文化の違いは，非常に大きな負荷となるに違いない．これに関しては経験がモノを言う部分だと思うので，今回のような本を参考にして知識を詰め込んでおいたり，実際に留学された経験者が身近にいて相談できる状態というのは非常に大事だと感じた．

参考

拙著『研究留学術』の紹介

まとめ：RNA-Seqの数理―生成モデルによる発現量推定：アーカイブ - Wolfeyes Bioinformatics beta

前回はCount Based Modelの基礎やRPKMとの関連を解説したが，これらはすべてシングルエンドを前提とした発現量推定の方法だった．そこで今回は，現在主流となっているペアエンドに対応できるようモデルの一部を拡張し，その際に必要になるフラグメント長の分布を考慮した有効配列長の考え方について紹介する．

ペアエンドリード

ペアエンドによるシーケンスは，illuminaなどで現在もっとも広く使われているNGSの手法だ．ペアエンドのショートリードは，ゲノムやmRNAを特定の長さのフラグメントに細かく分割したあと，その両端にアダプターを付けて左右同時にシーケンスすることによって得られる．この方法によって，単純にシーケンスのデータ量が2倍になるだけでなく，ある程度離れた距離の2つの配列の関係性も同時に得られる．つまり，元々ひとつの同じフラグメントからシーケンスされたのだから，そのフラグメントの長さの範囲内で2つのショートリードは必ず存在していたはずだという情報である．これらは発現量推定において転写物のアイソフォームを確認するために用いられるほか，主にゲノムアセンブリやトランスクリプトームアセンブリにおいて，コンティグと呼ばれるアセンブルされた断片をつなぎ合わせるスキャフォールディングにおいて有効な実験手法になっている．

RPKMからFPKMへ

シングルエンドではRPKMとして計算していた発現量の指標をペアエンドに拡張するのはとても簡単で，2本のマッピングされたショートリードのペアを1本のフラグメントの本数として見なし，転写物ごとにその本数数えるだけで良い．これをFPKM（Fragments per kilobase of exon per million mapped sequence reads）という．RPKMのReadsの部分をFragments，すなわちペアエンドの2つのショートリードに置き換えただけで，全体のデータ量と転写物の配列長で正規化する部分は変わっていない．ただし，マッピングする配列が実質2本になったということで，両方ともマッピングされるという以外にも，片方だけマッピングされたといった特殊な場合が出てくる．このような特殊なケースはマッピングするソフトウェアの設定やFPKMの計算方法に依存するものの，基本的には2本の配列がフラグメント長の間を置いてマッピングされるかどうかによってFPKMの値を求める．

Count Based Modelにおける変更点

ではCount Based Modelにペアエンドのショートリードを対応させよう．RPKMからFPKMに拡張したときと同様に発現量の導出部分に手を加える必要はないのだが，一つだけペアエンドを考慮する必要のある変数がある．それが転写物の有効配列長 $\tilde{l_t}$ だ．

転写物の有効配列長Effective length of a transcriptsは，ある転写物 $l_t$ から長さ $m$ のショートリードが得られるときの開始点の場合の数であり，シングルエンドの場合は $\tilde{l_t} = l_t - m + 1$ で表される．これはショートリードに長さがあるぶんだけフラグメントが生じる開始点として考慮できる長さは短くなるということだ．

では，今回のペアエンドリードの場合はどう考えればいいだろうか．基本的にはシングルエンドの場合と同様に，長さぶんだけ差し引いた値が有効配列長となる．長さといっても，ペアエンドは元々一つのフラグメントからシーケンスされたショートリードなので，フラグメントの長さを用いる．

しかしここで一つ問題点が出てくる．それは，フラグメントの長さが実は一定ではないということだ．mRNAなどを断片化してある長さのフラグメントに選定する作業は実験的に行なわれるため，どうしてもフラグメントの長さに幅ができてしまう．そのため，ショートリードによって，ペアエンドの間の距離が短いものや長いものが出てくる．この現象は実際にゲノムにペアエンドリードを貼り付けて間の距離を測れば確認することができる．以下の図は，Cufflinksの論文中で実際に求められたフラグメント長の分布を表している（Supplementary Figure 1.）．この場合では，200bpを中心にしてフラグメント長にばらつきがあり，おおよそ正規分布の形をしていることがわかる．山の左側が途中で切断されているのは，ショートリードの長さが75bpのためである．

この図からわかるように，フラグメント長はショートリードによってばらつきがあるため，有効配列長を求める際に一意に値を決めることができない．

フラグメント長の分布関数を使って有効配列長を求める

この問題を解消するための手段として，フラグメント長の分布関数を使うという方法がある．

転写物 $t$ の長さを $l_t$ とする．ショートリードのフラグメントの長さを $i$ としたときのフラグメント長の分布関数を $F(i)$ とし， $\sum_{i=1}^{\infty}F(i)=1$ となるように確率分布に近似する．このときの転写物 $t$ の有効配列長 $\tilde{l_t}$ は，以下のように表される．

$\tilde{l_t} = \sum_{i=1}^{l_t}F(i)(l_t-i+1)$

この式の解釈としては， $(l_t-i+1)$ の部分はシングルエンドの有効配列長の求め方と同じ方法で，ショートリードの長さがフラグメントの長さに変わっただけである．そこに $F(i)$ という分布関数を掛けあわせ，長さ $1$ から転写物の長さ $l_t$ まで足しあわせた値が，ペアエンドの有効配列長となる．つまり，分布の重み付きの配列長として表現しているということだ．

この考え方を模式的に表した図を，以下に示した．ここでは青色の直線で表した転写物 $t$ に黒色の直線で表したショートリードをマッピングしている．太い部分が実際にシーケンスされたショートリード，細い部分がシーケンスされていないフラグメントの部分で，全体で一つのフラグメントを表している．フラグメント長 $i$ を長さ $1$ から $l_t$ までずらしたときのそれぞれの分布は，右に示したような正規分布のような形をとり，これを $F(i)$ としている（実際は離散的）．そして，それぞれのフラグメント長 $i$ における有効配列長 $(l_t-i+1)$ に分布関数 $F(i)$ を掛けあわせたものを全て足し合わせることにより，長さ $l_t$ の転写物における有効配列長 $l_t$ を求めることができる．

Cufflinksとの関連

さて，今回示した数式 $\tilde{l_t} = \sum_{i=1}^{l_t}F(i)(l_t-i+1)$ は，Cufflinksの論文のOnline Methodsに出てくる数式だ．ただしそちらの方では $l_t$ が $l(t)$ となっているが，表記を前回の記事と統一しているだけで表していることは変わりない．

また，Cufflinksのソフトウェアの実行においても，このフラグメント長の分布の種類やパラメータを具体的に指定することができる．

Cufflinks RNA-Seq analysis tools - User’s Manual

例えば，–frag-len-meanはフラグメント長の分布の平均，–frag-len-std-devではフラグメント長の分布の標準偏差を指定することができる．ただし，現在のバージョンでは，

Note: Cufflinks now learns the fragment length mean for each SAM file, so using this option is no longer recommended with paired-end reads.

http://cufflinks.cbcb.umd.edu/manual.html

のように自動的に推定されるようになっており，ユーザが特に意識することなく実行できるようになっている．

まとめ

このように，ペアエンドにおいてはフラグメント長の分布を考慮することによって，有効配列長をより正確に記述した．今回のペアエンドリードの導入ではモデル自体は変更しなかったものの，きちんと実験環境やシーケンスの癖を数式としてモデルに反映することができている．

参考

まとめ：RNA-Seqの数理―生成モデルによる発現量推定：アーカイブ - Wolfeyes Bioinformatics beta

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で使われているような技術も導入しつつ，より複雑なモデルへと進んでいきたい．

参考

まとめ：RNA-Seqの数理―生成モデルによる発現量推定：アーカイブ - Wolfeyes Bioinformatics beta

生成モデルの話を始める前に，まずはRNA-Seqの発現量としてよく知られているRPKMという指標からスタートし，生成モデルの話へと繋げていこう．

RNA-Seqの発現量推定の基本

RNA-Seqで得られるデータは，ある長さの配列（ショートリード）とクオリティスコアのセットして表現される．それをリファレンスとなるゲノム配列の中から探しだして，読まれた配列が何の転写物由来かを見つける作業をマッピングという．このマッピングをNGSから出力された数千万本/数億本のショートリードに対して行い，ある遺伝子の部位にどのくらい貼り付いたかをカウントすることで特定の転写物がどれくらい発現していたかを推定するという方法が，RNA-Seqの発現量推定では基本となる．

この図はRNA-Seqの実験の流れを模式的に表している．まず対象となるトランスクリプトームであるmRNAを取り出し，断片化やアダプターの付与を経てシーケンスライブラリを作成する．次にNGSによって断片が大量にシーケンスされ，いわゆるATCGからなる文字列で構成されたショートリードが得られる．そしてゲノム配列にマッピングすることによって，RNAの発現量が推定される．この図では，ゲノム配列の各塩基ごとにどれくらいの本数のショートリードが貼り付いたかをカウントすることによって，RNA expression levelを図示している．Nucleotide positionによって特定の領域の発現量が欠損していたり遺伝子の末端で発現量が少なくなっているのには，マッピングの際に生じる幾つかの理由がある．まず，図中のJunction readsで示されている点線の箇所は，イントロンのスプライシングが起きた領域であり，mRNAの配列をゲノム配列にマッピングするために起こる現象である．そのため，イントロンの領域だけ発現量が無いように見える．また，Poly-A tailが配列の一部として混入することでショートリード自体がゲノム配列にマッピングできない場合もあり，末端の発現量が落ちる原因の一つになっている．

では，このようにして求められた値を転写物ごとの発現量として一つの値で表現するには，どうすれば良いのか．ここで，RPKMという考え方が登場する．

RPKMの考え方

RNA-Seqの発現量の値として最初に考えられ今でも広く使われているRPKMという指標は，配列の本数を数えるという方法で転写物ごとの発現量を求めている．RPKMはreads per kilobase of exon per million mapped sequence readsの略称で，マッピングされたショートリードの数をエキソンの長さとシーケンサで読まれた配列の総数で正規化した値である．RPKMの計算方法を式で表すと以下のようになる．

$RPKM \ for \ transcript \ t = 10^9 \times\frac{X_t}{l_t N}$

ここで， $X_t$ は転写物 $t$ にマッピングされたショートリードの本数， $l_t$ は転写物 $t$ の長さ(bp)， $N$ はマッピングされたショートリードの総数を表している．

さて，このRPKMという発現量の考え方の裏には以下のような仮定がある．

配列数が長い転写物からはショートリードが多く読まれる
全体で読まれるショートリードの数が多いと，マッピングされるショートリードも多くなる

NGSによるRNA-Seqでは，特定の量のサンプルの中に含まれるトランスクリプトームをランダムにシーケンスするため，出てくるデータというのはあくまでも相対的な量となってしまう．そのため，このような仮定を置くことで転写物の発現や実験サンプルなどの違いを補正している．RPKMはある遺伝子の領域内に含まれるショートリードの本数を数え，転写物の長さと全体のデータ数で割って10⁹ を掛けることによって簡単に求められるため，計算のしやすさと直感的な分かりやすさがある．

RPKMと生成モデルとの関係

このように，RPKMはマッピングされた配列の本数を数えるという単純な方法で発現量を推定している．一方で，これから考えていくことになる生成モデルによる発現量推定は，データが観測される事象を確率的に推定することでデータが与えられたときのトランスクリプトームの発現量というものを求める．これらは一見して全く違うもののように見えるのだが，実際には最も単純な生成モデルはRPKMと基本的な考え方が同じになる．つまり，RPKMも単純な生成モデルによる推定方法も，発現量はマッピングされたショートリードの本数と，転写物の長さおよびマッピングされたショートリードの総数の逆数によって決まるということだ．

それを確かめるために，次回はCount Based Modelと呼ばれる生成モデルを作って発現量推定を定式化し，それがRPKMの計算式と同じ意味を持つということを示してみよう．