生命科学基礎論 (第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日(金)
© Copyright 2024 ExpyDoc