Transcript 第3回 - 京都大学

分子生物情報学(3)
確率モデル（隠れマルコフモデル）に
基づく配列解析
阿久津 達也
京都大学 化学研究所
バイオインフォマティクスセンター
内容



最尤法、ベイズ推定、MAP推定
隠れマルコフモデルによる推定
文脈自由文法によるRNA二次構造予測
バイオインフォマティクスにおけ
る確率統計

重要なのはデータからのモデル（もしくは
パラメータ）の推定



最尤法
ベイズ推定
最大事後確率推定（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
出力確率
（開始状態=
終了状態= 0)
A : 0 .1
0 .5
0.5
ek(b)

2
A : 0 .2
B : 0 .8
0.6
0.7
B : 0 .9
A : 0 .7
B : 0 .3
3
HMMにおける基本アルゴリズム

Viterbiアルゴリズム




出力記号列からパ
ラメータを推定
Learning（学習）
2
BABBBAB
0.3
出力記号列から状
態列を推定
Parsing（構文解析）
Baum-Welchアルゴ
リズム
(EMアルゴリズム）

0.4
1
A : 0 .1
0 .5
0.5
0.6
A : 0 .2
B : 0 .8
0.7
B : 0 .9
A : 0 .7
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,
4,6,5,3,6,6,1,2
0.95
0.9
→不正サイコロ
6,1,5,3,2,4,6,3,
1: 1/6
0.05
1: 1/10
2,2,5,4,1,6,3,4
2: 1/6
2: 1/10
3: 1/6
3: 1/10
→公正サイコロ
4: 1/6
4: 1/10
5: 1/6
5: 1/10
6,6,3,6,5,6,6,1,
0.1
6: 1/6
6: 1/2
5,4,2,3,6,1,5,2
→途中で公正サイ
不 正 サ イコロ
公 正 サ イコロ
コロに交換
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) を計算
v

(i ) 
k
max
π
i 1
a 0π1  e π j ( x j ) a π j π j 1
j 1
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(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 | θ ) 
t
y

t
P ( y | x ,θ ) log P ( y | x , θ)
y
t
Q (θ | θ )とおくと、
右辺第１項を
log P ( x | θ )  log P ( x | θ ) 
t
Q (θ | θ )  Q (θ | θ ) 
t
t
t
 P ( y | x, θ
t
t
y
最後の項は相対エント
) log
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
2
a 31
3
3
f 1 (i ) 
(f1
f
2
f
3
e (x )
( i  1) a 
( i  1) a 
( i  1) a )
1
i
11
21
31
Viterbi と前向きアルゴリズムの比較
n
Pr( x , π | θ)   a π
i 1

i 1
e
π π,x
i
i
i
a 11
Viterbiアルゴリズム
max  Pr( x , π |θ) 
1
Forwardアルゴリズム
  Pr( x , π |θ)
π

e 1,B
B
a 12
a 21
π

A
e 1,A
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 }
後向きアルゴリズム
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
k
k
1
a 11
a 12
2
k
0l
1
b
2
a 13
3
kl
l
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))
HMMに対するEMアルゴリズム
（Baum-Welchアルゴリズム）
が使われる回数の期待
A kl : a
E k ( b ) : 文字 b が状態
値
x
kl
A kl 
j
Ek
(b ) 
j
P(x )

j
: j 番目の配列
k から現れる回数の期待
1

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
(b )
k

M
a
k 0
b
 P (π | x, θ
Q (θ | θ ) 
t
M
l 1
t
A kl ( π )
および
kl
) log P ( x , π | θ )
より、
π
M
P (π | x, θ )   
 k 1 b

Q (θ | θ ) 
t
t
π
M
E
M
E

k 1
ここで
 p
i
log
e
q
k
( b ) log
は
i
(b ) 
k
( b , π ) log
M
e
k
(b ) 
q

k
(b ) 
 A
kl
log
a
i
E (b ) ,
a
 E (b ' )
k
k
kl

A
 A
kl
l'

k  0 l 1
p の時、最大より、
i
b'
e
M
k  0 l 1
b
i
k
M
kl
kl

A kl ( π ) log a kl 

配列アライメント

２個もしくは３個以上の配列の類似性の判定に利用





２個の場合：ペアワイズアライメント
３個以上の場合：マルチプルアライメント
文字間の最適な対応関係を求める（最適化問題）
配列長を同じにするように、ギャップ記号（挿入、欠失に対応）を挿入
入力配列が定数個（実用上は３個まで）の場合は動的計画法で多項式時間
で最適解を計算可能、それ以外の場合はNP困難
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
講義のまとめ(HMM)

HMMによる配列解析





最尤推定、ベイズ推定、MAP推定
隠れマルコフモデル(HMM)
Viterbiアルゴリズム
EMアルゴリズム
Baum-Welchアルゴリズム、

前向きアルゴリズム、後向きアルゴリズム
プロファイルHMM
参考文献



阿久津、浅井、矢田 訳： バイオインフォマティクス －確率モデルに
よる遺伝子配列解析―、医学出版、2001