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

生命科学基礎論
(第9回)
阿久津 達也
京都大学 化学研究所
バイオインフォマティクスセンター
内容



最尤法、ベイズ推定、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(θ )
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
6が3回続けて出た場合の事後確率
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→l)
0.4
0.3
1
akl

出力確率
(開始状態=
終了状態= 0)
0.5
0.5
ek(b)

2
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
出力記号列から状
態列を推定
Parsing(構文解析)
Baum-Welchアルゴ
リズム
(EMアルゴリズム)

0.4
出力記号列からパ
ラメータを推定
Learning(学習)
2
0.3
1
0.5
0.5
0.6
A: 0.2
B: 0.8
0.7
A: 0.1
B: 0.9
A: 0.7
B: 0.3
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,
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
公
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.05
不
0.95
公
0.05
0.1
不
0.9
0.95
0.1
0.9
不
0.9
不
(i  1)  max{e公 (1)  0.95 v公 (i), e公 (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
右辺第1項を Q(θ | θ t )とおくと、
log P( x | θ)  log P( x | θ t ) 
P ( y | x, θ t )
Q(θ | θ )  Q(θ | θ )   P( y | x, θ ) log
P ( y | x, θ )
y
t
t
t
t
最後の項は相対エント ロピーで常に正なので 、
log P( x | θ)  log P( x | θ t )  Q(θ | θ t )  Q(θ t | θ t )
よって、 θ t 1  arg maxQ(θ | θ 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が増大しなくなるまで、2,3を繰り返す。
前向きアルゴリズム
配列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)
後向きアルゴリズム
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 kl e ( x )b (i  1)

k0
k
k
i 1
l
l
k
P( x)  k a0l el ( x1) bl (1)
1
2
3
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
1
i 1
1
12
2
i 1
2
13
3
i 1
3
11
HMMに対するEMアルゴリズム
(Baum-Welchアルゴリズム)
j
x : j番目の配列
Akl: a
E k (b) : 文字bが状態kから現れる回数の期待 値
が使われる回数の期待 値
kl
Akl
Ek

j
(b) 
1
P( x
j

)
i
f (i) akl el ( xi 1) bl (i  1)
1
 P(
j
j
j
k

i|
j
j
b
j
f (i) bk (i )
x ) x 
j
j
k
j
パラメータの更新式
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
配列アライメント

2個もしくは3個以上の配列の類似性の判定に利用





2個の場合:ペアワイズアライメント
3個以上の場合:マルチプルアライメント
文字間の最適な対応関係を求める(最適化問題)
配列長を同じにするように、ギャップ記号(挿入、欠失に対応)を挿入
入力配列が定数個(実用上は3個まで)の場合は動的計画法で多項式時間
で最適解を計算可能、それ以外の場合はNP困難
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
参考文献

参考文献


阿久津、浅井、矢田 訳: バイオインフォマティクス -確率モデルによる遺伝子
配列解析―、医学出版 (2000)
レポート課題のための参考WWWページ

アミノ酸配列データ取得


ゲノムネット(http://www.genome.ad.jp/dbget/dbget.links.html)

アミノ酸配列データ: SwissProt

タンパク質立体構造データ: PDB
構造予測

CAFASP3参照(http://www.cs.bgu.ac.il/~dfischer/CAFASP3/)

GTOP(http://spock.genes.nig.ac.jp/~genome/gtop-j.html)

PHD(http://www.embl-heidelberg.de/predictprotein/predictprotein.html)
レポート課題



インターネット上で利用可能な立体構造予測ソフ
ト(2次構造予測でも可)を2種類以上利用し、得
られた結果について比較、考察せよ。ただし、各
サーバーに負荷をかけすぎないようにテストデー
タ(アミノ酸配列)は3種類以下とすること。
提出先:10号館事務室のレポート提出BOX
提出期限:6月20日(金)