数理科学特別講義 バイオインフォマティクスにおける 確率モデル

Presentation on theme: "数理科学特別講義 バイオインフォマティクスにおける 確率モデル"— Presentation transcript:

1 数理科学特別講義 バイオインフォマティクスにおける 確率モデル
数理科学特別講義 バイオインフォマティクスにおける 確率モデル　 阿久津　達也 京都大学　化学研究所 バイオインフォマティクスセンター

2 講義内容 I （基礎編） 分子生物学の基本事項 バイオインフォマティクスに必要な確率統計の基礎 配列アラインメントとスコア関数
隠れマルコフモデル EMアルゴリズム

3 講義内容 II（応用編） 確率文脈自由文法とRNA二次構造予測 進化系統樹構成法 タンパク質立体構造予測 遺伝子ネットワーク推定 実習

4 実習内容 配列検索プログラム マルチプルアラインメント PSI-BLASTによる検索 立体構造予測

5 ゲノム研究と情報科学 ゲノム解析計画の急速な進展 情報解析の必要性 DNA配列：A,C,G,Tの４文字からなる文字列
既に数十種の微生物の全塩基配列が決定 線虫（１０００細胞）、しょうじょうばえも決定 ヒトも概要配列が公開済み 情報解析の必要性 DNA配列⇔プログラムのオブジェクトコード 意味の解析が必要 DNA配列：A,C,G,Tの４文字からなる文字列 様々な情報技術が利用可能

6 遺伝子と蛋白質 遺伝情報の流れ 遺伝子 ゲノム タンパク質 DNA⇒RNA⇒タンパク DNA配列中で直接的に 機能する部分
染色体全体（半数体） 遺伝情報の総体 タンパク質 アミノ酸（２０種類）の鎖

7 DNAとアミノ酸 DNAはA,C,G,Tの４文字の並び DNAは二重ラセン構造⇒相補鎖 塩基：DNA１文字、 残基：アミノ酸１文字
　（アミノ酸は２０種類）

8 アミノ酸と蛋白質 アミノ酸：２０種類 蛋白質：アミノ酸の鎖（短いものはペプチドと呼ばれる）

9 側鎖の例

10 アミノ酸コード表 Ala A アラニン Leu L ロイシン Arg R アルギニン Lys K リシン Asn N アスパラギン Met
メチオニン Asp D アスパラギン酸 Phe F フェニルアラニン Cys C システイン Pro P プロリン Gln Q グルタミン Ser S セリン Glu E グルタミン酸 Thr T トレオニン Gly G グリシン Trp W トリプトファン His H ヒスチジン Tyr Y チロシン Ile I イソロイシン Val V バリン

11 バイオインフォマティクスにおける確率統計
重要なのはデータからのモデル（もしくはパラメータ）の推定 最尤法 ベイズ推定 最大事後確率推定（MAP）

12 最尤推定 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 で尤度は最大

13 ベイズ推定とMAP推定 ベイズ推定：尤度とモデル（パラメータ）の事前確率から、ベイズの定理により、事後確率を推定 最大事後確率（MAP）推定
P(D|θ)P(θ)を最大化するθを計算 P(θ)が一様分布なら最尤推定と同じ

14 不正サイコロのベイズ推定 公正サイコロと不正サイコロ ６が３回続けて出た場合の事後確率 公正：P(i|公正)=1/6
不正：P(6|不正)=1/2,P(i|不正)=1/10 for i≠6 P(公正)=0.99, P(不正)=0.01 ６が３回続けて出た場合の事後確率

15 アラインメント バイオインフォマティクスの最重要技術 ２個もしくは３個以上の配列の類似性の判定に利用
文字間の最適な対応関係を求める（最適化問題） ２個の配列長を同じにするように、ギャップ記号（挿入、欠失に対応）を挿入

16 スコア行列（置換行列） 残基間（アミノ酸文字間）の類似性を表す行列 PAM250, BLOSUM45 など

17 ギャップ無しアラインメントのスコア スコア：P(x,y|M)とP(x,y|R)の対数オッズ比
P(x,y|R): ランダムモデルRにおいて、配列x,yが独立に生起する確率 qa:文字aのRにおける出現確率 P(x,y|M):　一致モデルMにおいて、配列x,yが生起する確率 pab :文字ペアa,bのMにおける出現確率

18 ギャップペナルティ 線形スコア アファインギャップスコア ギャップの確率的解釈 γ(g)=-gd (gはギャップの長さ、 ｄは定数）
　　 ｄは定数） アファインギャップスコア γ(g)=-d-e(g-1)　（d:ギャップ開始ペナルティ, e:ギャップ伸長ペナルティ） ギャップの確率的解釈 P(gap)=f(g)Πqxi γ(g)=log(f(g)) ギャップ長の確率の対数

19 ペアワイズ・アラインメント 配列が２個の場合でも可能なアラインメントの個数は指数オーダー
しかし、スコア最大となるアラインメント（最適アラインメント）は動的計画法により、O(mn)時間で計算可能（m,n:入力配列の長さ）

20 動的計画法による最適アラインメント 最長経路計算による最適アライメント DP (動的計画法)による 最長経路の計算 行列からの経路の復元は、
F(m,n)からmaxで＝となっている F(i,j)を逆にたどることに行う （トレースバック）

21 局所アラインメント 配列の一部のみに共通部分があることが多い 　　⇒共通部分のみのアラインメント （スコア：ギャップ-1、 　　　置換-1、一致1）

22 配列検索の実用プログラム（１） O(mn):mは数百だが、nは数ＧＢにもなる ⇒実用的アルゴリズムの開発
　⇒実用的アルゴリズムの開発 FASTA:短い配列（アミノ酸の場合、1,2文字、DNAの場合、4-6文字）の完全一致をもとに対角線を検索し、さらにそれを両側に伸長し、最後にＤＰを利用 BLAST:固定長（アミノ酸では3, DNAでは11）の全ての類似単語のリストを生成し、ある閾値以上の単語ペアを探し、それをもとに両側に伸長させる。ギャップは入らない。

23 配列検索の実用プログラム（２） SSEARCH: 局所アラインメント（Smith-Watermanアルゴリズム）をそのまま実行
PSI-BLAST: ギャップを扱えるように拡張したBLASTを繰り返し実行。「BLASTで見つかった配列からプロファイルを作り、それをもとに検索」という作業を繰り返す。

24 スコアのベイズ論的解釈 P(M|x,y)をベイズの定理に基づき導出 P(M):２個の配列に関連性がある事前確率
P(R)=1-P(M):ランダムモデルの事前確率 以下の式より尤度比にオッズ比を掛け合わせたものを０と比較すべきであることがわかる

25 スコア行列の導出 基本的には頻度の比の対数をスコアとする BLOSUM行列
既存のスコア行列を用いて多くの配列のアラインメントを求め、ギャップ無しの領域（ブロック）を集める 残基がL％以上一致しているものを同一クラスタに集める 同じクラスタ内で残基aが残基bにアラインされる頻度Aabを計算 qa=∑b Aab / ∑cd Acd, pab=Aab / ∑cd Acd　を求め、　　　ｓ（a,b)=log(pab/qaqb) としたのち、スケーリングし近傍の整数値に丸める

26 マルチプルアラインメント S(mi) = -∑cia log pia （cia= i列におけるaの出現回数,
３本以上の配列が与えられた時、長さが同じで、かつ、スコアが最適となるように各配列にギャップを挿入したもの スコアづけ　（全体スコアは基本的に各列のスコアの和:∑S(mi)） 最小エントロピースコア S(mi) = -∑cia log pia　　　（cia= i列におけるaの出現回数, 　　　　 pia = i列におけるaの生起確率） SPスコア(Sum-of-Pairs) S(mi)=∑k＜ｌ s(mik,mil) （mik = i列, k行目の文字）

27 多次元DPによる マルチプルアラインメント
N個の配列に対するマルチプルアラインメント N次元DPによりO(2NnN)時間（各配列の長さはO(n)を仮定） 例：N=3

28 実用的マルチプルアラインメント法 ヒューリスティックアルゴリズムの開発 プログレッシブアラインメント 逐次改善法
N次元DPは（N=4ですら） 　　　非実用的 一般にはNP困難 プログレッシブアラインメント 近隣結合法などを用いて　案内木を作る 類似度が高い節点から低い節点へという順番で、配列対配列、配列対プロファイル、プロファイル対プロファイルのアラインメントを順次計算 逐次改善法 「配列を一本取り除いては、アラインメントしなおす」を繰り返す

29 マルチプルアラインメントの 近似アルゴリズム(1)
スコアに三角不等式を仮定し、最小化問題として定義 アルゴリズム ∑k≠iS(xk,xi)が最小となるxkを求める S(xk,xi)をもとにギャップを適切に挿入し、マルチプルアラインメントA を構成

30 マルチプルアラインメントの 近似アルゴリズム(2)
定理 A のSPスコア　≦ (2-2/N)・最適解の SPスコア 証明 N・Starのスコア ≦ 2・最適解のスコア (図１でNスターにより、各辺を２回ずつカバー) A のスコア ≦ 　　 (N-１)・Starのスコア 　 (図２でｘ１に接続しない辺のスコアの合計はスター(N-2)個分以下)

31 局所マルチプルアラインメント (Local Multiple Alignment)
複数配列と長さLが与えられた時、スコア最大となるように各配列から長さLの部分列を抽出 モチーフ（共通の性質を持つ遺伝子などに共通の部分パターン）抽出などに有用

32 相対エントロピースコアのもとでの 局所マルチプルアラインメント
相対エントロピースコアの定義 fj(a): （モチーフ領域の）j列目におけるaの出現頻度 p(a): aの出現頻度（事前確率） L: モチーフ領域の長さ 実用的アルゴリズム Gibbsサンプリング, EM NP困難

33 Gibbs サンプリング 各配列xjからランダムに部分配列tjを選ぶ 1個の配列xiをランダムに選ぶ
xiの部分列ti’を　　　　 に比例する確率で選ぶ ti をti’ でおきかえる ステップ2-4を十分な回数だけ繰り返す ( ti[j]: 部分列ti のj列目の文字 )

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

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

36 時々いかさまをするカジノ サイコロの出目だけが観測可能、どちらのサイコロを振っているかは観測不可能
サイコロの出目から、どちらのサイコロを振っているかを推定 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 →途中で公正サイ 　　　　コロに交換

37 Viterbi アルゴリズム(1) 観測列（出力配列データ） x=x1…xLと状態列π=π1…πLが与えられた時、その同時確率は
　P(x,π)=a0 π1Πeπi (xi)aπiπi 但し、πL+1=0 xが与えられた時の、最も尤もらしい状態列は π*=argmaxπ　P(x,π) 例：どちらのサイコロがいつ使われたかを推定

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

39 Viterbiアルゴリズム(3)

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

41 EMアルゴリズムの導出

42 EMアルゴリズムの一般形 初期パラメータ Θ0 を決定。t=0とする。
Q(θ|θt)=∑P(y|x, θt) log P(x,y|θ)　を計算。 Q(θ|θt)を最大化するθ*を計算し、 θt+1 = θ* とする。t=t+1とする。 Qが増大しなくなるまで、２，３を繰り返す。

43 前向きアルゴリズム 配列xの生成確率 P(x)=∑P(x,π) を計算 Viterbiアルゴリズムと類似
fk(i)=P(x1…xi,πi=k) 　をDPにより計算

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

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

46 Baum-WelchのEMによる解釈

47 プロファイルHMM(1) 配列をアラインメントするためのHMM 一致状態(M)、欠失状態(D)、挿入状態(I)を持つ

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