ゲノム統合データベースからの知識発見

配列解析アルゴリズム特論
配列アライメントI
阿久津達也
京都大学化学研究所
バイオインフォマティクスセンター

講義予定

１１月５日配列アライメントI（担当：阿久津）
講義内容
Smith-waterman アルゴリズム




１１月１２日配列アライメントII（担当：山口）



糖尿病モデル動物等での解析例等を用いた解説
１１月２６日機械モデルと学習手法（担当：馬見
塚）


マルチプルアライメント，E-Value
PSI-BLAST
１１月１９日 QTL解析（担当：中谷）


Hash法を用いた高速検索
Gibbsサンプリング
HMM、有限混合モデル、EMアルゴリズムなど
１２月３日文字列照合アルゴリズム（担当：坂内）

文字列照合アルゴリズム、索引構造、最適パターン探索
配列アライメント




バイオインフォマティクスの
最重要技術の一つ
２個もしくは３個以上の配列
の類似性の判定に利用
文字間の最適な対応関係を
求める（最適化問題）
配列長を同じにするように、
ギャップ記号（挿入、欠失に
対応）を挿入
A L G F G S L Y G
A L G G V S V G
A L G F G
A L G
S L Y G
G V S V
G
スコア行列（置換行列）

残基間（アミノ酸文字間）の類似性を表す行列

PAM250, BLOSUM45 など
A
A
R
N
D
C
Q
E
G
H
I
L
K
M
F
P
S
T
W
Y
V
5
-2
-1
-2
-1
-1
T
W
Y
V
-2 -1 -2 -1 -1 -1 0 -2 -1 -2 -1 -1 -3 -1 1 0
7 -1 -2 -4 1 0 -3 0 -4 -3 3 -2 -3 -3 -1 -1
-1 7 2 -2 0 0 0 1 -3 -4 0 -2 -4 -2 1 0
-2 2 8 -4 0 2 -1 -1 -4 -4 -1 -4 -5 -1 0 -1
-4 -2 -4 13 -3 -3 -3 -3 -2 -2 -3 -2 -2 -4 -1 -1
1 0 0 -3 7 2 -2 1 -3 -2 2 0 -4 -1 0 -1
R
N
D
C
Q
E
G
H
I
L
K
M
F
P
-3
-3
-4
-5
-5
-1
-2
-1
-2
-3
-3
-1
0
3
-3
-4
-1
-3
BLOSUM50 スコア行列
（置換行列）の一部分
S
スコア行列の導出


基本的には頻度の比の対数をスコアとする
BLOSUM行列




既存のスコア行列を用いて多くの配列のアライメント
を求め、ギャップ無しの領域（ブロック）を集める
残基がL％以上一致しているものを同一クラスタに
集める
同じクラスタ内で残基aが残基bにアラインされる頻
度Aabを計算
qa=∑b Aab / ∑cd Acd, pab=Aab / ∑cd Acd を求め、
ｓ
（a,b)=log(pab/qaqb) としたのち、スケーリングし近傍
の整数値に丸める
ペアワイズ・アライメント


配列が２個の場合でも可能なアラインメントの個数は
指数オーダー
しかし、スコア最大となるアライメント（最適アライメン
ト）は動的計画法により、O(mn)時間で計算可能
（m,n:入力配列の長さ）
入力配列
AGCT, ACGCT
アラインメント
AGCT ACGCT
スコア -3
AG - CT
ACGCT
1
最適アラ
インメント
A - GCT
ACGCT
3
- AGC - - T
AC - - GCT
-5
（同じ文字の時1、違う文字の時-1,ギャップ1文字-1）
ギャップペナルティ
ギャップ（g =3)


線形コスト -gd
A L G
L Y G
 g：ギャップ長
A V G V S D L
G
 ｄ：ギャップペナルティ
 この図の例では、ペナルティ= -3d
アフィンギャップコスト –d – e(g-1)
 ｄ：ギャップ開始ペナルティ
 e：ギャップ伸張ペナルティ
 この図の例では、ペナルティ= -d - 2e
 よく利用されるペナルティ (d,e)=(12,2),(11,1)
動的計画法による大域アライメント(1)
(Needleman-Wunschアルゴリズム)



入力文字列から格子状グラフを構成
アライメントと左上から右下へのパスが一対一対応
最長経路＝最適アライメント
G
G
F
5
-5
K
-2
-5
Y
-5
7
D
1
-6
アライメント
-7
-7
GKY
D
G F V D
GK Y D
V
D
-1
-2
-2
-2
1
0
-4
4
-7
-7
-7
-7
スコア
-7
GF V D
GKY D
-7
G
F V D
5 -7 +7
-7 +4 = 2
-7 -7 -1 +0
-7 -7 = -29
-7 -7 -5 -7
-7 -7 -7 = -47
動的計画法による大域アライメント(2)
F(0,0)
G
F(1,0)
=-d
K
F(2,0)
=-2d
=0
F (0, j )   jd , F (i,0)  id
G
F(0,1)
=-d
DP (動的計画法)による
最長経路(スコア)の計算
F(i-1, j-1)
F(i, j-1)
s(K,F)
F
 F (i  1, j  1)  s ( xi , y j )

F (i, j )  max
F (i  1, j )  d

F (i, j  1)  d

-d
F(0,2)
=-2d
F(i-1, j)
-d
F(i, j)
行列からの経路の復元は、
F(m,n)からmaxで＝となっている
F(i,j)を逆にたどることに行う
（トレースバック）
局所アライメント(1)
(Smith-Watermanアルゴリズム)




配列の一部のみ共通部分があることが多い
⇒共通部分のみのアラインメント
x1x2 … xm, y1y2 … yn を入力とする時、スコアが最大とな
る部分列ペア xixi+1 … xk, yjyj+1 … yh を計算
例えば、HEAWGEH と GAWED の場合、
AWGE
A W －E
というアライメントを計算
大域アライメントを繰り返すとO(m3n3)時間
⇒Smith-WatermanアルゴリズムならO(mn)時間
局所アライメント(2)
動的計画法
の式
0
F ( i  1, j  1)  s ( x y )

,
F ( i , j )  max
F ( i  1, j )  d
F ( i , j  1)  d
（最大のF(i,j)から
トレースバック）
局所アライメント(3)

局所アライメントの正当性の証明（下図）

局所アライメントの定義：x1x2 … xm, y1y2 … yn を入力とする時、
スコアが最大となる部分列ペア xixi+1 … xk, yjyj+1 … yh を計算
0
F ( i  1, j  1)  s ( x y )

,
F ( i , j )  max
F ( i  1, j )  d
F ( i , j  1)  d
maxF ( i , j )}
0
0
0
（一部の辺は
省略）
アフィンギャップコストによる
アライメント
F ( i  1, j )  d
Ix ( i , j )  max
Ix ( i  1, j )  e
F ( i , j  1)  d
Iy ( i , j )  max
Iy ( i , j  1)  e
F ( i  1, j  1)  s ( x , y )

F ( i , j )  max
Ix ( i , j )

Iy ( i , j )



三種類の行列を用いる動的
計画法によりO(mn)時間
Smith-Watermanアルゴリ
ズムとの組み合わせが広く
利用されている
Ix (i, j)
Iy (i, j)
F (i, j)
任意ギャップコストによる
アライメント

動的計画法（右式）によ
り、O(n 3 )時間
（ただし、m=O(n)とす
る）

F ( i  1, j  1)  s ( x , y )

F ( i , j )  max  max F ( k , j )   ( i  k )
 0,..., 1
 max F ( i , k )   ( j  k )
 0,... 1
ギャップコストと計算時間の関係
線形

ギャップコストが
凸の時：
O(n 2α(n))時間
ただし、α(n)は
アッカーマン関数
の逆関数（実用的
には定数と同じ）
O(n 2 )時間
2
凸 O(n α(n))時間
アフィン
任意
O(n 2 )時間
3
O(n )時間
ペアワイズ・アライメントの計算量に
関するその他の結果

線形領域アライメント





スコア計算だけなら簡単
トレースバックが難しい
しかし、分割統治法により
O(n)領域が可能(右図)
O(n2)時間の改善は可能か？
⇒O(n2/logn)が可能
s(xi,yj)が疎(O(n)程度)な場合
(Sparse DP)
⇒O(n logn)程度で可能
計算時間
C×mn
C×(mn/2)
C×(mn/4)
配列検索の実用プログラム(1)
O(mn):mは数百だが、nは数ＧＢにもなる
⇒実用的アルゴリズムの開発
 FASTA:短い配列（アミノ酸の場合、1,2文字、DNAの
場合、4-6文字）の完全一致をもとに対角線を検索
し、さらにそれを両側に伸長し、最後にＤＰを利用。
 BLAST:固定長（アミノ酸では3, DNAでは11）の全て
の類似単語のリストを生成し、ある閾値以上の単語
ペアを探し、それをもとに両側に伸長させる。ギャッ
プは入らない。伸長の際に統計的有意性を利用。
 BLAST,PSI-BLASTの詳細はおそらく次回に説明。

配列検索の実用プログラム(2)
FASTA
A
BLAST
Query
・・・ A A F D M F D A D G G ・・・
C
A
T
G
A
C
類似ワード
G
A
MFD MFE MFN
T
MYD MYE MYN
G
・・・
A
Query
T
( ktup=2 )
・・・ A A F D M F D A D G G ・・・
・・・ E A F S M F E K D G D ・・・
Database
BLAT(The Blast-Like Alignment Tool)
[Genome Res. 12:656]



BLAST: データベース配
列を検索
BLAT: 質問配列を検索
データベース配列のイン
デックスをメモリ上に配
置⇒高速化

HIT
HIT
検索
メモリ
ただし、前処理が必要
前処理
データベース配列
BLATの確率解析：１ヒット完全一致

以下の二つを解析



(I)
(I) 相同領域を見逃さない確率
(II) 相同でない領域のワード
がヒットする個数の期待値
K(=5)文字が一致する確率：
この２個のワードが一致しない確率：
M
K
1-M
K
記号の定義






K: ワード長(5-20程度)
M: 相同領域中の塩基の一致
率(>90%)
H: 相同領域の長さ(50-200)
G: データベース配列(ゲノム配
列)の長さ(3Gbase)
Q: 質問配列の長さ
A: 塩基の種類(4)
(II)
K T
相同領域中に一致ワードが無い確率：
(T = H/K = 20/5 = 4ワード)
(1-M )
相同領域中に一致ワードがある確率：
1-(1-M )
質問配列中のワード数：
K T
Q-K+1
ゲノム配列（G/Kワード）
ワードヒット回数の期待値：
（Q-K+1)*(G/K)*(１/A)
K
BLATの確率解析：１ヒット近似一致

各ワード内で１文字まで
のミスマッチを許す



解析法は基本的に１ヒッ
ト完全一致の場合と同じ


(I) 相同領域を見逃さない
確率
(II) 相同でない領域の
ワードがヒットする個数の
期待値
MK、(1/A)K
を右図のよう
に置き換える
計算時間が難点
(I)
確率：
MK
確率：
M K のかわりに
K*(1-M)*M K-1
M K ＋ K*(1-M)*M K-1 を用いる
(II)
(1/A)
K
のかわりに
(1/A) K + K*(1-(1/A)) *(1/A) K-1
を用いる
BLATの確率解析：Nヒット完全一致
(I)
(II)
n=2
1ワードマッチの期待値
F1 = （Q-K+1)*(G/K)*(１/A) K
W/Kワード
n=3
相同領域中でちょうど
ｎ回ヒットする確率：
Pn = p1 n * (1-p1)T-n * TCn
K
( p1 = M )
N回ヒットする確率：
P = PN + PN+1 + ... + PT
1ワードの後、W文字以内に１個以上
のマッチがある確率
K
S = 1- (1-(1/A) )
W/K
ワード間の間隔がW文字以上離れず
に存在するNワードの期待値
FN = F1 * S
（DNA配列解析：K=11,N=2を利用）
N-1
計算式および実データによるBLATの評価
１ヒット完全一致
M
81%
85%
91%
95%
F
１ヒット近似一致
M
２ヒット完全一致
K=9
K=11 K=13
K=13 K=15 K=17
0.833
0.945
0.998
1.000
0.607
0.808
0.981
0.999
0.373
0.594
0.912
0.994
81%
85%
91%
95%
0.880
0.971
1.000
1.000
0.721
0.900
0.996
1.000
0.526
0.767
0.979
1.000
635783 32512
1719
F
68775
4284
267
M
K=9
K=11
81%
85%
91%
95%
0.508
0.762
0.980
1.000
0.220
0.460
0.884
0.993
F
27
0.1
( Q=500, G=3Gbase )
ESTアライメント
K
１ヒット完全
１ヒット近似
２ヒット完全
３ヒット完全
計算時間（sec.)
ヒト/マウス・アライメント
F
14 399
22
0.3
11
0.1
9
0.0
( M=95%, Q=100, P=99% )
K
１ヒット完全
F
7 13,078,962
１ヒット近似 12 275,671
２ヒット完全 6
237,983
３ヒット完全 5
109,707
( M=86%, Q=100, P=99% )
(10,000EST vs. ヒトゲノム配列)
6
7
K
N
10
2
3.9
35.6
680.1
10
3
3.2
21.4
348.2
11
11
2
3
2.4
2.3
8.1
6.5
92.4
61.8
2×10 2×10 2×10
8
PatternHunter

[Bioinformatics, 18:440]
ＰａｔｔｅｒｎHunter



BLASTの精度かつMegaBlastの高速性を主張
連続した文字をseedとせず、飛び飛びの文字
をseedとする
111010010100110111で1の位置のみを考慮
Query
・・・ A C G T A A A A C C C C G G G G T T ・・・
1 1 1 0 1 0 0 1 0 1 0 0 1 1 0 1 1 1
・・・ A C G A A T G A A C A G G G A G T T ・・・
Database
速度向上の理由
Blast
A C G T A C G T
1 1 1
1 1 1
1 1 1
A C G A C C A G
1
p
p2
PatternHunter
A C G
1 0 1
1 0
1
A C G
T
0
1
0
A
A
0
0
1
C
C
1
0
0
C
G T
1
0 1
A G
1
p3
p2
局所マルチプルアライメント
(Local Multiple Alignment)


複数配列と長さLが与えられた時、スコア最大と
なるように各配列から長さLの部分列を抽出
モチーフ（共通の性質を持つ遺伝子などに共通
の部分パターン）抽出などに有用
Sequence 1
Sequence 2
Sequence 3
A A T C G G T
A A T C C G T
A T T C G G A
相対エントロピースコアのもとでの
局所マルチプルアライメント

相対エントロピースコアの定義





fj(a): （モチーフ領域の）j列目におけるaの出現頻度
p(a): aの出現頻度（背景確率）
L: モチーフ領域の長さ
f j (a)
L
score   j 1a f j (a) log
p(a)
実用的アルゴリズム
 Gibbsサンプリング, EM
NP困難
Gibbs サンプリング
1.
2.
3.
4.
5.
各配列xjからランダム
に部分配列tjを選ぶ
1個の配列xiをランダム
に選ぶ
L f (t 'i[ j ])
j

xiの部分列ti’を j 1 p(t'i[ j])
に比例する確率で選ぶ
ti をti’ でおきかえる
ステップ2-4を十分な回
数だけ繰り返す
( ti[j]: 部分列ti のj列目の文字 )
Step 1
x1
t1
x2
x3
t2
t3
Step 2-3
Probabilistic Search
x2
t2’
Step 4
x1
x2
x3
t1
t2’
t3
講義のまとめ（配列アライメント
I）

動的計画法によるペアワイズアライメント





高速配列検索アルゴリズム


大域アライメント
局所アライメント(Smith-Watermanアルゴリズム)
アフィンギャップコストを用いたアライメント
線形領域アライメント
BLAT と PatternHunter
局所マルチプルアライメント(Gibbsサンプリング)

Download Report