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

(1)

第

₅

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

アラインメントアルゴリズム

₍₃₎

(2)

アラインメント

CA-AGACATTTTAC

CATACAC--TTTAC

** * ***

CAAGACATTTTAC

CATACACTTTAC

置換、挿入、欠損を考慮して塩基配列あるいは、

アミノ酸配列の似た部分をそろえること

ギャップ”-”を挿入する

(3)

A

C

G

A

C

G

A-CG

AACG

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

seq1

seq2

(4)

直前のノード

MS(i,j)

MS(i,j-1)

MS(i-1,j)

MS(i-1,j-1)

+GP

+SC(seq1(i),seq2(j))

+GP

MS(i,j)を位置(i,j)における最大スコアと定義

(5)

最大スコアの再帰的定義

⎪

⎩

⎪⎪

⎨

⎧

⎪

⎩

⎪

⎨

⎧

−

+

−

+

−

+

−

=

))

1 (

2 ),

1 (

)

1 ,

1 (

)

,

1 (

)

1 ,

(

max

0 or

0 if

0 )

,

(

j

seq

i

seq

SC

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)

最大スコア方向の情報からアラインメントを求める

a_seq1: A-CG

a_seq2: ATCG

G

A

T

C

_G

A

C

seq１

seq２

0

3

1 -1

_-2

₂

7

0

1 -1

-2

0

4

2 a_seq1:

a_seq2:

G

C

-T

A

(8)

void board_to_alignment(char a_seq1[], char a_seq2[], int m, int n){

**/* a_seqにアラインメント結果を格納する。m, nは２つの入力配列のそれぞれの長さ */**

int i, j, p;

**/* i,jが現在注目対象となっている２つの配列の塩基の位置、pがアラインメント格納位置**

char tmp;

for(p = 0, i = m, j = n; i != 0 || j != 0;){

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] = '¥0';

a_seq2[p] = '¥0';

strrev(a_seq1);

strrev(a_seq2);

}

C言

(9)

アラインメントの呼び出し

_(C

言語

_言語

₎

main(){

static char a_seq1[100], a_seq2[100];

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は２つの入力配列のそれぞれの長さ */**

**my(@a_seq1, @a_seq2); /* アラインメント結果を格納する配列 */**

my($i, $j, $p);

**/* $i, $jが現在注目対象となっている２つの配列の塩基の位置、$pがアラインメント格納位置 */**

my($tmp);

for($p = 0, $i = $m, $j = $n;

$i != 0 || $j != 0;){

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

}

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

0 1

a_seq1: G

C

a_seq2: G

C

G

A

T

C

_G

A

C

seq１

seq２

0

3

1 -1

_-2

₂

7

0

1 -1

-2

0

4

2 0 1 2 3

seq1: A C G

seq2: A T C G

i = 1

p = 1

j = 2

(13)

変数の中味の具体例 (2)

0 1 2

a_seq1: G C

-a_seq2: G C

T

G

A

T

C

_G

A

C

seq１

seq２

0

3

1 -1

_-2

₂

7

0

1 -1

-2

0

4

2 0 1 2 3

seq1: A C G

seq2: A T C G

i = 1

p = 2

j = 1

(14)

結局どこが動的計画法だったのか？

z

アラインメントの問題を、各位置について、以下

のように分割

– 配列１と２の塩基をマッチさせる場合

– 配列１にGAPを入れる場合

– 配列２にGAPを入れる場合

z

各ノードに以下の途中経過を記録し、高速化

– そのノードまでの最大スコア

– (そのノードで最大スコアになるときの直前のノード)

(15)

ローカルアラインメント

₍₁₎

cagtagctgatcgatgctagctgat

|||| ||||||||||

tttatgat—gatgctagctacactac

配列１

配列２

本当に似ている部

切り捨てたい部

アラインメントを開始

したいポイント

(16)

ローカルアラインメント

₍₂₎

G

A

T

C

_G

A

C

seq１

seq２

0

3

1 -1

-2

2

7

0

1 -1

-2

0

4

2 負のスコアに未来はない？

(17)

ローカルアラインメント

₍₃₎

⎪

⎩

⎪

⎨

⎧

⎪

⎩

⎪

⎨

⎧

−

+

−

+

−

+

=

))

1 (

2 ),

1 (

)

1 ,

1 (

)

,

1 (

)

1 ,

(

0 max

0 or

0 if

0 )

,

(

j

seq

i

seq

SC

j

i

MS

GP

j

i

MS

GP

j

i

MS

j

i

j

i

MS

(18)

最大スコアと方角の記録

₍₁₎

G

T

A

G

_C

C

A

seq１

seq２

0

3

1

0

1

6

4

0

(19)

最大スコアと方角の記録

₍₂₎

G

T

A

G

_C

C

A

seq１

seq２

0

3

1

0

1

6

4

0 a_seq1: GA

a_seq2: GA

a_seq1: AG

a_seq2: AG

ローカルアラインメントはMS(i,j）が最も高い位置から始める

(20)

より根拠のあるアラインメントのスコア体系

z

ゲノム配列には生物学的にどのような変異が生じやすいかを考えて

アラインメントのスコア体系を構築すべき

z

DNAでどのような変異が起こりやすいか

– スペーサ

– 繰り返し配列

– プロモータ

– タンパク質のコード領域(コドンの第1～3塩基)

– Etc.

z

タンパク質の配列でどのような変異が起こりやすいかの方が考慮す

べき要素が減り、考えやすい

– 疎水性のアミノ酸は他の疎水性のアミノ酸に変化しやすい

– 正電荷を持つアミノ酸は他の正電荷を持つアミノ酸に変化しや

い

– Etc.

(21)

基本的な考え方：

同一祖先由来と思われるタンパク質を収集し、置換率を考える

例：ロイシンジッパーを含むタンパク質

ACA2_YEAST|P40535 KRARLLERNRIAASKCRQRKKVAQLQLQKEFNEIKDENRIL

AP1_HUMAN|P05412 KAERKRMRNRIAASKCRKRKLERIARLEEKVKTLKAQNSEL

K → K

4箇

K → L

1箇

K → E

1箇

置換率に基づいてスコア体系を決め

より正確には、アミノ酸aとbがアラインメントされたときのスコアは、

ム抽出で得られる確率

の組み合わせがランダ

と

生じる確率

先からの置換によって

の組み合わせが同一祖

と

b

a

b

a

log

(22)

標準的に使用されているスコア体系

z

PAM

– Dayhoff et al. 1978

– タンパク質の全長の類似性

に基づいて導出された

z

BLOSUM

– Henikoff and Henikoff,

1992

– タンパク質の高い保存性を

持つ部分配列に基づいて導

出された

– BLASTなど有名なツールに

も使用されている

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 0 -1 -1 0 -2 -1 -1 -1 -1 -2 -1 1 0 -3 -2 0 -2 -1 0 -4 R -1 5 0 -2 -3 1 0 -2 0 -3 -2 2 -1 -3 -2 -1 -1 -3 -2 -3 -1 0 -1 -4 N -2 0 6 1 -3 0 0 0 1 -3 -3 0 -2 -3 -2 1 0 -4 -2 -3 3 0 -1 -4 D -2 -2 1 6 -3 0 2 -1 -1 -3 -4 -1 -3 -3 -1 0 -1 -4 -3 -3 4 1 -1 -4 C 0 -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 0 0 -3 5 2 -2 0 -3 -2 1 0 -3 -1 0 -1 -2 -1 -2 0 3 -1 -4 E -1 0 0 2 -4 2 5 -2 0 -3 -3 1 -2 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4 G 0 -2 0 -1 -3 -2 -2 6 -2 -4 -4 -2 -3 -3 -2 0 -2 -2 -3 -3 -1 -2 -1 -4 H -2 0 1 -1 -3 0 0 -2 8 -3 -3 -1 -2 -1 -2 -1 -2 -2 2 -3 0 0 -1 -4 I -1 -3 -3 -3 -1 -3 -3 -4 -3 4 2 -3 1 0 -3 -2 -1 -3 -1 3 -3 -3 -1 -4 L -1 -2 -3 -4 -1 -2 -3 -4 -3 2 4 -2 2 0 -3 -2 -1 -2 -1 1 -4 -3 -1 -4 K -1 2 0 -1 -3 1 1 -2 -1 -3 -2 5 -1 -3 -1 0 -1 -3 -2 -2 0 1 -1 -4 M -1 -1 -2 -3 -1 0 -2 -3 -2 1 2 -1 5 0 -2 -1 -1 -1 -1 1 -3 -1 -1 -4 F -2 -3 -3 -3 -2 -3 -3 -3 -1 0 0 -3 0 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 0 -1 0 0 0 -1 -2 -2 0 -1 -2 -1 4 1 -3 -2 -2 0 0 0 -4 T 0 -1 0 -1 -1 -1 -1 -2 -2 -1 -1 -1 -1 -2 -1 1 5 -2 -2 0 -1 -1 0 -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 0 -3 -3 -3 -1 -2 -2 -3 -3 3 1 -2 1 -1 -2 -2 0 -3 -1 4 -3 -2 -1 -4 B -2 -1 3 4 -3 0 1 -1 0 -3 -4 0 -3 -3 -2 0 -1 -4 -3 -3 4 1 -1 -4 Z -1 0 0 1 -3 3 4 -2 0 -3 -3 1 -1 -3 -1 0 -1 -3 -2 -2 1 4 -1 -4 X 0 -1 -1 -1 -2 -1 -1 -1 -1 -1 -1 -1 -1 -1 -2 0 0 -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

(23)

CCCAAAAAAAAAAGGG

CCCAAAAA---GGG

CCCAAAAAAAAAAGGG

CCCA-A-A-A-A-GGG

どちらが良いアラインメント？

CCCAAAAAAAAAAGGG

CCCAAAAAGGG

塩基の挿入、欠損はまとめて起こりやすい

より根拠のあるアラインメントのスコア体系

₍₂₎

(24)

より複雑な

_GAP

のペナルティ

_{のペナルティ}

z

GAP open penalty –GAPの始まりに対して付けられる

ペナルティー (ex. -5)

z

GAP extension penalty –連続してGAPが入っている

領域の２番目以降のGAPに対して付けられるペナル

ティー (ex. -1)

CCCAAAAAAAAAAGGG

CCCAAAAA---GGG

|

-5

CCCAAAAAAAAAAGGG

CCCA-A-A-A-A-GGG

| | | | |

-5-5-5-5-5

ＧＡＰペナルティーが一定でない…再帰式の変更が必要

--25

-9

(25)

もはやGAPペナルティーは一定ではない…

MS(i,j)

MS(i,j-1)

MS(i-1,j)

MS(i-1,j-1)

+SC(seq1(i),seq2(j))

+GP

_open

?,+GP

_ext

?

+GP

_open

?,+GP

_ext

?

MS(i,j)を位置(i,j)における最大スコアと定義

(26)

最大スコアの場合分け

z

MS(i, j)を位置(i, j)における最

大スコアと定義

z

MS

_H

(i, j)を直前のノードが位

置(i, j-1)だったときの位置(i, j)

における最大スコアと定義

MS(i,j)

MS

_V

(i,j)

MS

_H

(i,j)

MS

_D

(i,j)

⎪

⎩

⎪

⎨

⎧

=

)

,

(

)

,

(

)

,

(

)

,

(

j

i

MS

j

i

MS

j

i

MS

j

i

MS

D

V

H

(27)

２つ前まで見れば、経路のペナルティーは一定

MS

_H

(i,j)

MS

_H

(i,j-1) + GP

_ext

MS

_D

(i,j-1)

+GP

_open

(i-1,j-2)

(i,j-2)

⎩

⎨

⎧

+

−

+

−

=

ext

open

)

1 ,

(

)

1 ,

(

max

GP

j

i

MS

GP

j

i

MS

H

D

H

(28)

アラインメントの準最適解を求めるもの

z

BLAST (Altschul SF et al. 1990)

基本的にGAPなし

高速

z

FASTA (Lipman DJ et al. 1985)

GAPを扱うことができる

そこそこ速い

※パッケージ中のssearchは最適解を求める

z

BLAT (Kent J et al. 2002)

(29)

BLAST

z

wordの抽出

– Queryをword(nアミノ酸単位）に分割

– データベース中のwordに一致する部分の検出

z

Hitからの拡張

z

Hit部分の統合

(30)

マルチプルアラインメント

z

2本の配列のアラインメントをペアワイズアラ

インメントという

z

3本以上の配列のアラインメントをマルチプル

アラインメントという

z

ツールとしては、CLUSTALW(Thompson

JD et al. 1994)などが有名

AC-GT-A

ACCAT-A

GC--T-A

GC--TCA

(31)

マルチプルアラインメントを行う手法

z

ツリーベース法

z

シミュレーテッドアニーリング法

z

隠れマルコフモデル

ACGTA

ACCATA

AC-GTA

ACCATA

AC-GT-A

ACCAT-A

GC--T-A

GC--TCA

GCTA

GCTCA

GCT-A

GCTCA

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

第

第

5

5

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

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

アラインメントアルゴリズム

アラインメントアルゴリズム

(3)

(3)

アラインメント

アラインメント

CA-AGACATTTTAC

CATACAC--TTTAC

** * ** *****

CAAGACATTTTAC

CATACACTTTAC

置換、挿入、欠損を考慮して塩基配列あるいは、

アミノ酸配列の似た部分をそろえること

ギャップ”-”を挿入する

A

A

C

G

A

C

G

A-CG

AACG

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

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

seq1

seq2

直前のノード

直前のノード

MS(i,j)

MS(i,j-1)

MS(i-1,j)

MS(i-1,j-1)

+GP

+SC(seq1(i),seq2(j))

+GP

MS(i,j)を位置(i,j)における最大スコアと定義

最大スコアの再帰的定義

最大スコアの再帰的定義

⎪

⎪

⎩

⎪⎪

⎨

⎧

⎪

⎩

⎪

⎨

⎧

−

−

+

−

−

+

−

+

−

=

=

=

))

1

(

2

),

1

(

1

(

)

1

₅

₅

₍₃₎

₍₃₎

** * ***

_(C

_(C

_言語

₎

₎

**/* そのノードに至ったときの最大スコア*/**

**/* そのノードに到達して最大スコアになったときの**

**直前のノードの方角*/**

**/* MAX_SEQ_LEN は扱える最大の配列長さ*/**