奈良女子大集中講義 - Kyoto University Bioinformatics

奈良女子大集中講義
バイオインフォマティクス (6)
モチーフ発見・隠れマルコフモデル
阿久津達也
京都大学化学研究所
バイオインフォマティクスセンター
講義予定
• 9月5日
–
–
–
–
分子生物学概観
分子生物学データベース
配列アラインメント
実習１（データベース検索と配列アラインメント）
• 9月6日
–
–
–
–
モチーフ発見
隠れマルコフモデル
カーネル法
進化系統樹推定
• 9月7日
–
–
–
–
–
RNA二次構造予測
タンパク質立体構造予測
ネットワーク推定
スケールフリーネットワーク
実習２（構造予測）
内容
•
•
•
•
•
•
配列モチーフ
最尤推定、ベイズ推定、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   Lj1a f j (a) log

実用的アルゴリズム
 Gibbsサンプリング, EM
f j (a)
p(a)
相対エントロピースコア
score 
L
 j 1a
f j (a) log
f j (a)
p(a)
ti
A A T C G G T
s1
s2
s3
•
•
•
•
A A T C C G T
A T T C G G A
L=7
f1(A)=1, f1(C)=f1(G)=f1(T)=0
f2(A)=2/3, f2(T)=1/3, f2(C)=f2(G)=0, …
p(A)=p(C )=p(G)=p(T)=1/4
Gibbs サンプリング
1. 各配列 xj からランダム
に部分配列 tj を選ぶ
2. 1個の配列 xi をランダム
に選ぶ
３. xi の部分列 ti’ を
Step 1
x1
x2
x3
Step 2-3
Probabilistic Search
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
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アルゴリズ
ム
0.4
BABBBAB
– 出力記号列から
パラメータを推定
– 学習
0.3
1
0.5
0.7
A: 0.2
B: 0.8
0.6
A: 0.1
B: 0.9
A: 0.7
B: 0.3
0.5
– 出力記号列から
状態列を推定
– 構文解析
• Baum-Welchアル
ゴリズム
(EMアルゴリズム）
2
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 アルゴリズム(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) を計算
v (i)  max
k
π
i 1
a0π1  eπ j( x j)a π j π j 1
j 1
• vk(i)は以下の式に基づき動的計画法で計算
v (i  1)  e ( x
l
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
公
v
公
0.95
公
0.05
不
公
0.05
0.1
不
0.9
0.95
0.1
0.9
(i  1)  max{ e公 (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, θ)
両辺に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 maxQ(θ | θ t ) とすれば尤度は増大
θ
EMアルゴリズムの一般形
1. 初期パラメータ Θ0 を決定。t=0とする
2. Q(θ|θt)=∑P(y|x, θt) log P(x,y|θ) を計算
3. Q(θ|θt)を最大化するθ*を計算し、 θt+1 =
θ* とする。t=t+1とする
4. 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 ) a k 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) a 31)
Viterbi と前向きアルゴリズムの比較
Pr( x, π | θ)   aπ
n
i 1
e
π πi, xi
i 1 i
a11
• Viterbiアルゴリズム
max  Pr( x,π|θ) 
π
• Forwardアルゴリズム
 Pr( x,π|θ) 
π
e1,A
1
e1,B
A
B
a12
a21
a22
2
e2,A
e2,C
A
C
Π for “BACA”=
{ 1121, 1122, 1221,1222 }
後向きアルゴリズム
• 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
2
3
kl
0l
a11 1
a12
a13
2
3
i 1
l
l
1
k
l
b (i) 
(a e ( x ) b (i  1) 
a e ( x ) b (i  1) 
a e ( x ) b (i  1))
1
1
i 1
1
12
2
i 1
2
13
3
i 1
3
11
HMMに対するEMアルゴリズム
（Baum-Welchアルゴリズム）
j
Akl: a が使われる回数の期待値 x : j番目の配列
E k (b) : 文字bが状態kから現れる回数の期待値
kl
Akl  
j
Ek
(b) 
1
P( x
j
f

)
j
k
i
1
 P(
j
(i) akl el ( xi 1) bl (i  1)
j

j
j
j
f k (i) bk (i )
x ) x 
j
i|
j
j
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
b'
l'
kl
配列アラインメント
• ２個もしくは３個以上の配列の類似性の判定に利用
– ２個の場合：ペアワイズアラインメント
– ３個以上の場合：マルチプルアラインメント
• 文字間の最適な対応関係を求める（最適化問題）
• 配列長が同じになるよう、ギャップ記号を挿入
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
まとめ
• HMMによる配列解析
–
–
–
–
–
最尤推定、ベイズ推定、MAP推定
隠れマルコフモデル(HMM)
Viterbiアルゴリズム
EMアルゴリズム
Baum-Welchアルゴリズム
• 前向きアルゴリズム、後向きアルゴリズム
– プロファイルHMM
• 参考文献
– 阿久津、浅井、矢田訳：バイオインフォマティクス－確率モデルによる
遺伝子配列解析―、医学出版、2001
– 阿久津：バイオインフォマティクスの数理とアルゴリズム、共立出版、
2007

Download Report