第 5 回バイオインフォマティクスアルゴリズム アラインメントアルゴリズム (3) 慶應義塾大学先端生命科学研究所
アラインメント 置換 挿入 欠損を考慮して塩基配列あるいは アミノ酸配列の似た部分をそろえることギャップ - を挿入する CAAGACATTTTAC CATACACTTTAC CA-AGACATTTTAC CATACAC--TTTAC ** * ** *****
アラインメントはグラフで表現できる seq2 A A C G A seq1 C G A-CG AACG
直前のノード 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)
最大スコアの再帰的定義最大スコアの再帰的定義 + + + = = = 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
動的計画法を使う (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;
最大スコア方向の情報からアラインメントを求める seq2 A T C G A a_seq1: G a_seq2: G C C - T A A seq1 C 3 1-1 -2 1 4 2 G -1-2 2 7 a_seq1: A-CG a_seq2: ATCG
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 言
アラインメントの呼び出し (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);
sub board_to_alignment($$){ my($m, $n) = @_; /* $m, $n は 2 つの入力配列のそれぞれの長さ */ my(@a_seq1, @a_seq2); /* アラインメント結果を格納する配列 */ 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"){ #????????????????????????? } } @a_seq1 = reverse(@a_seq1); @a_seq2 = reverse(@a_seq2); } return (join("", @a_seq1), join("", @a_seq2)); Perl
アラインメントの呼び出し (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);
変数の中味の具体例 (1) seq2 i = 1 A A T C G 1 2 3 seq1: A C G seq2: A T C G seq1 C 3 1 1-1 -2 4 2 j = 2 p = 1 G -1-2 2 7 1 a_seq1: G C a_seq2: G C
変数の中味の具体例 (2) seq2 i = 1 A A T C G 1 2 3 seq1: A C G seq2: A T C G seq1 C 3 1 1-1 -2 4 2 j = 1 p = 2 G -1-2 2 7 1 2 a_seq1: G C - a_seq2: G C T
結局どこが動的計画法だったのか? アラインメントの問題を 各位置について 以下のように分割 配列 1 と 2 の塩基をマッチさせる場合 配列 1 に GAP を入れる場合 配列 2 に GAP を入れる場合 各ノードに以下の途中経過を記録し 高速化 そのノードまでの最大スコア ( そのノードで最大スコアになるときの直前のノード )
ローカルアラインメント (1) 配列 1 配列 2 cagtagctgatcgatgctagctgat tttatgat gatgctagctacactac 切り捨てたい部 本当に似ている部 アラインメントを開始したいポイント
ローカルアラインメント (2) seq2 A T C G A seq1 C 3 1-1 -2 1 4 2 G -1-2 2 負のスコアに未来はない? 7
ローカルアラインメントローカルアラインメント (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
最大スコアと方角の記録 (1) seq2 T A G C C seq1 A 3 3 1 1 G 1 6 4
最大スコアと方角の記録 (2) ローカルアラインメントは MS(i,j) が最も高い位置から始める seq2 T A G C seq1 C A 3 a_seq1: GA a_seq2: GA G 3 1 1 1 6 4 a_seq1: AG a_seq2: AG
より根拠のあるアラインメントのスコア体系 ゲノム配列には生物学的にどのような変異が生じやすいかを考えてアラインメントのスコア体系を構築すべき DNA でどのような変異が起こりやすいか スペーサ 繰り返し配列 プロモータ タンパク質のコード領域 ( コドンの第 1~3 塩基 ) Etc. タンパク質の配列でどのような変異が起こりやすいかの方が考慮すべき要素が減り 考えやすい 疎水性のアミノ酸は他の疎水性のアミノ酸に変化しやすい 正電荷を持つアミノ酸は他の正電荷を持つアミノ酸に変化しやい Etc.
基本的な考え方 : 同一祖先由来と思われるタンパク質を収集し 置換率を考える 例 : ロイシンジッパーを含むタンパク質 ACA2_YEAST P4535 KRARLLERNRIAASKCRQRKKVAQLQLQKEFNEIKDENRIL AP1_HUMAN P5412 KAERKRMRNRIAASKCRKRKLERIARLEEKVKTLKAQNSEL K K K L K E 4 箇 1 箇 1 箇 置換率に基づいてスコア体系を決めより正確には アミノ酸 aとbがアラインメントされたときのスコアは aとbの組み合わせが同一祖先からの置換によって生じる確率 log aとbの組み合わせがランダム抽出で得られる確率
標準的に使用されているスコア体系 A R N D C Q E G H I L K M F P S T W Y V B Z X * A 4-1 -2-2 -1-1 -2-1 -1-1 -1-2 -1 1-3 -2-2 -1-4 R -1 5-2 -3 1-2 -3-2 2-1 -3-2 -1-1 -3-2 -3-1 -1-4 N -2 6 1-3 1-3 -3-2 -3-2 1-4 -2-3 3-1 -4 D -2-2 1 6-3 2-1 -1-3 -4-1 -3-3 -1-1 -4-3 -3 4 1-1 -4 C -3-3 -3 9-3 -4-3 -3-1 -1-3 -1-2 -3-1 -1-2 -2-1 -3-3 -2-4 Q -1 1-3 5 2-2 -3-2 1-3 -1-1 -2-1 -2 3-1 -4 E -1 2-4 2 5-2 -3-3 1-2 -3-1 -1-3 -2-2 1 4-1 -4 G -2-1 -3-2 -2 6-2 -4-4 -2-3 -3-2 -2-2 -3-3 -1-2 -1-4 H -2 1-1 -3-2 8-3 -3-1 -2-1 -2-1 -2-2 2-3 -1-4 I -1-3 -3-3 -1-3 -3-4 -3 4 2-3 1-3 -2-1 -3-1 3-3 -3-1 -4 L -1-2 -3-4 -1-2 -3-4 -3 2 4-2 2-3 -2-1 -2-1 1-4 -3-1 -4 K -1 2-1 -3 1 1-2 -1-3 -2 5-1 -3-1 -1-3 -2-2 1-1 -4 M -1-1 -2-3 -1-2 -3-2 1 2-1 5-2 -1-1 -1-1 1-3 -1-1 -4 F -2-3 -3-3 -2-3 -3-3 -1-3 6-4 -2-2 1 3-1 -3-3 -1-4 P -1-2 -2-1 -3-1 -1-2 -2-3 -3-1 -2-4 7-1 -1-4 -3-2 -2-1 -2-4 S 1-1 1-1 -1-2 -2-1 -2-1 4 1-3 -2-2 -4 T -1-1 -1-1 -1-2 -2-1 -1-1 -1-2 -1 1 5-2 -2-1 -1-4 W -3-3 -4-4 -2-2 -3-2 -2-3 -2-3 -1 1-4 -3-2 11 2-3 -4-3 -2-4 Y -2-2 -2-3 -2-1 -2-3 2-1 -1-2 -1 3-3 -2-2 2 7-1 -3-2 -1-4 V -3-3 -3-1 -2-2 -3-3 3 1-2 1-1 -2-2 -3-1 4-3 -2-1 -4 B -2-1 3 4-3 1-1 -3-4 -3-3 -2-1 -4-3 -3 4 1-1 -4 Z -1 1-3 3 4-2 -3-3 1-1 -3-1 -1-3 -2-2 1 4-1 -4 X -1-1 -1-2 -1-1 -1-1 -1-1 -1-1 -1-2 -2-1 -1-1 -1-1 -4 * -4-4 -4-4 -4-4 -4-4 -4-4 -4-4 -4-4 -4-4 -4-4 -4-4 -4-4 -4 1 BLOSUM62 PAM Dayhoff et al. 1978 タンパク質の全長の類似性に基づいて導出された BLOSUM Henikoff and Henikoff, 1992 タンパク質の高い保存性を持つ部分配列に基づいて導出された BLAST など有名なツールにも使用されている
より根拠のあるアラインメントのスコア体系 (2) 塩基の挿入 欠損はまとめて起こりやすい CCCAAAAAAAAAAGGG CCCAAAAAGGG CCCAAAAAAAAAAGGG CCCA-A-A-A-A-GGG CCCAAAAAAAAAAGGG CCCAAAAA-----GGG どちらが良いアラインメント?
より複雑な GAP のペナルティ GAP open penalty GAP の始まりに対して付けられるペナルティー (ex. -5) GAP extension penalty 連続して GAP が入っている領域の 2 番目以降の GAP に対して付けられるペナルティー (ex. -1) CCCAAAAAAAAAAGGG CCCA-A-A-A-A-GGG -5-5-5-5-5 CCCAAAAAAAAAAGGG CCCAAAAA-----GGG - -5-25 -9 GAP ペナルティーが一定でない 再帰式の変更が必要
もはや 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?
最大スコアの場合分け 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)
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
アラインメントの準最適解を求めるもの BLAST (Altschul SF et al. 199) 基本的に GAP なし高速 FASTA (Lipman DJ et al. 1985) GAP を扱うことができるそこそこ速い パッケージ中の ssearch は最適解を求める BLAT (Kent J et al. 22) BLAST より精度が落ちるとされるが とても高速
BLAST word の抽出 Query を word(n アミノ酸単位 ) に分割 データベース中の word に一致する部分の検出 Hit からの拡張 Hit 部分の統合
マルチプルアラインメント 2 本の配列のアラインメントをペアワイズアラインメントという 3 本以上の配列のアラインメントをマルチプルアラインメントという ツールとしては CLUSTALW(Thompson JD et al. 1994) などが有名 AC-GT-A ACCAT-A GC--T-A GC--TCA
マルチプルアラインメントを行う手法 ツリーベース法 ACGTA ACCATA GCTA GCTCA AC-GTA ACCATA GCT-A GCTCA AC-GT-A ACCAT-A GC--T-A GC--TCA シミュレーテッドアニーリング法 隠れマルコフモデル