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

Size: px

Start display at page:

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

ひろとみょうだに
9 years ago
Views:

1 第 5 回バイオインフォマティクスアルゴリズムアラインメントアルゴリズム (3) 慶應義塾大学先端生命科学研究所

2 アラインメント置換挿入欠損を考慮して塩基配列あるいはアミノ酸配列の似た部分をそろえることギャップ - を挿入する CAAGACATTTTAC CATACACTTTAC CA-AGACATTTTAC CATACAC--TTTAC ** * ** *****

3 アラインメントはグラフで表現できる seq2 A A C G A seq1 C G A-CG AACG

4 直前のノード MS(i,j) を位置 (i,j) における最大スコアと定義 MS(i-1,j-1) MS(i-1,j) +SC(seq1(i),seq2(j)) +GP MS(i,j-1) +GP MS(i,j)

5 最大スコアの再帰的定義最大スコアの再帰的定義 = = = 1)) 2( 1), 1( ( 1) 1, ( ) 1, ( 1), ( max or if ), ( j seq i SC seq j i MS GP j i MS GP j i MS j i j i MS

6 動的計画法を使う (C 言語 ) MS(i,j) 各ノード (i,j) について (1) 最大スコア MS(i,j) と (2) 方角 direction(i,j) を記録しておく struct { int score; /* そのノードに至ったときの最大スコア */ int direction; /* そのノードに到達して最大スコアになったときの直前のノードの方角 */ } score_path_matrix[max_seq_len][max_seq_len]; /* MAX_SEQ_LEN は扱える最大の配列長さ */ score_path_matrix[i][j].score = max_score; score_path_matrix[i][j].direction = max_dir;

7 最大スコア方向の情報からアラインメントを求める seq2 A T C G A a_seq1: G a_seq2: G C C - T A A seq1 C G a_seq1: A-CG a_seq2: ATCG

8 void board_to_alignment(char a_seq1[], char a_seq2[], int m, int n){ /* a_seq にアラインメント結果を格納する m, n は 2 つの入力配列のそれぞれの長さ */ int i, j, p; /* i,j が現在注目対象となっている 2 つの配列の塩基の位置 p がアラインメント格納位置 char tmp; } for(p =, i = m, j = n; i!= j!= ;){ switch(score_path_matrix[i][j].direction){ case V: a_seq1[p] = seq1[ -- i]; a_seq2[p] = GAPM; p ++; break; case H:??????????????????? case D:??????????????????? } } a_seq1[p] = ' '; a_seq2[p] = ' '; strrev(a_seq1); strrev(a_seq2); C 言

9 アラインメントの呼び出し (C 言語 ) main(){ static char a_seq1[1], a_seq2[1]; strcpy(seq1, "acg"); strcpy(seq2, "aacg"); seq1_len = strlen(seq1); seq2_len = strlen(seq2); printf("%d n", find_max_score(seq1_len, seq2_len)); } board_to_alignment(a_seq1, a_seq2, seq1_len, seq2_len); printf("%s n%s n", a_seq1, a_seq2);

10 sub board_to_alignment($$){ my($m, $n) /* $m, $n は 2 つの入力配列のそれぞれの長さ */ /* アラインメント結果を格納する配列 */ my($i, $j, $p); /* $i, $j が現在注目対象となっている 2 つの配列の塩基の位置 $p がアラインメント格納位置 */ my($tmp); for($p =, $i = $m, $j = $n; $i!= $j!= ;){ my $d = $score_path_matrix[$i]->[$j]->{"direction"}; if($d eq "V"){ $a_seq1[$p] = substr($seq1, --$i, 1); $a_seq2[$p] = $::GAPM; $p ++; } elsif($d eq "H"){ #????????????????????????? } elsif($d eq "D"){ #????????????????????????? } = = reverse(@a_seq2); } return Perl

11 アラインメントの呼び出し (Perl) $seq1 = "acg"; $seq2 = "atcg"; my $seq1_len = length($seq1); my $seq2_len = length($seq2); printf("%d n", find_max_score($seq1_len, $seq2_len)); my($a_seq1, $a_seq2) = board_to_alignment($seq1_len, $seq2_len); printf("%s n%s n", $a_seq1, $a_seq2);

12 変数の中味の具体例 (1) seq2 i = 1 A A T C G seq1: A C G seq2: A T C G seq1 C j = 2 p = 1 G a_seq1: G C a_seq2: G C

13 変数の中味の具体例 (2) seq2 i = 1 A A T C G seq1: A C G seq2: A T C G seq1 C j = 1 p = 2 G a_seq1: G C - a_seq2: G C T

14 結局どこが動的計画法だったのか? アラインメントの問題を各位置について以下のように分割配列 1 と 2 の塩基をマッチさせる場合配列 1 に GAP を入れる場合配列 2 に GAP を入れる場合各ノードに以下の途中経過を記録し高速化そのノードまでの最大スコア ( そのノードで最大スコアになるときの直前のノード )

15 ローカルアラインメント (1) 配列 1 配列 2 cagtagctgatcgatgctagctgat tttatgat gatgctagctacactac 切り捨てたい部本当に似ている部アラインメントを開始したいポイント

16 ローカルアラインメント (2) seq2 A T C G A seq1 C G 負のスコアに未来はない? 7

17 ローカルアラインメントローカルアラインメント (3) (3) = = = 1)) 2( 1), 1( ( 1) 1, ( ) 1, ( 1), ( max or if ), ( j seq i SC seq j i MS GP j i MS GP j i MS j i j i MS

18 最大スコアと方角の記録 (1) seq2 T A G C C seq1 A G 1 6 4

19 最大スコアと方角の記録 (2) ローカルアラインメントは MS(i,j) が最も高い位置から始める seq2 T A G C seq1 C A 3 a_seq1: GA a_seq2: GA G a_seq1: AG a_seq2: AG

20 より根拠のあるアラインメントのスコア体系ゲノム配列には生物学的にどのような変異が生じやすいかを考えてアラインメントのスコア体系を構築すべき DNA でどのような変異が起こりやすいかスペーサ繰り返し配列プロモータタンパク質のコード領域 ( コドンの第 1~3 塩基 ) Etc. タンパク質の配列でどのような変異が起こりやすいかの方が考慮すべき要素が減り考えやすい疎水性のアミノ酸は他の疎水性のアミノ酸に変化しやすい正電荷を持つアミノ酸は他の正電荷を持つアミノ酸に変化しやい Etc.

21 基本的な考え方 : 同一祖先由来と思われるタンパク質を収集し置換率を考える例 : ロイシンジッパーを含むタンパク質 ACA2_YEAST P4535 KRARLLERNRIAASKCRQRKKVAQLQLQKEFNEIKDENRIL AP1_HUMAN P5412 KAERKRMRNRIAASKCRKRKLERIARLEEKVKTLKAQNSEL K K K L K E 4 箇 1 箇 1 箇置換率に基づいてスコア体系を決めより正確にはアミノ酸 aとbがアラインメントされたときのスコアは aとbの組み合わせが同一祖先からの置換によって生じる確率 log aとbの組み合わせがランダム抽出で得られる確率

22 標準的に使用されているスコア体系 A R N D C Q E G H I L K M F P S T W Y V B Z X * A R N D C Q E G H I L K M F P S T W Y V B Z X * BLOSUM62 PAM Dayhoff et al タンパク質の全長の類似性に基づいて導出された BLOSUM Henikoff and Henikoff, 1992 タンパク質の高い保存性を持つ部分配列に基づいて導出された BLAST など有名なツールにも使用されている

23 より根拠のあるアラインメントのスコア体系 (2) 塩基の挿入欠損はまとめて起こりやすい CCCAAAAAAAAAAGGG CCCAAAAAGGG CCCAAAAAAAAAAGGG CCCA-A-A-A-A-GGG CCCAAAAAAAAAAGGG CCCAAAAA-----GGG どちらが良いアラインメント?

24 より複雑な GAP のペナルティ GAP open penalty GAP の始まりに対して付けられるペナルティー (ex. -5) GAP extension penalty 連続して GAP が入っている領域の 2 番目以降の GAP に対して付けられるペナルティー (ex. -1) CCCAAAAAAAAAAGGG CCCA-A-A-A-A-GGG CCCAAAAAAAAAAGGG CCCAAAAA-----GGG GAP ペナルティーが一定でない再帰式の変更が必要

25 もはや GAP ペナルティーは一定ではない MS(i,j) を位置 (i,j) における最大スコアと定義 MS(i-1,j-1) MS(i-1,j) +SC(seq1(i),seq2(j)) +GP open?,+gp ext? MS(i,j-1) MS(i,j) +GP open?,+gp ext?

26 最大スコアの場合分け MS(i, j) を位置 (i, j) における最大スコアと定義 MS H (i, j) を直前のノードが位置 (i, j-1) だったときの位置 (i, j) における最大スコアと定義 MS D (i,j) MS H (i,j) MS V (i,j) MS(i,j) MS( i, j) = MS MS MS H V D ( i, ( i, ( i, j) j) j)

27 2 つ前まで見れば経路のペナルティーは一定 (i-1,j-2) MS H = MS max MS D H ( i, j ( i, j 1) + GP 1) + GP open ext MS D (i,j-1) (i,j-2) +GP open MS H (i,j) MS H (i,j-1) + GP ext

28 アラインメントの準最適解を求めるもの BLAST (Altschul SF et al. 199) 基本的に GAP なし高速 FASTA (Lipman DJ et al. 1985) GAP を扱うことができるそこそこ速いパッケージ中の ssearch は最適解を求める BLAT (Kent J et al. 22) BLAST より精度が落ちるとされるがとても高速

29 BLAST word の抽出 Query を word(n アミノ酸単位 ) に分割データベース中の word に一致する部分の検出 Hit からの拡張 Hit 部分の統合

30 マルチプルアラインメント 2 本の配列のアラインメントをペアワイズアラインメントという 3 本以上の配列のアラインメントをマルチプルアラインメントというツールとしては CLUSTALW(Thompson JD et al. 1994) などが有名 AC-GT-A ACCAT-A GC--T-A GC--TCA

31 マルチプルアラインメントを行う手法ツリーベース法 ACGTA ACCATA GCTA GCTCA AC-GTA ACCATA GCT-A GCTCA AC-GT-A ACCAT-A GC--T-A GC--TCA シミュレーテッドアニーリング法隠れマルコフモデル

bioinfo-a10s-4_align

bioinfo-a10s-4_align 4 (2) ATCAATCGATCGATC ATCCATTGAACCATC ATCAATCGATCGATC ATCAATCGAGGAGGATCGATC ATCAATCGATCGATC ATCAACGATC - CAAGACATTTTAC CATACACTTTAC CA-AGACATTTTAC CATACAC--TTTAC ** * ** ***** seq2 A A C G A seq1 C G A-CG