Transcript P - 京都大学

生命情報学基礎論 （４）
隠れマルコフモデル
阿久津 達也
京都大学 化学研究所
バイオインフォマティクスセンター
講義予定













４月１４日（月）: 生命情報学の基盤
４月２１日（月）： 配列の比較と相同性検索
４月２８日（月）： 進化系統樹推定
５月１２日（月）： 隠れマルコフモデル
５月１９日（月）： タンパク質立体構造予測
５月２６日（月）、６月２日（月）： カーネル法
６月９日（月）： 生物情報ネットワークの構造解析
６月１６日（月）: 遺伝子ネットワークの解析と制御（田村）
６月２３日（月）： 代謝ネットワークの堅牢性（田村）
６月３０日（月）： 木の編集距離（田村）
７月７日（月）： タンパク質相互作用予測（林田）
７月１４日（月）： タンパク質複合体予測（林田）
７月１７日（木）： 生物データの圧縮による比較（林田）
内容
配列モチーフ
 最尤推定、ベイズ推定、MAP推定
 隠れマルコフモデル(HMM)
 Viterbiアルゴリズム
 EMアルゴリズム
 Baum-Welchアルゴリズム



前向きアルゴリズム、後向きアルゴリズム
プロファイルHMM
配列モチーフ
モチーフ発見

配列モチーフ : 同じ機能を持つ遺伝子配列などに見
られる共通の文字列パターン
正規表現など文法表現を用いるもの
例： ロイシンジッパーモチーフ
L-x(6)-L-x(6)-L-x(6)-L
ジンクフィンガーモチーフ
C-x(2,4)-C-x(3)-[LIVMFYWC]-x(8)-H-x(3,5)-H
人間にとってわかりやすいが表現力が弱い
確率的な表現法を用いるもの
重み行列（プロファイル）
HMM (隠れマルコフモデル)
人間にとってわかりにくいが
一般に表現力は高い
A
3 .8
-3 .5
1 .2
2 .3
C
1 .5
1 .3
-0 .3
-4 .6
G
-1 .5
-2 .9
4 .2
3 .1
T
0 .2
-4 .1
3 .7
-1 .3
A
C
A
G
G
s co re = 3 .8 + 1 .3 + 4 .2 + 3 .1
C
モチーフの例
• ジンクフィンガーモチーフ
C-x(2,4)-C-x(3)-[LIVMFYWC]-x(8)-H-x(3,5)-H
• ロイシンジッパーモチーフ
L-x(6)-L-x(6)-L-x(6)-L
局所マルチプルアラインメント


複数配列と長さ L が与えられた時、スコア最大と
なるように各配列から長さ L の部分列を抽出
モチーフ発見などに有用
Se qu en c e 1
Se qu en c e 2
Se qu en c e 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: モチーフ領域の長さ
score   Lj 1  a

f j ( a ) log
f j (a )
p (a )
実用的アルゴリズム
 Gibbsサンプリング, EMアルゴリズム
Gibbs サンプリング
1. 各配列 xj からランダム
に部分配列 tj を選ぶ
2. 1個の配列 xi をランダム
に選ぶ
３. xi の部分列 ti’ を
L

j 1
S te p 1
x1
t1
x2
x3
t2
t3
S te p 2 - 3
f j ( t 'i[ j ])
P rob ab ilis tic Se arch
x2
t2 ’
p (t 'i[ j ])
に比例する確率で選ぶ
4. ti をti’ でおきかえる
5. ステップ2-4を十分な回
数だけ繰り返す
( ti[j]: 部分列ti のj列目の文字 )
S te p 4
x1
x2
x3
t1
t2 ’
t3
最尤推定、ベイズ推定、MAP推定
最尤推定

P(D|θ) （尤度）


最尤法


モデルパラメータ θ のもとでのデータ D の出現確率
P(D|θ) を最大化する θ を選ぶ
例




コインを5回投げて、表が3回出た後、裏が2回出た
p(表)=a, p(裏)=1-a とすると、P(D|θ)=a3(1-a)2
a=3/5の時、 P(D|θ) は最大
一般に表が出る頻度を f とすると a=f で尤度は最大
ベイズ推定とMAP推定

ベイズ推定：尤度とモデル（パラメータ）の事前確率
から、ベイズの定理により、事後確率を推定
P (θ | D ) 
P ( D | θ ) P (θ )
P(D)
ただし、

P(D) 

θ'
P ( D | θ' ) P ( θ' ) （ θ が連続値の時）
最大事後確率（MAP）推定


P(D|θ)P(θ) を最大化する θ を計算
P(θ) が一様分布なら最尤推定と同じ
不正サイコロのベイズ推定

公正サイコロと不正サイコロ




公正： P(i|公正)=1/6
不正： P(6|不正)=1/2, P(i|不正)=1/10 for i≠6
P(公正)=0.99, P(不正)=0.01
６が３回続けて出た場合の事後確率
P ( 不正 | 666 ) 
P ( 666 | 不正 ) P ( 不正 )
P ( 666 )
3

( 0 . 5 ) ( 0 . 01 )
( 0 . 5 ) ( 0 . 01 )  ( ) ( 0 . 99 )
3
1
6
3
 0 . 21
隠れマルコフモデル
隠れマルコフモデル(HMM)


HMM≒有限オートマトン＋確率
定義



出力記号集合Σ
状態集合
S={1,2,…n}
遷移確率(k→ｌ)
0.4
0.3
akl

1
出力確率
ek(b)

（開始状態=
終了状態= 0)
2
A : 0 .1
0 .5
0.5
A : 0 .2
B : 0 .8
0.6
0.7
B : 0 .9
A : 0 .7
B : 0 .3
3

HMMにおける基本アルゴリズ
ム
Viterbiアルゴリズ
ム




出力記号列から
パラメータを推定
学習
2
BABBBAB
0.3
出力記号列から
状態列を推定
構文解析
Baum-Welchアル
ゴリズム
(EMアルゴリズム）

0.4
1
A : 0 .1
0 .5
0.7
A : 0 .7
0.5
A : 0 .2
0.6
B : 0 .8
B : 0 .9
B : 0 .3
3
23 12 1 31
2
1
0.4
0.3
3
BABBBAB
ABBABBAAB
BBBABBABAB
BAABBBBA
2
1
A : 0 .1
0.5
0.5
A : 0 .2
B : 0 .8
0.6
0.7
B : 0 .9
A : 0 .7
B : 0 .3
3
時々いかさまをするカジノ





サイコロの出目だけが観測可能、どちらのサイコロを振ってい
るかは観測不可能
サイコロの出目から、どちらのサイコロを振っているかを推定
6,2,6,6,3,6,6,6,
0.95
0.9
4,6,5,3,6,6,1,2
→不正サイコロ
1: 1/6
0.05
1: 1/10
2: 1/6
2: 1/10
6,1,5,3,2,4,6,3,
3: 1/6
3: 1/10
2,2,5,4,1,6,3,4
4: 1/6
4: 1/10
5: 1/6
5: 1/10
0.1
→公正サイコロ
6: 1/6
6: 1/2
6,6,3,6,5,6,6,1,
不 正 サ イコロ
公 正 サ イコロ
5,4,2,3,6,1,5,2
→途中で公正サイコロに交換
Viterbiアルゴリズム
Viterbiアルゴリズム(1)



観測列（出力配列データ） x=x1…xLと状態列π=π1…πLが与えら
れた時、その同時確率は
P(x,π)=a0 π1Πeπi (xi)aπiπi+1 但し、πL+1=0
x が与えられた時、最も尤もらしい状態列は π*=argmaxπ P(x,π)
例：どちらのサイコロがいつ使われたかを推定
0.95
x1= 4
x2= 3
x3= 2
公
0.05
max P ( x 1 x 2
x
不
, π)  P ( x 1 x 2
3
0.95
公
0.05
0.1
0.5
0.9
公
0.05
0
不
π
0.5
0.1
0.95
公
0
0.1
0.9
x
不
0.9
不
1
1
1
,
公公公
)

0
.
5


0
.
95


0
.
95

6
6
6
3
Viterbiアルゴリズム(2)


x から、π*=argmaxπ P(x,π) を計算
そのためには x1…xi を出力し、状態 k に至る確率
最大の状態列の確率 vk(i) を計算
 i

v ( i )  max   e π j ( x j ) a π j 1 π j 
π
 j 1

k

vk(i)は以下の式に基づき動的計画法で計算
v
( i  1) 
l
e (x
l
i 1
)
max ( v
k
k
( i ) a kl )
Viterbiアルゴリズム(3)
0 .95
i -1
i
i+ 1
4
6
1
公
0 .05
公
0.1
( i  1)  max{
公
不
e
公
公
0 .05
0 .1
0 .9
公
0 .95
0.05
不
v
0.95
0 .1
0.9
(1)  0 . 95  v 公 ( i ),
不
e
0 .9
公
不
(1)  0 . 1  v 不 ( i ) }
EMアルゴリズム
EM(Expectation Maximization)アルゴリズム

「欠けているデータ」のある場合の最尤推定のため
の一般的アルゴリズム
x : 観測データ、
y : 欠けているデータ、
θ : パラメータ集合
目標：
log P ( x|θ )  log
 P ( x,y| θ ) の最大化
y


最大化は困難であるので、反復により尤度を単調
増加させる（θtよりθt+1を計算）
HMMの場合、「欠けているデータ」は状態列
EMアルゴリズムの導出
log P ( x | θ )  log P ( x , y | θ )  log P ( y | x , θ )
t
両辺に P ( y | x , θ )をかけて
log P ( x | θ ) 

y についての和をとり、
P ( y | x , θ ) log P ( x , y | θ )   P ( y | x ,θ ) log P ( y | x , θ )
t
t
y
右辺第１項を
y
t
Q ( θ | θ )とおくと、
log P ( x | θ )  log P ( x | θ ) 
t
Q (θ | θ )  Q (θ | θ ) 
t
t
t

t
t
P ( y | x , θ ) log
y
最後の項は相対エント
P ( y | x, θ )
P ( y | x, θ )
ロピーで常に正なので
、
log P ( x | θ )  log P ( x | θ )  Q ( θ | θ )  Q ( θ | θ )
t
よって、
θ
t 1
t
t
t
 arg max Q ( θ | θ ) とすれば尤度は増大
t
θ
EMアルゴリズムの一般形
1.
2.
3.
4.
初期パラメータ θ0 を決定。t=0とする
Q(θ|θt)=∑P(y|x, θt) log P(x,y|θ) を計算
Q(θ|θt)を最大化するθ*を計算し、
θt+1 = θ* とする。t=t+1とする
Qが増大しなくなるまで、２，３を繰り返す
前向きアルゴリズム
配列 x の生成確率
P(x)=∑P(x,π) を計算
 Viterbiアルゴリズム
と類似
 fk(i)=P(x1…xi,πi=k)
をDPにより計算

f 0 ( 0 )  1, f k ( 0 )  0
f l (i ) 
e ( x )
l
i
f k ( i  1) a kl
k
P ( x) 
1
2

k
a 11
a 21
f k (L) a k0
1
f 1 (i ) 
2
a 31
3
3
(f1
f
2
f
3
e (x )
( i  1) a 
( i  1) a 
( i  1) a )
1
i
11
21
31
後向きアルゴリズム
bk(i)=
P(xi+1…xL|πi=k)
をDPにより計算
 P(πi=k|x) =
fk(i)bk(i)/P(x)

b
(L) 
a
b ( i )   a e ( x ) b ( i  1)
P ( x )   a e ( x ) b (1)
k0
k
kl
k
k
0l
k
1
a 11
a 12
2
1
2
a 13
3
l
b
1
k
1
l
(i ) 
( a 11 e 1 ( x i  1) b 1 ( i  1) 
) b 2 ( i  1) 
a e (x
i 1
a e (x
i 1
12
3
i 1
l
13
2
3
) b 3 ( i  1))
Viterbi と前向きアルゴリズムの比較
n
P ( x, π | θ )   a π
i 1

i 1
e
π i π i, xi
a 11
Viterbiアルゴリズム
max
 P ( x , π |θ ) 
π

Forwardアルゴリズム
  P ( x , π |θ )
π

A
e 1,A
1
e 1,B
B
a 12
a 21
a 22
2
e 2,A
e 2,C
A
C
Π fo r “B A C A ”=
{ 1 12 1 , 1 12 2 , 1 22 1 ,12 22 }
HMMに対するEMアルゴリズム
（Baum-Welchアルゴリズム）
値
x
A kl : a が使われる回数の期待
E k ( b ) : 文字 b が状態 k から現れる回数の期待
j
: j 番目の配列
kl
A kl  
j
Ek
(b ) 
1
j
P(x )

j

値
j
j
f k ( i ) a kl e l ( x i  1) b l ( i  1)
j
i
1
j
P(x )

j
j
f k (i ) b k (i )
i | x  b 
j
j
パラメータの更新式
aˆ

kl
A
 A
eˆ k ( b ) 
kl
l'
kl '
E (b )
 E (b ' )
k
b'
k
Baum-WelchのEMによる解釈
  e
E k ( b ,π )
M
P ( x, π | θ ) 
k 1
k
(b )

M
M
a
k 0
b
Q (θ | θ ) 
t

l 1
A kl ( π )
および
kl
t
P ( π | x , θ ) log P ( x , π | θ )
より、
π
M

t
P (π | x, θ )   
 k 1 b

Q (θ | θ ) 
t
π
M
E k ( b , π ) log
M
E

k 1
ここで
 p
i
M
k
( b ) log
e
k
(b ) 
log
i
e
q
は
i
(b ) 
k
q


k  0 l 1
A kl ( π ) log

a kl 

M
 A
k  0 l 1
b
e k (b ) 
M
kl
log
a
kl
p の時、最大より、
i
i
E (b ) ,
a
(
b
'
)
E
k
kl

k
A
 A
kl
l'
kl
b'
π の時、 a が使われる回数
( b ,  ) : パスが π の時、 b が状態 k から現れる回数
A kl ( ): パスが
Ek
kl
配列は１個のみを仮定
プロファイルHMM
配列アラインメント

２個もしくは３個以上の配列の類似性の判定に利用




２個の場合：ペアワイズアラインメント
３個以上の場合：マルチプルアラインメント
文字間の最適な対応関係を求める（最適化問題）
配列長が同じになるよう、ギャップ記号を挿入
HBA_HUM AN
HBB_HUM AN
MYG_PHY CA
GLB5_PE TMA
LGB2_LU PLU
GLB1_GL YDI
VGAHA GEY
VNVDE V
VEADV AGH
VYSTY ETA
FNANI PKH
IAGAD NGAGV
HBA_HUM AN
HBB_HUM AN
MYG_PHY CA
GLB5_PE TMA
LGB2_LU PLU
GLB1_GL YDI
V
V
V
V
F
I
G
E
Y
N
A
A
A
S
A
G
A
D
H
N
D
T
N
N
A
V
V
Y
I
G
G
D
A
E
P
A
E
E
G
T
K
G
Y
V
H
A
H
V
プロファイルHMM (1)


配列をアラインメントするためのHMM
タンパク質配列分類やドメイン予測などに有用

例：ドメインの種類ごとにHMMを作る


PFAM(http://pfam.wustl.edu/)
一致状態(M)、欠失状態(D)、挿入状態(I)を持つ
D
D
D
I
I
I
I
B E G IN
M
M
M
E ND
プロファイルHMM (2)
マルチプル
アラインメント
M M ・
プロファイル
HMM
・
・ M
こ うもり
A G － － － C
ラット
A － A G － C
ネコ
A G － A A －
ハエ
－ － A A A C
ヤギ
A G － － － C
D
D
D
I
I
I
I
B E G IN
M
M
M
E ND
プロファイルHMM (3)


各配列ファミリーごとに HMM を作成
スコア最大のHMMのファミリーに属すると予測
w in !
K no w n Se q. (Tra ining D a ta )
Sco re = 1 9 .5
EM
class 1
HMM1
V ite rbi
EM
class 2
N e w Se q.
(Te s t Da ta )
V ite rbi
HMM2
Sco re = 1 5 .8
まとめ

配列モチーフ



局所マルチプルアラインメント
Gibbsサンプリング
HMMによる配列解析




最尤推定、ベイズ推定、MAP推定
隠れマルコフモデル(HMM)
Viterbiアルゴリズム
Baum-Welchアルゴリズム



EMアルゴリズムに基づく
前向きアルゴリズム、後向きアルゴリズム
プロファイルHMM