RNA-Seqの数理―生成モデルによる発現量推定（3）ペアエンドリードとフラグメント長の分布関数

まとめ：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

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