京都大学 化学研究所 バイオインフォマティクスセンター 生命科学基礎論 (第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推定 ベイズ推定：尤度とモデル（パラメータ）の事前確率から、ベイズの定理により、事後確率を推定 最大事後確率（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 ６が３回続けて出た場合の事後確率

隠れマルコフモデル(HMM) ek(b) HMM≒有限オートマトン＋確率 定義 出力記号集合Σ 状態集合 S={1,2,…n} akl 出力確率 ek(b) （開始状態= 終了状態= 0)

HMMにおける基本アルゴリズム Viterbiアルゴリズム Baum-Welchアルゴリズム (EMアルゴリズム） 出力記号列から状態列を推定 Parsing（構文解析） Baum-Welchアルゴリズム 　　(EMアルゴリズム） 出力記号列からパラメータを推定 Learning（学習）

時々いかさまをするカジノ サイコロの出目だけが観測可能、どちらのサイコロを振っているかは観測不可能 サイコロの出目から、どちらのサイコロを振っているかを推定 6,2,6,6,3,6,6,6, 4,6,5,3,6,6,1,2 →不正サイコロ 6,1,5,3,2,4,6,3, 2,2,5,4,1,6,3,4 →公正サイコロ 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,π) 例：どちらのサイコロがいつ使われたかを推定

Viterbiアルゴリズム(2) xから、π*=argmaxπ P(x,π) を計算 そのためにはx1…xiを出力し状態kに至る確率最大の状態列の確率 vk(i) を計算 vk(i)は以下の式に基づき動的計画法で計算

Viterbiアルゴリズム(3)

EM(Expectation Maximization) アルゴリズム 「欠けているデータ」のある場合の最尤推定のための一般的アルゴリズム 最大化は困難であるので、反復により尤度を単調増加させる（θtよりθt+1を計算） HMMの場合、「欠けているデータ」は状態列

EMアルゴリズムの導出

EMアルゴリズムの一般形 初期パラメータ Θ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により計算

後向きアルゴリズム bk(i)= P(xi+1…xL|πi=k) 　をDPにより計算 P(πi=k|x) = fk(i)bk(i)/P(x)

HMMに対するEMアルゴリズム （Baum-Welchアルゴリズム）

Baum-WelchのEMによる解釈

配列アライメント ２個もしくは３個以上の配列の類似性の判定に利用 文字間の最適な対応関係を求める（最適化問題） ２個の場合：ペアワイズアライメント ３個以上の場合：マルチプルアライメント 文字間の最適な対応関係を求める（最適化問題） 配列長を同じにするように、ギャップ記号（挿入、欠失に対応）を挿入 入力配列が定数個（実用上は３個まで）の場合は動的計画法で多項式時間で最適解を計算可能、それ以外の場合はNP困難

プロファイルHMM(1) 配列をアライメントするためのHMM タンパク質配列分類やドメイン予測などに有用 PFAM(http://pfam.wustl.edu/) 一致状態(M)、欠失状態(D)、挿入状態(I)を持つ

プロファイルHMM(2) マルチプル アラインメント プロファイル HMM

参考文献 参考文献 レポート課題のための参考WWWページ 阿久津、浅井、矢田　訳：　バイオインフォマティクス －確率モデルによる遺伝子配列解析―、医学出版 (2000) レポート課題のための参考WWWページ アミノ酸配列データ取得 ゲノムネット（http://www.genome.ad.jp/dbget/dbget.links.html） アミノ酸配列データ: SwissProt タンパク質立体構造データ： PDB 構造予測 CAFASP３参照（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）

レポート課題 インターネット上で利用可能な立体構造予測ソフト（２次構造予測でも可）を２種類以上利用し、得られた結果について比較、考察せよ。ただし、各サーバーに負荷をかけすぎないようにテストデータ（アミノ酸配列）は３種類以下とすること。 提出先：１０号館事務室のレポート提出ＢＯＸ 提出期限：６月２０日(金)