bioinfo-a10s-4_align

Size: px

Start display at page:

Download "bioinfo-a10s-4_align"

かおりあざみ
5 years ago
Views:

1 4 (2)

2 ATCAATCGATCGATC ATCCATTGAACCATC ATCAATCGATCGATC ATCAATCGAGGAGGATCGATC ATCAATCGATCGATC ATCAACGATC

3 - CAAGACATTTTAC CATACACTTTAC CA-AGACATTTTAC CATACAC--TTTAC ** * ** *****

4 seq2 A A C G A seq1 C G A-CG AACG

5 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) GP (i,j) SC(seq1(i-1),seq2(i-1))

6 + + + = = = 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

7 A T Match = +3 Mismatch = -3 GAP = -2 G A MS(i,j)

8 (i,j) (1) MS(i,j) (2) direction(i,j) MS(i,j)

9 seq A T C G A seq C G

10 seq A T C G seq A C a_seq1: GCAA a_seq2: GC-A G a_seq1: AACG a_seq2: A-CG

11 (2) atcga atga (1)MS(i,j) (2)direction(i,j)

12 (Perl) (i,j) (1) MS(i,j) (2) direction(i,j) MS(i,j) $score_path_matrix[$i]->[$j]->{"score" = $max_score; $score_path_matrix[$i]->[$j]->{"direction" = $max_dir;

13 #!/usr/bin/perl w use strict; use vars qw($seq1 *::Gap_Penalty = -5; sub find_max_score($$){ my($i, $j) my %score_tmp; my $max_dir; my $max_score; if($i == ){ $max_score = ; $max_dir = "H"; elsif($j == ){ $max_score = ; $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; $score_path_matrix[$i]->[$j]->{"score" = $max_score; $score_path_matrix[$i]->[$j]->{"direction" = $max_dir; return $max_score; Perl

14 (C (C ) MS(i,j) struct { int score; /* */ int direction; /* */ score_path_matrix[max_seq_len][max_seq_len]; /* MAX_SEQ_LEN */ (i,j) (1) MS(i,j) (2) direction(i,j) score_path_matrix[i][j].score = max_score; score_path_matrix[i][j].direction = max_dir;

15 #include <stdio.h> #define H /* */ #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){ if(i == ){ max_score = ; max_dir = H; else if(j == ){ max_score = ; 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]; /* */ score_path_matrix[i][j].score = max_score; score_path_matrix[i][j].direction = max_dir; return max_score; /* */ C

16 seq A T C G A a_seq1: a_seq2: G G C C - T A A seq C G a_seq1: A-CG a_seq2: ATCG

17 void board_to_alignment(char a_seq1[], char a_seq2[], int m, int n){ int i, j, 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:??????????????????? for(i =, j = p - 1;i < j;i ++, j--){ tmp = a_seq1[i]; a_seq1[i] = a_seq1[j]; a_seq1[j] = tmp; tmp = a_seq2[i]; a_seq2[i] = a_seq2[j]; a_seq2[j] = tmp; a_seq1[p] = ' '; a_seq2[p] = ' '; C

18 (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);

19 sub board_to_alignment($$){ my($m, $n) my($i, $j, $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

20 (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);

21 n O(2 2n /n 1/2 ) O(n 2 ) O(n 2 )

ﾀｲﾄﾙ

ﾀｲﾄﾙ第 3 回バイオインフォマティクスアルゴリズムアラインメントアルゴリズム慶應義塾大学先端生命科学研究所基本的なアルゴリズム設計技法再帰法分割統治法動的計画法分割統治法とは? (1) 例 :(18, 37, 21, 14, 7, 12, 19, 6) を小さいものから順に並べ替えることを考える最小値 2 番目に小さい値 3 番目に小さい値を求めていくやり方では数列の長さを n とすると