最尤推定法を利用した異数倍数体の遺伝子発現量定量化手法の開発
理学専攻 情報科学コース

1

恵はじめに網羅的に遺伝子発現量を得る方法として，RNA-seqがある（図 1）．RNA-seq は遺伝子の転写物であるmRNA を大量に読み，既知の遺伝子配列やゲノム配列と対応（アラインメント）させる．対応する配列の本数を数える事で遺伝子の発現量を計測する．近年のシーケンサでは，連続して読める配列（リード）の長さは 100 塩基程度と短いが，1 千万を超える本数が容易に計測可能となったことで，RNA-seq が実現可能となった．一方 RNA-seq には欠点もあり，対象種が持つ複数の遺伝子がほぼ同じ配列を有する場合，リードがどの遺伝子由来か特定できず，発現量の測定を誤る可能性がある．

近縁種：

A.thaliana

1 AAGGTTCC	

既知の遺伝子の配列	既知の遺伝子の配列と	リードが同じ並びの場所を探す（アラインメント）	

同じ場所に貼り付いた	リードの数を数える	
↑３ := 遺伝子発現量 	AAGGTTCC	3 AAGGTTCC	AAGGTTCC	AAGGTTCC	
遺伝子	
転写回数分の転写物が	存在する状態（3回分↑）	

C. flexuosaの遺伝子X	の変異が入った

親：C.amara	  遺伝子Xʼ	の変異が入った

親：C.	  hirsuta	  遺伝子X”	

C.amara由来の遺伝子X	
C. hirsuta由来の遺伝子X	

2本	

C.amara由来	

2本	子：C.	  ﬂexuosa	  遺伝子X”	

C. hirsuta由来	

A.thalianaの遺伝子Xの領域	

遺伝子Xʼ	図 3: 手法の概要

図 2: 実験対象の種
AAGGTTCC	 3	
AAGGTTCC	 2	AAGGTTCC	 1	AAGGTTCC	 0	2 AAGGTTCC	

遺伝子Xの変異が入った

次世代シーケンサで読む	
ひとつのデータをリードと呼ぶ	

図 1: RNA-seq の概要
本研究は，この RNA-seq の利点と欠点を上手に利用することで，異数倍数体における遺伝子発現の定量を行う．異数倍数体とは，異なる二種が交配し，子が両親のゲノムを持つ種で，生物の新規機能獲得プロセスの解明に有用であると考えられているが，両親のゲノム両方とも働くのか，働くとすれば同程度に働くのか，など未知の点が多い．本研究では，これらのゲノムに含まれる遺伝子の発現を同定する手法を確立する．
異数倍数体の例として，本研究では Cardamine 属のC.ﬂexuosa を用いる．この種は C.amara と C.hirsutaを親種に持つ．親種は，同属で進化的距離が近いため，互いに数塩基のみが異なる，ほぼ同一遺伝子集合を有している（図 2）．この種に対し RNA-seq を行うと，上記に述べた欠点によって，遺伝子の発現を定量化することが困難であるので，この問題点を解決する．

2

山田

手法

RNA-seq では，リードをアラインメントするゲノムが必要だが，C.amara と C.hirsuta のゲノムは未知のため，Cardamine 属の近縁で，ゲノムや遺伝子が既知の A.thaliana を利用する．図 2 のように，A.thaliana
は遺伝子 X，C.amara は遺伝子 X’，C.hirsuta は遺伝
子 X”，C.ﬂexuosa は X’，X を持つとすると，X，X’，X は互いに類似した配列を持つ．子である C.ﬂexuosaの X’，X” 由来のリードを考えると，類似の配列を持
つ C.amara の X’ と C.hirsuta の X” の配列両方に対応する．RNA-seq では，このような異なる遺伝子由来の配列が同一の遺伝子に対応することが問題となるが，今回は親種間の遺伝子配列に数塩基の違いにがある事を利用していずれの親由来であるか分離する．

2.1

変異の種類と求める方法

A.thaliana のゲノムと C.ﬂexuosa のリードの間の変異には，親種特有の変異，祖先種から親種へ進化する過程で起きた変異（以下，それぞれ親由来の変異，祖先由来の変異）の２種類がある．図 2 の，C.ﬂexuosa のX’，X”上のｘ印の，青と赤の x 印が前者で，緑の x 印が後者を示す．図 2 の X’，X”のリードを，A.thalianaの X にアラインメントすると図 3 のように表せる．x 印
は図 2 の同じ色のｘ印に相当し，x 印の赤は C.hirsuta
由来の X”の変異と同じで，青は C.amara 由来，緑は祖先由来の変異とみなせる．A.thaliana の X と異なる部分に注目することで，C.ﬂexuosa の X’，X”の変異が得られる．
図示で色が付けられているが，実データでは変異の場所と変異の塩基の種類が分かるのみなので，変異が，親由来か祖先由来かを見分ける方法が必要になる．
2.2

遺伝子発現量の算出

変異の種類に注目すると，リードの組み合わせは２種類ある．親種特有の変異（赤と青の x 印）を見ると，リードの本数は赤と青で 2 本と 2 本に，祖先由来の変異（緑のｘ印）を見ると，緑だけで 4 本と 0 本となる．
各変異が，どちらに由来するものかは分からないので，実際には各変異位置での変異の有無の比率がわかるのみである．祖先由来の変異に関わるリードの本数の比は 0 に近く，親由来の変異に関わるリードの本数の比は，他の親特有の変異の比に近いと考えられる．この特徴を用いて，変異が親由来か祖先由来かを見分ける．
ここでは特定の変異に着目して親種の比率を求めたが1 点のみの推定ではシークエンスエラーによる影響も多く確度が低い．そこで本研究では，リード毎にいずれの親由来であるかを決定し，その本数を遺伝子全体にわたって数える事で計算の精度を向上させる．

2.3

アルゴリズム次の 4 つのステップからなる．
ステップ１では，アラインメントを行い，祖先種の遺伝子配列からの変異を全て抽出する．ステップ２では，抽出した変異から，親由来の変異を求める．ステップ

３では，各リードを親種別に分ける．ステップ４では，片親由来と判定されたリード数を数え，両親のリード数の比と遺伝子の発現量から，親ごとの遺伝子発現量を求める．以下詳細を述べる．
ステップ１：子のリードを祖先種のゲノムに，いくつかの変異を許容して全リードをアラインメントする．以降，遺伝子ごとに計算する．変異が遺伝子上のk 箇所にあるとき，変異を有する位置の集合を M ={M1 , M2 , ..., Mk } で表す．R1i , R2i を，Mi の位置を含むリードを Mi の塩基の種類ごとに分けたリードの集合とする．ここでは変異の有無のみを考えるので，2種類が考えられる．R1i の要素数を |R1i | と表す．
ステップ２：変異 Mi に着目した時，この変異は親由来か祖先由来の 2 種類の可能性がある．RNA-seq はランダムに配列を選択する手法なので，変異 Mi の状態が起こる確率 P (Mi ) は，ポアソン分布を用いてモデル化できる．ポアソン分布は平均 λ, 確率変数 k のとき，λk e−λ /k! と書ける．λp を親由来の平均，λa を祖先由来の平均とすると，Mi が親由来である確率は，|R1 |P (Mi ) = λp i e−λp /|R1i |! ，祖先由来である確率は，|R1i | −λaP (Mi ) = λa e/|R1i |! と表せる．λa は任意の値 t（t <= 1, 実験では 0.1）のとき，λa = t×(|R1i |+|R2i |)で計算する．それぞれの変異は独立な事象であるとす
ると変異 M の尤度 P (M) は P (M1 ) × P (M2 ) × · · · ×
P (Mk ) で表せるので，この尤度が最大になるように，すべての変異を親由来か祖先由来かに分ければ良い．
この時，親由来と考えられる変異セットを MP（ ⊆ M）とする．
予め親種の平均 λp がわかっている場合には，簡単にMP を決定することができるが，λp は求めたい変数であり不明である．つまり，P (M) を最大化するには，最悪全ての考えうる MP 及び λp を試さねばならない．
これには膨大な時間を要する上に，現在考えているのは単一の遺伝子なので，実際にはこの計算を 2 万以上の遺伝子に対して行わなければならない．そこで，現実的な時間で計算できるよう，以下の工夫を行った．
（１）MP に対して，親由来を表すポアソン分布の平均λp は，大体 MP に属する変異の平均に近いと考えられるので，次の方法で計算する．I = {i | Mi ∈ MP }とすると，k = |R1i | のとき，∑()i∈I |R1i |∑λp = (|R1i | + |R2i |) × ∑i∈I |R1i | +i∈I |R2i |（２）|R1i |/|R2i | が小さいもの程，λp を小さくする可能性が高いので，この比率が 1 に近いものから順にMP に入れて，尤度が最大となる組み合わせを計算した．最大尤度を取る I を I とする．
ステップ３：ステップ 2 までで求めた結果を元に，各リードに親ラベル A もしくは B を割り当てる．親ラベルは，|R1i |/(|R1i | + |R2i |) が大きい変異 i ∈ I から順に割り振る．一度親ラベルが決まったリードはラベルを変更せず，次の候補へ移る．候補がなくなるまで繰り返す．
この割当を初期値とし，リード数の多い i ∈ I から順に，割り振ったラベル A, B を反転させ，尤度が大きくなるかを調べる．これは，ラベルが親種を取り違えていないかどうかを確認するための手順である． 以上の手法で，リードごとに由来する親種を特定する．

ス テップ ４：RNA-seq に よ る 遺 伝 子 発 現 量 は
RPKM[1] と呼ばれる．ステップ３で求まった，両親それぞれのリード数の比と RPKM の積が今回求めたい両親由来別の遺伝子発現量となる．

実データを用いた実験の結果

3

次世代シーケンサ SOLiD5 で観測した C.ﬂexuosa のリード約 6 千万本に対し，本手法を適用した．アライ
ンメント先に A.thaliana（TAIR10）[4] のゲノムを選び，ゲノム上の反復配列を省くため RepeatMasker[2]を用い，アラインメントソフトには SHRiMP[3] を用いた．

数	遺伝子の配列	

図 4: 遺伝子上の親由来の変異を持つリードの分布
図 4 は，一つの遺伝子に対し本手法を適用した結果を示している．横軸に遺伝子上の位置，縦軸に塩基ごとのリードの本数を表す．色分けは，リードごとに判定した親種に対応し，赤と青が片親，緑は由来が不明のリードである．この遺伝子の両親種ごとの遺伝子発現量は，それぞれ 5.93, 6.56 だった．グラフからは，親種ごとの発現量を求めるのに使われるリードが，ほぼ遺伝子の範囲上を覆っていて全体的な変異の傾向から発現量を求められていることが分かる．

4まとめ

本研究では，異数倍数体の C.ﬂexuosa の各遺伝子について，親種を区別して遺伝子発現量を計測する手法を開発した．本研究では，対象種である異数倍数体とその親種のゲノムと遺伝子が同定されていない状態で，近縁種との遺伝子セットが類似しているが完全には一致しないという特徴を使う事で，親種別の遺伝子発現量を求める手法を開発した．

謝辞
本研究にあたり，ご指導くださいました東京工業大
学情報理工学研究科 瀬々潤准教授と，貴重なデータとご助言を賜りましたチューリッヒ大学 清水健太郎
教授，清水理恵研究員に深く感謝いたします．

参考文献[1] Mortazavi, et al. Mapping and quantifying mammalian transcriptomes by rna-seq. Nat Methods,
5(7):621-628 (2008).
[2] Smit, AFA, Hubley, R & Green, P. RepeatMasker
Open-3.0. 1996-2004. http://www.repeatmasker.org[3] Rumble, S.M. et al. SHRiMP: accurate mappingof short color-space reads. PLOS Comput. Biol. 5,e1000386 (2009).
[4] http://www.arabidopsis.org/