アラインメントはグラフで表現できる

Size: px

Start display at page:

Download "アラインメントはグラフで表現できる"

おさむかやぬま
5 years ago
Views:

1 2009 年度春学期バイオインフォマティクスアルゴリズム隠れマルコフモデル (1) 慶應義塾大学環境情報学部

2 この世の現象の多くは偶然性に支配されている? 決定的でない現象の例 : コインを投げたときの裏表サイコロの目天気 ( 降水確率 ) 野球の打者の安打? ( 打率 ) 偶然性にはどのような法則があるのか?

3 確率モデルとは? 確率を取り入れた模型抽象化部分的にはでたらめに見える全体的には秩序がありそれは確率のルールによって決まるサイコロが出る目の確率モデル : P(X) = 1/6

4 確率分布確率モデルを記述する上で有用な概念確率変数 X が取り得る各々の値に対してその起こりやすさを関数 P として記述するサイコロの目を X とすれば 1 が出る確率は P(X=1) = 1/6 より一般的には P(X) = 1/6 コインをトスしたときに出る面を X とすれば表が出る確率 P(X=head) = ½ Etc.

5 確率分布のパラメータ (1) 確率分布の大まかな形 ( 性質?) は決まっているが一意な形が定まっていない場合がある形を一意に定める部分を定数とする式が作られるパラメータ P( 向き = 上 ) = p P( 向き = 下 ) = 1 p

6 確率分布のパラメータ (2) 確率変数が 1 つのときは P(X=i) = p i で多くの確率分布を表現できるサイコロの例では p 1 =1/6, p 2 =1/6, p 3 =1/6, p 4 =1/6, p 5 =1/6, p 6 =1/6 となるゆがんだサイコロでは? パラメータが未知の場合はパラメータの推定が必要

7 確率モデルの構築例 : 画鋲が出る向きモデルの構築 P( 向き = 上 ) = p P( 向き = 下 ) = 1 p パラメータの推定何度も画鋲を投げて向きが上下になる割合を観測する例えば1000 回投げた結果 600 回上向きになったとすれば推定値は P ˆ ( 向き = 上 ) = 0.6 とするのが妥当?

8 試行結果の独立性同じ試行を繰り返すと n-1 回の試行で出た目を X 1, X 2,.. X n-1 として P(X n ) は過去 n-1 回の試行に関係なく 1/6? 6,2,5,1,

9 毎日の天気を試行とすると X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 10 X 11 X 12 天気の移り変わりが穏やかな季節なら n 日目の天気をX n として同じ天気は続きやすい? 曇りの後は雨になりやすい?

10 以前の試行結果が次の試行結果に影響する一般的に雨の確率が20% だとしても前の日が曇りならその確率は 40% になる? 式で表すなら P(X n = Rainy) = 0.2 P(X n = Rainy X n-1 = Cloudy) = 0.4 前の試行結果 X 1,X 2, X n-1 に依存して次の試行結果の確率 P(X n X 1,X 2, X n-1 ) が決まる確率モデルをマルコフモデルと呼ぶ 1つ前の試行結果 X n-1 のみに依存して次の試行結果の確率 P(X n X n-1 ) が決まるモデルを一次のマルコフモデルと呼ぶ

11 一次のマルコフモデルの表現方法条件付確率 (S:Sunny, C:Cloudy, R:Rainy) P(X n =S X n-1 =S)=0.4,P(X n =S X n-1 =C)=0.4, P(X n =S X n-1 =C)=0.2 P(X n =C X n-1 =C)=0.4,P(X n =C X n-1 =S)=0.2, P(X n =C X n-1 =R)=0.4 P(X n =R X n-1 =R)=0.4,P(X n =R X n-1 =C)=0.5, P(X n =R X n-1 =S)= X n X n 状態遷移図行列

12 どんな法則で塩基配列の出現パターンが決まるのか? ggagctgcagcccgaccgcggggaggacgccatcgccgcctgcttcctcatcaactgcct ctacgagcagaacttcgtgtgcaagttcgcgcccagggagggcttcatcaactacctcac gagggaagtgtaccgctcctaccgccagctgcggacccagggctttggagggtctgggat ccccaaggcctgggcaggcatagacttgaaggtacaaccccaggaacccctggtgctgaa ggatgtggaaaacacagattggcgcctactgcggggtgacacggatgtcagggtagagag gaaagacccaaaccaggtggaactgtggggactcaaggaaggcacctacctgttccagct gacagtgactagctcagaccacccagaggacacggccaacgtcacagtcactgtgctgtc caccaagcagacagaagactactgcctcgcatccaacaaggtgggtcgctgccggggctc tttcccacgctggtactatgaccccacggagcagatctgcaagagtttcgtttatggagg ctgcttgggcaacaagaacaactaccttcgggaagaagagtgcattctagcctgtcgggg tgtgcaaggcccctccatggaaaggcgccatccagtgtgctctggcacctgtcagcccac ccagttccgctgcagcaatggctgctgcatcgacagtttcctggagtgtgacgacacccc caactgccccgacgcctccgacgaggctgcctgtgaaaaatacacgagtggctttgacga gctccagcgcatccatttccccagtgacaaagggcactgcgtggacctgccagacacagg actctgcaaggagagcatcccgcgctggtactacaaccccttcagcgaacactgcgcccg ctttacctatggtggttgttatggcaacaagaacaactttgaggaagagcagcagtgcct cgagtcttgtcgcggcatctccaagaaggatgtgtttggcctgaggcgggaaatccccat

13 ゲノム全体は A,C,G,T が均等に出現するランダム配列か? 大腸菌 Accession Number: NC_ Length of Sequence : A Content : (24.62%) T Content : (24.59%) G Content : (25.37%) C Content : (25.42%) Others : 0 (0.00%) AT Content : 49.21% GC Content : 50.79% マイコプラズマ菌 (M. genitalium) Accession Number: NC_ Length of Sequence : A Content : (34.57%) T Content : (33.74%) G Content : (15.91%) C Content : (15.78%) Others : 0 (0.00%) AT Content : 68.31% GC Content : 31.69% Output of G-language GAE (Arakawa et al. 2003)

14 特定の部分での塩基の出現パターンの法則は? ゲノム全体コード領域非コード領域シグナル配列部分

15 核酸配列アミノ酸配列のシグナル開始コドン終止コドン翻訳開始シグナル (SD, Kozak) プロモータ配列タンパクのモチーフ

16 リボソーム Methionine fmet-trnaf 16S rrna mrna AUUCCUCC AUG AGGAGG 開始コドン Shine-Dalgarno sequence 16S rrna の 3 末端は Shine-Dalgarno 配列と対合する

17 SD 配列のパターンは A,G が多いけど曖昧 ttacagagtacacaacatcc atg aaacgcatta [thr operon leader peptide:thrl] aaggtaacgaggtaacaacc atg cgagtgttga [aspartokinase I] atggaagttaggagtctgac atg gttaaagttt [homoserine kinase:thrb] cacgagtactggaaaactaa atg aaactctaca [threonine synthase:thrc] aatgataaaaggagtaacct gtg aaaaagatgc [hypothetical protein:b0005] atttcctgcaaggactggat atg ctgattctta [hypothetical protein:yaaa] gtttaaagagaaatactatc atg acggacaaat [transaldolase B:talB]

18 タンパク質のモチーフ (Helix-loop-Helix) これもアミノ酸の偏りはあるけど曖昧 1 50 YNP2_CAEEL_6-57 AKR..NARER TRVHTVNQAF.LVLKQHLPS.LR Q9W7E6_ KREMVNAKER LRIRNLNTMF.SRLKRMLPL MQ O55208_ RRRVANAKER ERIKNLNRGF.AKLKALVPF LP ASH3_MOUSE_ IRK.RNERER QRVKCVNEGY.ARLRRHLPE.D O76488_ ARR..NARER NRVKQVNDGF.NALRRHLPA.SVVAALS..... AST5_DROME_27-91 IRR..NARER NRVKQVNNGF.SQLRQHIPA.AVIADLSN.... AST4_DROME_ QRR..NARER NRVKQVNNSF.ARLRQHIPQ.SIITDLT YNP2_CAEEL_ Q FTKR.VSKLR ILNAAITYID TLL Q9W7E6_ PDKK.PSKVD TLKAATEYIR LLL O55208_ QSRK.PSKVD ILKGATEYIQ ILG ASH3_MOUSE_ YLEKRLSKVE TLRAAIKYIS YLQ O76488_ GGARRG SGKK.LSKVD TLRMVVEYIR YLQ AST5_DROME_ GRRGIGPG ANKK.LSKVS TLKMAVEYIR RLQ AST4_DROME_ KGGGRG PHKK.ISKVD TLRIAVEYIR RLQ

19 配列パターンのゆらぎ結局ゲノム配列のどこを見ても A,C,G,T が均等にランダムというところはあまりなさそうだシグナル配列もカッチリ決まっているわけではなくゆらぎがあるこれら塩基配列の出現の法則を何とかエレガントに表せないか? 確率モデ

20 塩基配列の出現を表す確率モデルは? モデルの構築 : 各塩基がそれぞれ一定確率で出現 P(A) = P(C) = P(G) = P(T) = ¼ 種によって塩基の偏りがある変異しやすい方向 trna の性質二次構造およびその回避その他たくさんの理由 P(A) = P(T) = 0.34, P(C) = P(G) = 0.16 ではどうか?( マイコプラズマ菌 ) 大腸菌 Accession Number: NC_ Length of Sequence : A Content : (24.62%) T Content : (24.59%) G Content : (25.37%) C Content : (25.42%) Others : 0 (0.00%) AT Content : 49.21% GC Content : 50.79% マイコプラズマ菌 (M. genitalium) Accession Number: NC_ Length of Sequence : A Content : (34.57%) T Content : (33.74%) G Content : (15.91%) C Content : (15.78%) Others : 0 (0.00%) AT Content : 68.31% GC Content : 31.69% Output of G-language GAE (Arakawa et al. 2003)

21 ある塩基の出現率は前の塩基に影響される 1 番目の塩基 2 番目の塩基 2 番目の塩基の出現率 a c A g t a c C g t a c G g t 終止コドンがコード領域中にないから? CpGがTpGに変異? ( 高等生物 ) 前の塩基を考慮するような確率モデルが必要 T a c g t 大腸菌 G-language GAE(Arakawa et al. 2003) 使用

22 マルコフモデルとは? 状態 + 状態間の遷移確率前の状態に依存して確率的に次の状態が決まるゲノム配列を例にとってみると

23 マルコフモデル 0.6 T 1.0 A 0.3 G 1.0 A 0.4 C 0.7 ACTA が観測される確率は?

24 単純なマルコフモデル b1 b2 P(b2 b1) P(A A) A P(C A) P(A C) P(C G) C P(C C) A a c g t a P(G G) P(A T) P(A G) P(G A) G P(G C) P(T G) P(G T) P(C T) P(T A) T P(T C) P(T T) C G c g t a c g t a T c g t 0.298

25 塩基配列の観測確率 (1) = = + = ) ( ) ( ) ( ) ( ) ( ) ( i j j j j i i k k t k P k k t k k t k k t k P L P(k 1 ) を初期状態が k 1 である確率 t(k l) を状態 k から状態 l への遷移確率として経路 k 1,k 2,k 3,,k i を通過する確率は A G T C A

26 塩基配列の観測確率 (2) P(A A) P(C C) A G P(C A) P(A C) P(C G) C P(A T) P(C T) P(T C) P(A G) P(G A) P(T A) P(G C) P(T G) P(G T) P(G G) P(T T) T 状態とそこから出力される塩基の種類が 1 対 1 対応の場合長さ i の塩基配列 x 1,x 2,x 3,,x i が観測される確率は P( x ) P( x x ) P( x x ) L P( x 1 = P( x j = i 1) j= 2 2 P( x 1 j x j 1 ) 3 2 i x i 1 )

27 どこまで前の状態が影響するか? 0 次のマルコフモデル一次のマルコフモデル n 次のマルコフモデル多くのゲノムの非コード領域は 1 次のマルコフモデルで近似されることが多い

28 隠れマルコフモデルとは? 状態とその間をつなぐ遷移確率がある状態から出力される記号は確率的に決まる

29 隠れマルコフモデル A:0.2 T:0.2 C:0.2 G: A:0.3 T:0.2 C:0.4 G: A:0.3 T:0.4 C:0.2 G: A:0.1 T:0.2 C:0.4 G:0.3 A:0.3 T:0.2 C:0.1 G:0.4

30 隠れマルコフモデルは以下の 2 つで定義される t(k l) 状態 k から状態 l への遷移確率 e l (b) 状態 l において記号 b が出力される確率上記モデルから長さ i の塩基配列 x 1,x 2,x 3,,x i が出力される

31 隠れマルコフモデルを使って配列パターンをモデル化したい構造をどうする? 状態遷移確率記号出力確率をどう決める? 与えられた配列の出力確率をどう求める? 与えられた配列の通った経路をどう求める?

32 隠れマルコフモデルの構造の一例欠損挿入 start end A:0.7 C:0.1 G:0.1 T:0.1 A:0.0 C:0.7 G:0.2 T:0.1 A:0.0 C:0.9 G:0.1 T:0.0 A:0.0 C:0.1 G:0.1 T:0.8

33 与えられた配列の出力確率は? 経路が分かっている場合 A:0.2 T:0.2 C:0.2 G:0.4 P( ATCC ) =0.3 x 0.4 x 0.4 x 0.7 x 0.4 x 1.0 x A:0.3 T:0.2 C:0.4 G: A:0.3 T:0.4 C:0.2 G: A:0.1 T:0.2 C:0.4 G:0.3 A:0.3 T:0.2 C:0.1 G:0.4

34 与えられた配列の出力確率は? 経路が分かっていない場合 Forward アルゴリズム Backward アルゴリズム

35 与えられた配列の通った経路はどうやって求める? 通った経路は確率的にしか決まらない最大の確率を与えるような経路は Viterbi のアルゴリズムで求める

36 Viterbi のアルゴリズム A:0.4 C:0.1 G:0.3 T:0.2 A:0.2 C:0.2 G:0.3 T:0.3 ATGA が観測されたときどの経路を通った可能性が最も高いか? A:0.1 C:0.1 G:0.1 T: A:0.1 C:0.3 G:0.5 T: A:0.7 C:0.1 G:0.1 T:0.1

37 v(l, i) の定義 v(l, i) 状態 l において塩基配列 x 1,x 2,..x i が観測されているときの最大確率 1.0 A:0.4 C:0.1 G:0.3 T: A:0.2 C:0.2 G:0.3 T: A:0.7 T:0.1 C:0.1 G: A:0.1 C:0.1 G:0.1 T: A:0.1 C:0.3 G:0.5 T: 今 atga が観測されたとして v(3,2) は状態 3 の時点で at が観測されている最大確率

38 k 1 e l (x i ) k 2 e kmax (x i-1 ) t(k max l) l v(l, i) k max p max = v(k max, i-1) k 4 状態 l につながる直前の状態 k v(l,i) を与える直前の状態を k max そこまでの確率を p max とすると v(l,i)=p max t(k max l) e l (x i ) 各々の k に対して t(k l) は一定

39 Viterbi のアルゴリズム = = = = )) ( 1), ( ( ) max ( 0 0 and 1if 0 0 and 0 if ), ( l k t i k v b e l i l i i l v k i l

40 最大確率と直前のノードの記録 k 2 k 3 k 1 各状態 l と i 番目について (1) 最大スコア v(l,i) と (2) 方角 ptr(l,i) を記録しておく l v(l,i)

41 atga はどの経路を通った? A:0.4 C:0.1 G:0.3 T:0.2 A:0.2 C:0.2 G:0.3 T: A:0.1 C:0.1 G:0.1 T: A:0.1 C:0.3 G:0.5 T: A:0.7 T:0.1 C:0.1 G:0.1

42 Viterbi のテーブル v(l,i) 配列の i 番目状態 l

43 Viterbi のテーブル ptr(l, i) 配列の i 番目状態 l 経路 :

44 演習問題下記マルコフモデルで agga が出力されるときの最大確率およびその経路を求めよ v(l,i), ptr(l,i) を記録すること 1.0 A:0.4 C:0.1 G:0.3 T: A:0.2 C:0.2 G:0.3 T: A:0.1 C:0.1 G:0.1 T: A:0.1 C:0.3 G:0.5 T: A:0.7 T:0.1 C:0.1 G:0.1

45 演習問題解答 v(l, i)

46 演習問題解答 ptr(l, i)

47 実装 (C 言語 ) 塩基配列 0 static char x[l]; strcpy(x, "atga"); 状態遷移確率 1.0 A:0.4 C:0.1 G:0.3 T: A:0.1 C:0.1 G:0.1 T: static double t[k][k] = { /* */ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 ], /* 0 */ 0.0, 0.0, 0.8, 0.2, 0.0, 0.0 ], /* 1 */ 0.0, 0.0, 0.5, 0.0, 0.5, 0.0 ], /* 2 */ 0.0, 0.0, 0.0, 0.0, 0.7, 0.3 ], /* 3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 ], /* 4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ] /* 5 */ }; A:0.2 C:0.2 G:0.3 T: A:0.1 C:0.3 G:0.5 T: 記号出力確率状態のタイプ static int type[k] = ( "TYPE_S", "TYPE_N", "TYPE_N", "TYPE_N", "TYPE_N", "TYPE_N" ); 5 A:0.7 C:0.1 G:0.1 T:0.1 static double e[k][nuc_num] = { /* a c g t */ [ 0.00, 0.00, 0.00, 0.00 ], /* 0 */ [ 0.40, 0.10, 0.30, 0.20 ], /* 1 */ [ 0.10, 0.10, 0.10, 0.70 ], /* 2 */ [ 0.20, 0.20, 0.30, 0.30 ], /* 3 */ [ 0.10, 0.30, 0.50, 0.10 ], /* 4 */ [ 0.70, 0.10, 0.10, 0.10 ] /* 5 */ }; K: 状態の数 L: 配列の長さ NUC_NUM: 塩基の種類数

48 Viterbi 関数の実装例 (C 言語 ) double viterbi(int l, int i){ int k, k_max; double p, p_max; switch(type[l]){ case TYPE_S: if(i <= 0){ return 0.0; } else { return LOG0; } case TYPE_N: if(i <= 0){ return LOG0; } } } for(p_max = LOG0-1, k = 0;k < K;k ++){ if(tr[k][l] > 0.0)p = viterbi(k, i - 1) + log(tr[k][l]); else p = LOG0; if(p > p_max){ p_max = p; k_max = k; } } return p_max + log(e[l][ (int)cton(x[i - 1]) ]);

49 実装 (Perl) 塩基配列 0 = split("", "atga"); 状態遷移確率 A:0.4 C:0.1 G:0.3 T: A:0.1 C:0.1 G:0.1 T: = ( # [ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 ], # 0 [ 0.0, 0.0, 0.8, 0.2, 0.0, 0.0 ], # 1 [ 0.0, 0.0, 0.5, 0.0, 0.5, 0.0 ], # 2 [ 0.0, 0.0, 0.0, 0.0, 0.7, 0.3 ], # 3 [ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 ], # 4 [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ] # 5 ); A:0.2 C:0.2 G:0.3 T: A:0.1 C:0.3 G:0.5 T:0.1 = ( "TYPE_S", "TYPE_N", "TYPE_N", "TYPE_N", "TYPE_N", "TYPE_N" ); 記号出力確率 5 A:0.7 C:0.1 G:0.1 = ( # a c g t [ 0.00, 0.00, 0.00, 0.00 ], # 0 [ 0.40, 0.10, 0.30, 0.20 ], # 1 [ 0.10, 0.10, 0.10, 0.70 ], # 2 [ 0.20, 0.20, 0.30, 0.30 ], # 3 [ 0.10, 0.30, 0.50, 0.10 ], # 4 [ 0.70, 0.10, 0.10, 0.10 ] # 5 );

50 Viterbi 関数の実装例 (Perl) sub viterbi($$){ my($l, $i) my($k, $k_max); my($p, $p_max); } if($type[$l] eq "TYPE_S"){ if($i <= 0){ return 0.0; } else { return $::LOG0; } } elsif($type[$l] eq "TYPE_N"){ if($i <= 0){ return $::LOG0; } for($p_max = $::LOG0-1, $k = 0;$k < $::K;$k ++){ if($tr[$k]->[$l] > 0.0){ $p = viterbi($k, $i - 1) + log($tr[$k]->[$l]); } else { $p = $::LOG0; } if($p > $p_max){ $p_max = $p; $k_max = $k; } } return $p_max + log($e[$l]->[ $cton{ $x[$i - 1] } ]); }

第４回バイオインフォマティクスアルゴリズム実習

第４回バイオインフォマティクスアルゴリズム実習第 5 回バイオインフォマティクスアルゴリズムアラインメントアルゴリズム (3) 慶應義塾大学先端生命科学研究所アラインメント置換挿入欠損を考慮して塩基配列あるいはアミノ酸配列の似た部分をそろえることギャップ - を挿入する CAAGACATTTTAC CATACACTTTAC CA-AGACATTTTAC CATACAC--TTTAC ** * ** ***** アラインメントはグラフで表現できる