生命情報学基礎論 （４）
隠れマルコフモデル
阿久津 達也
京都大学 化学研究所
バイオインフォマティクスセンター
講義予定
４月１３日（月）: 生命情報学の基盤
４月２０日（月）： 配列の比較と相同性検索
４月２７日（月）： 進化系統樹推定
５月１日（金）： 隠れマルコフモデル
５月１１日（月）： 遺伝子ネットワークの解析と制御（田村）
５月１８日（月）： タンパク質立体構造予測
５月２５日（月）、６月１日（月）： カーネル法
６月８日（月）： 代謝ネットワークの堅牢性（田村）
６月１５日（月）： 生物情報ネットワークの構造解析
６月２２日（月）： 木の編集距離（田村）
６月２９日（月）： タンパク質相互作用予測（林田）
７月６日（月）： タンパク質複合体予測（林田）
７月１３日（月）： 生物データの圧縮による比較（林田）
内容
配列モチーフ
最尤推定、ベイズ推定、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
T
-1.5
-2.9
4.2
3.1
0.2
-4.1
3.7
-1.3
A
C
G
G
A
score= 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 の部分列を抽出
モチーフ発見などに有用
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: モチーフ領域の長さ
score Lj1a
f j (a)
f j (a) log
p(a)
実用的アルゴリズム
Gibbsサンプリング, EMアルゴリズム
Gibbs サンプリング
1. 各配列 xj からランダム
に部分配列 tj を選ぶ
2. 1個の配列 xi をランダム
に選ぶ
３. xi の部分列 ti’ を
Step 1
x1
x2
x3
Step 2-3
f j (t 'i[ j])
j 1 p(t 'i[ j ])
( ti[j]: 部分列ti のj列目の文字 )
t2
t3
L
に比例する確率で選ぶ
4. ti をti’ でおきかえる
5. ステップ2-4を十分な回
数だけ繰り返す
t1
Probabilistic Search
x2
t2’
Step 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(θ )
P(θ | D)
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(不正)
P(不正 | 666)
P(666)
(0.5) 3 (0.01)
0.21
3
1 3
(0.5) (0.01) ( 6 ) (0.99)
隠れマルコフモデル
隠れマルコフモデル(HMM)
HMM≒有限オートマトン＋確率
定義
出力記号集合Σ
状態集合
S={1,2,…n}
遷移確率(k→ｌ)
0.4
0.3
akl
1
出力確率
ek(b)
（開始状態=
終了状態= 0)
2
0.5
0.5
A: 0.2
B: 0.8
0.6
3
0.7
A: 0.1
B: 0.9
A: 0.7
B: 0.3
HMMにおける基本アルゴリズ
ム
Viterbiアルゴリズ
ム
BABBBAB
出力記号列から
状態列を推定
構文解析
Baum-Welchアル
ゴリズム
(EMアルゴリズム）
0.4
出力記号列から
パラメータを推定
学習
2
0.3
1
0.5
0.7
A: 0.7
B: 0.3
0.5
A: 0.2
B: 0.8
0.6
A: 0.1
B: 0.9
3
2312131
2
1
0.4
0.3
3
BABBBAB
ABBABBAAB
BBBABBABAB
BAABBBBA
2
1
0.5 0.7
0.5
A: 0.2
B: 0.8
0.6
3
A: 0.1
B: 0.9
A: 0.7
B: 0.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
公
0.05
不
0.9
0.5
0.1
公
x2=3
0.95
公
0.05
0
不
0.95
公
0.05
0.1
0.5
x3=2
0
0.1
0.9
不
0.9
不
max P( x1 x2 x3 , π) P( x1 x2 x3 , 公公公) 0.5 16 0.95 16 0.95 16
π
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) e ( x
l
l
)
i 1
max (v (i) a
k
k
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, θ )
両辺に P( y | x, θ t )をかけて yについての和をとり、
log P( x | θ ) P( y | x, θ t ) log P( x, y | θ ) P( y | x,θ t ) log P( y | x, θ )
y
y
右辺第１項を Q(θ | θ t )とおくと、
log P( x | θ ) log P( x | θ t )
t
P
(
y
|
x
,
θ
)
t
t
t
t
Q(θ | θ ) Q(θ | θ ) P( y | x, θ ) log
P ( y | x, θ )
y
最後の項は相対エント ロピーで常に非負なの で、
log P( x | θ ) log P( x | θ t ) Q(θ | θ t ) Q(θ t | θ t )
よって、 θ t 1 arg max Q(θ | θ t ) とすれば尤度は非減少
θ
p log( p / q ) p log( q / p ) log e p ((q / p ) 1) ( q p ) log e 0
i
i
i
i
i
i
i
i
i
i
i
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) el ( xi ) f k (i 1) akl
k
P( x) k f k ( L) ak 0
1 a11
2
3
a21
a31
1
2
3
f 1 (i ) e1 ( xi )
( f 1 (i 1) a11
f 2 (i 1) a 21
f 3 (i 1) a31)
後向きアルゴリズム
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)
k
k0
k
k
k
1
kl
0l
i 1
l
l
1
k
l
a11 1 b (i)
a12
(a e ( x ) b (i 1)
2
a e ( x ) b (i 1)
a13
3 a e ( x ) b (i 1))
1
2
3
1
i 1
1
12
2
i 1
2
13
3
i 1
3
11
Viterbi と前向きアルゴリズムの比較
P( x, π | θ ) aπ
n
i 1
e
π πi, xi
i 1 i
Viterbiアルゴリズム
max P( x, π|θ )
π
a11
Forwardアルゴリズム
P( x, π|θ )
π
e1,A
1
e1,B
A
B
a12
a21
a22
2
e2,A
e2,C
A
C
Π for “BACA”=
{ 1121, 1122, 1221,1222 }
HMMに対するEMアルゴリズム
（Baum-Welchアルゴリズム）
j
Akl: a が使われる回数の期待 値 x : j番目の配列
E k (b) : 文字bが状態 kから現れる回数の期待 値
kl
Akl
Ek
j
(b)
1
P( x
j
)
i
j
f k (i ) a kl el ( xi 1) bl (i 1)
j
1
P(
j
j
j
j
f k (i ) bk (i )
x ) x
j
i|
j
i
b
パラメータの更新式
ˆa A
A
kl
kl
l'
kl '
(b)
E
eˆ (b) (b' )
E
k
k
b'
k
Baum-WelchのEMによる解釈
M
P ( x, π | θ )
k 1
b
E k ( b ,π )
e (b)
k
M
M
Akl ( π ) および
akl
k 0 l 1
Q(θ | θ t ) P(π | x, θ t ) log P( x, π | θ ) より、
π
M M
M
Q (θ | θ ) P (π | x, θ ) E k (b, π ) log ek (b) Akl (π ) log a kl
π
k 0 l 1
k 1 b
t
t
M
M
M
E k (b) log ek (b) Akl log a kl
k 1 b
ここで
k 0 l 1
p log q は q p の時、最大より、
(b)
E
A
(
b
)
,
e
a
(
b
'
)
E
A
i
i
i
i
i
k
kl
kl
k
k
l'
kl
b'
Akl ( ): パスが πの時、 a kl が使われる回数
Ek (b, ) : パスが πの時、 b が状態 k から現れる回数
配列は１個のみを仮定
プロファイルHMM
配列アラインメント
２個もしくは３個以上の配列の類似性の判定に利用
２個の場合：ペアワイズアラインメント
３個以上の場合：マルチプルアラインメント
文字間の最適な対応関係を求める（最適化問題）
配列長が同じになるよう、ギャップ記号を挿入
HBA_HUMAN
HBB_HUMAN
MYG_PHYCA
GLB5_PETMA
LGB2_LUPLU
GLB1_GLYDI
VGAHAGEY
VNVDEV
VEADVAGH
VYSTYETA
FNANIPKH
IAGADNGAGV
HBA_HUMAN
HBB_HUMAN
MYG_PHYCA
GLB5_PETMA
LGB2_LUPLU
GLB1_GLYDI
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
BEGIN
M
M
M
END
プロファイルHMM (2)
マルチプル
アラインメント
プロファイル
HMM
M M ・ ・ ・ M
D
D
D
I
I
I
I
BEGIN
M
M
M
こうもり A G － － － C
ラット
A －A G －C
ネコ
A G －A A －
ハエ
－－A A A C
ヤギ
A G －－－C
END
プロファイルHMM (3)
各配列ファミリーごとに HMM を作成
スコア最大のHMMのファミリーに属すると予測
win !
Known Seq. (Training Data)
Score= 19.5
EM
class 1
HMM1
Viterbi
EM
class 2
New Seq.
(Test Data)
Viterbi
HMM2
Score= 15.8
まとめ
配列モチーフ
局所マルチプルアラインメント
Gibbsサンプリング
HMMによる配列解析
最尤推定、ベイズ推定、MAP推定
隠れマルコフモデル(HMM)
Viterbiアルゴリズム
Baum-Welchアルゴリズム
EMアルゴリズムに基づく
前向きアルゴリズム、後向きアルゴリズム
プロファイルHMM
© Copyright 2024 ExpyDoc
