第 3 回バイオインフォマティクスアルゴリズム アラインメントアルゴリズム 慶應義塾大学先端生命科学研究所
基本的なアルゴリズム設計技法 再帰法 分割統治法 動的計画法
分割統治法とは? (1) 例 :(18, 37, 21, 14, 7, 12, 19, 6) を小さいものから順に並べ替えることを考える 最小値 2 番目に小さい値 3 番目に小さい値 を求めていくやり方では数列の長さを n とすると 計算量は O(n 2 )
分割統治法とは? (2) この問題を 2 つの部分に分割 : (18, 37, 21, 14) を並べ替え それとは別に (7, 12, 19, 6) も並べ替える (14, 18, 21, 37) と (6, 7, 12, 19) に並べ替えられた これを併合して (6, 7, 12, 14, 18, 19, 21, 37) 計算量は O(nlog n)
分割統治法とは? (3) n 2 log n (18, 37, 21, 14, 7, 12, 19, 6) (18, 37, 21, 14), (7, 12, 19, 6) (18, 37), (21, 14), (7, 12), (19, 6) (18), (37), (21), (14), (7), (12), (19), (6) (18, 37), (14, 21), (7, 12), (6, 19) (14, 18, 21, 37), (6, 7, 12, 19) (6, 7, 12, 14, 18, 19, 21, 37)
分割統治法とは? (4) 再帰式で表すと 関数 msort( 数列 ): 数列が空なら それをそのまま返す数列中の数が 1 つなら それをそのまま返す msort_left = msort( 前半部分の数列 ) msort_right = msort( 後半部分の数列 ) merge(msort_left, msort_right) を返 最初からまともに全部並べようとするより 問題を分割して解き 後で統合した方が早く解けるケースがる
途中計算の記録による無駄な計算の省略 ( フィボナッチ数列 ) fb(5) 1 fb(5) 1 fb(4) 2 fb(3) 7 fb(4) 2 fb(3) 7 fb(3) 3 fb(2) 6 fb(2) 8 fb(1) 9 fb(3) 3 fb(2) 6 fb(2) 4 fb(1) 5 fb(2) 4 fb(1) 5 (Before) (After)
動的計画法 分割統治法による問題の分割 途中計算の記録による高速化
アラインメントアルゴリズム
入手可能なゲノム (DNA) 配列の例 種 ゲノムサイズ 遺伝子数 大腸菌 4.6M 4,000 酵母菌 15M 6,000 線虫 100M 14,000 ショウジョウバエ 170M 12,000 ヒト 3,000M 25,000
ゲノム配列からまず知りたいこと 遺伝子領域はどこ? 機能は? 実験で確かめるのは時間と労力と費用がかか
相同性検索 既に実験的に確かめられた配列からマッチするものを探す 機能予測にも広く使われている 入力配 類似配 MLSAQERAAI MLSFVERATL データベー
配列は完全にマッチするとは限らない 塩基置 ATCAATCGATCGATC ATCCATTGAACCATC ATCAATCGATCGATC 挿 ATCAATCGAGGAGGATCGATC ATCAATCGATCGATC 欠 ATCAACGATC 配列の置換をどのように処理するか配列の挿入 欠損をどのように処理するか
アラインメント 置換 挿入 欠損を考慮して塩基配列あるいは アミノ酸配列の似た部分をそろえることギャップ - を挿入する CAAGACATTTTAC CATACACTTTAC CA-AGACATTTTAC CATACAC--TTTAC ** * ** *****
アラインメントの使用用途 遺伝子領域予測 遺伝子の機能予測 保存部位 モチーフ推定 系統樹作成
全体と一部のアラインメント グローバルアラインメント全体を揃えるようなアラインメント cggtc----acgtcgagtg cg-tacgt acgacgg--- ** * *** ** ローカルアラインメント一部をなるべく揃えるようなアラインメント cggtcacgtcgagtg cgtacgtacgtcgg ****
ヒトと大腸菌は同一祖先を持つ!? 1810 1820 1830 1840 1850 1860 Human GCCCGTCGCTACTACCGATTGGA-TGGTTTAGTGAGGCCCTCGGATCGGCCCCGCCGGGG ::::::: : ::: :: ::: ::: :: : : : :: Ecoli GCCCGTCAC----ACC-ATGGGAGTGGGTTGCAAAAGAAGTAGG---------------- 1410 1420 1430 1870 1880 1890 1900 1910 1920 Human TCGGCCCACGGCCTGGCGGAGCGCTGAGAAGACGGTCGAACTTGACTATCTAGAGGAAGT : : :: : :: :: ::::: : : :: : ::::: : :: Ecoli TAGCTTAACCTTCGGGAGG-GCGCTTACCACTTTGTGATTCATGACT-----GGGGTG-- 1440 1450 1460 1470 1480 1490 1930 1940 1950 1960 Human AAAAGTCGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTA ::::::::::::::: ::::::: :::::::::: :::::: Ecoli --AAGTCGTAACAAGGTAACCGTAGGGGAACCTGCGGTTGGATCACCTCCTTA 1500 1510 1520 1530 1540 ヒトの 18S rrna と大腸菌の 16S rrna のアラインメン
アラインメントを定式化する 良いアラインメント = スコアが高いアラインメントとすれば どのようなスコア体系にすればよいか? そのスコア体系の中で最大のスコアおよび対応するアラインメントをどのように求めればよいか?
単純なスコア体系 (1) 塩基がマッチしたときのスコア ( > 0 ) 塩基がマッチしないときのスコア ( < 0 ) ギャップのペナルティー ( < 0 )
単純なスコア体系 (2) 各位置のスコアを合計するだけ ( スコアの独立性 : 各位置のスコアは その位置だけで決まる ) ** x** * at-catgc atgaat-c マッチスコア 5 + ミスマッチペナルティー 1 + ギャップペナルティー 2
最適なスコアを計算するには? アラインメントの表現方法を決める 最大スコアの数学的な定義 ( 再帰的定義 ) 定義式の効率の良い計算方法
アラインメントはグラフで表現できる seq2 A A C G A seq1 C G A-CG AACG
演習問題 (1) 以下の経路に対応するアラインメントを示せ T A C G C A C T G
アラインメントの計算量 (1) 全通りの経路をしらみつぶしに調べれば最適解が求まる アラインメントをかける配列の長さ n が長くなると計算時間が急速に増える
アラインメントの計算量アラインメントの計算量 (2) (2) [ ] = + n k k n k k n n n 1 2 2 )! (! )! (2!) ( )! (2 アラインメントの場合の数は AT-G A-CG A-TG AC-G とを同一とみなすなら n n n n π 2 2 2!) ( )! (2
計算量の改善 アラインメント最適スコアの再帰的定義 動的計画法の適用
最大スコアの数学的な定義を考える 再帰的定義的定義の復習 (1) 階乗の定義 1 n! = n x (n-1) x (n-2) x x 2 x 1 の部分の意味が不明瞭
最大スコアの数学的な定義を考える 再帰的定義的定義の復習 (2) 階乗の定義 2 n! = 1 n ( n 1)! for n 1 for n > 1 (1) (2) 階乗の定義の中に階乗がさらに登場することに注目 3! を求めると (2) より 3!=3x(3-1)! = 3x2! (2) より 3x2!=3x2x1! (1) より 3x2x1!=3x2x1=6
直前のノード MS(i,j) を位置 (i,j) における最大スコアと定義 MS(i-1,j-1) MS(i-1,j) +SC(seq1(i-1),seq2(j-1)) +GP MS(i,j-1) +GP MS(i,j) 但し seq(i) は配列 seqのi 番目の塩基 SC(a,b) は塩基 a,bがマッチしたときのスコア GPは一定 (i,j) を固定すれば SC(seq1(i-1),seq2(j-1)) も一定
最大スコアの最大スコアの再帰的再帰的定義定義 + + + = = = 1)) 1),seq2( (seq1( 1) 1, ( ) 1, ( 1), ( max 0 0 or 0 if ), ( j i SC j i MS GP j i MS GP j i MS j i j i MS
本当に最大のスコア? MS(i,j) を位置 (i,j) における最大スコアと定義 GP は一定 (i,j) を固定すれば SC(seq1(i),seq2(i)) も一定 MS(i-1,j-1) MS(i-1,j) MS( i, j 1) + GP MS = max MS( i 1, j) + GP MS( i 1, j 1) + SC(seq1( i 1),seq2( j 1)) +SC(seq1(i-1),seq2(j-1)) +GP AS = AS -1 + GP AS -1 MS(i, j- MS(i,j-1) AS -1 +GP MS AS MS(i,j) AS = AS -1 + GP MS(i, j-1) + MS
#include <stdio.h> #define H 0 /* 横方向 */ #define V 1 /* 縦方向 */ #define D 2 /* 斜め方向 */ #define NDIR 3 /* 方向数 */ #define Gap_Penalty -5 /* ギャップペナルティ */ static char seq1[max_seq_len], seq2[max_seq_len]; /* 塩基配列 */ int find_max_score(int i, int j){ main(){ strcpy(seq1, "acg"); strcpy(seq2, "aacg"); printf("%d n", find_max_score(3, 4)); } if(i == 0){ max_score = 0; max_dir = H;} else if(j == 0){ max_score = 0; max_dir = V; } else { score_tmp[v] = find_max_score(i - 1, j) + Gap_Penalty; /* (i-1,j) までの最適スコアにギャップペナルティを足したもの */ score_tmp[h] = find_max_score(i, j - 1) + Gap_Penalty; /* (i,j-1) までの最適スコアにギャップペナルティを足したもの */ score_tmp[d] = find_max_score(i - 1, j - 1) + score(seq1[i - 1], seq2[j - 1]); /* (i-1,j-1) までの最適スコアにアラインメントスコアを足したもの */ } max_dir = find_max_elem(score_tmp, NDIR); /* score_tmp[v], score_tmp[h], score_tmp[d] のうちで最大の方向 (V,H,D) を求める */ max_score = score_tmp[max_dir]; /* 上記で求めた方向に対応する最大スコアを求める */ } return max_score; /* 最大スコアを返す */ C 言
2 つの関数を定義しよう int score(char a, char b); 入力された 2 つの文字が同一ならマッチスコア (+10) 異なるならミスマッチスコア (-7) を与える int find_max_elem(int score_tmp[], int ndir); 長さ ndir( 実際は縦 横 斜で 3 つ ) の配列 score_tmp[] の中で最大の要素が何番目にあるかを返す C 言
各ノードの情報を表す構造体 MS(i,j) struct { int score; /* そのノードに至ったときの最大スコア */ int direction; /* そのノードに到達して最大スコアになったときの直前のノードの方角 */ } score_path_matrix[max_seq_len][max_seq_len]; /* MAX_SEQ_LEN は扱える最大の配列長さ */ C 言
#!/usr/bin/perl w use strict; use vars qw($seq1 $seq2); *::Gap_Penalty = -5; sub find_max_score($$){ my($i, $j) = @_; my %score_tmp; my $max_dir; my $max_score; $seq1 = "acg"; $seq2 = "atcg"; my $seq1_len = length($seq1); my $seq2_len = length($seq2); printf("%d n", find_max_score($seq1_len, $seq2_len)); } if($i == 0){ $max_score = 0; $max_dir = "H"; } elsif($j == 0){ $max_score = 0; $max_dir = "V"; } else { $score_tmp{"v"} = find_max_score($i - 1, $j) + $::Gap_Penalty; $score_tmp{"h"} = find_max_score($i, $j - 1) + $::Gap_Penalty; $score_tmp{"d"} = find_max_score($i - 1, $j - 1) + score(substr($seq1, $i-1, 1), substr($seq2, $j-1, 1)); $max_dir = find_max_elem(%score_tmp); $max_score = $score_tmp{$max_dir}; } return $max_score; Perl
2 つの関数を定義しよう sub score($$); 入力された 2 つの文字が同一ならマッチスコア (+10) 異なるならミスマッチスコア (-7) を与える sub find_max_elem(%); $_{ V },$_{ H },$_{ D } を比較し もっとも大きい値を持つキー ( V, H または D ) を返す Perl