アラインメントはグラフで表現できる

2009年度春学期バイオインフォマティクスアルゴリズム

隠れマルコフモデル (1)

慶應義塾大学

環境情報学部

この世の現象の多くは偶然性に

支配されている？

決定的でない現象の例：

z

コインを投げたときの裏表

z

サイコロの目

z

天気 (降水確率)

z

野球の打者の安打？ (打率)

偶然性にはどのような法則があるのか？

確率モデルとは？

サイコロが出る目の確率モデル：

P(X) = 1/6

確率を取り入れた模型、抽象化

・ 部分的にはでたらめに見える

・ 全体的には秩序があり、それは確率のルールによって決まる

確率分布

z

確率モデルを記述する上で有用な概念

z

確率変数Xが取り得る各々の値に対して、その

起こりやすさを関数Pとして記述する

–

サイコロの目をXとすれば1が出る確率は、P(X=1) =

1/6、より一般的にはP(X) = 1/6

–

コインをトスしたときに出る面をXとすれば、表が出る

確率P(X=head) = ½

–

Etc.

確率分布のパラメータ (1)

z

確率分布の大まかな形(性質？)は決まっているが、一意な形が定

まっていない場合がある

z

形を一意に定める部分を定数とする式が作られる

P(向き=上) = p

P(向き=下) = 1 – p

パラメータ

確率分布のパラメータ (2)

z

確率変数が１つのときは、 P(X=i) = p

_i

で多くの

確率分布を表現できる

–

サイコロの例では、 p

₁

=1/6, p

₂

=1/6, p

₃

=1/6, p

₄

=1/6,

p

₅

=1/6, p

₆

=1/6となる

–

ゆがんだサイコロでは？

z

パラメータが未知の場合はパラメータの推定が

必要

確率モデルの構築

z

モデルの構築

–

P(向き=上) = p

–

P(向き=下) = 1 – p

z

パラメータの推定

–

何度も画鋲を投げて、向きが上、

下になる割合を観測する

–

例えば1000回投げた結果、600

回上向きになったとすれば、推定

値は

とするの

が妥当？

例：画鋲が出る向き

6

.

0

)

(

ˆ

_向き

_{= 上}

₌

P

試行結果の独立性

同じ試行を繰り返すと …

z

_{n-1回の試行で出た}

目をX

₁

, X

₂

, .. X

_n-1

とし

て…

z

_P(X

_n

_{)は過去n-1回の}

試行に関係なく1/6?

6,2,5,1,…

毎日の天気を試行とすると…

天気の移り変わりが穏やかな季節なら…

n日目の天気をX

_n

として、

z

同じ天気は続きやすい？

z

曇りの後は雨になりやすい？

X

₁

X

₂

X

₃

X

₄

X

₅

X

₆

X

₇

X

₈

X

₉

X

₁₀

X

₁₁

X

₁₂

以前の試行結果が次の試行結果に

影響する

z

一般的に雨の確率が20%だとしても、前の日が曇りならその確率は

40%になる？

z

式で表すなら、

–

P(X

_n

= Rainy) = 0.2

–

P(X

_n

= Rainy | X

_n-1

= Cloudy) = 0.4

z

前の試行結果X

₁

,X

₂

,…X

_n-1

に依存して次の試行結果の確率

P(X

_n

|X

₁

,X

₂

,…X

_n-1

) が決まる確率モデルをマルコフモデルと呼ぶ

z

１つ前の試行結果X

_n-1

のみに依存して次の試行結果の確率P(X

_n

|X

_n-1

)

が決まるモデルを一次のマルコフモデルと呼ぶ

一次のマルコフモデルの表現方法

0.4 0.2 0.4 0.4 0.4 0.5 0.1 0.2 0.4

状態遷移図

0.4 0.4 0.2

0.2 0.4 0.4

0.1 0.5 0.4

行列

X

_n-1

X

_n

条件付確率

(S:Sunny, C:Cloudy, R:Rainy)

P(X

_n

=S|X

_n-1

=S)=0.4,P(X

_n

=S|X

_n-1

=C)=0.4, P(X

_n

=S|X

_n-1

=C)=0.2

P(X

_n

=C|X

_n-1

=C)=0.4,P(X

_n

=C|X

_n-1

=S)=0.2, P(X

_n

=C|X

_n-1

=R)=0.4

P(X

_n

=R|X

_n-1

=R)=0.4,P(X

_n

=R|X

_n-1

=C)=0.5, P(X

_n

=R|X

_n-1

=S)=0.1

どんな法則で塩基配列の

出現パターンが決まるのか？

ggagctgcagcccgaccgcggggaggacgccatcgccgcctgcttcctcatcaactgcct

ctacgagcagaacttcgtgtgcaagttcgcgcccagggagggcttcatcaactacctcac

gagggaagtgtaccgctcctaccgccagctgcggacccagggctttggagggtctgggat

ccccaaggcctgggcaggcatagacttgaaggtacaaccccaggaacccctggtgctgaa

ggatgtggaaaacacagattggcgcctactgcggggtgacacggatgtcagggtagagag

gaaagacccaaaccaggtggaactgtggggactcaaggaaggcacctacctgttccagct

gacagtgactagctcagaccacccagaggacacggccaacgtcacagtcactgtgctgtc

caccaagcagacagaagactactgcctcgcatccaacaaggtgggtcgctgccggggctc

tttcccacgctggtactatgaccccacggagcagatctgcaagagtttcgtttatggagg

ctgcttgggcaacaagaacaactaccttcgggaagaagagtgcattctagcctgtcgggg

tgtgcaaggcccctccatggaaaggcgccatccagtgtgctctggcacctgtcagcccac

ccagttccgctgcagcaatggctgctgcatcgacagtttcctggagtgtgacgacacccc

caactgccccgacgcctccgacgaggctgcctgtgaaaaatacacgagtggctttgacga

gctccagcgcatccatttccccagtgacaaagggcactgcgtggacctgccagacacagg

actctgcaaggagagcatcccgcgctggtactacaaccccttcagcgaacactgcgcccg

ctttacctatggtggttgttatggcaacaagaacaactttgaggaagagcagcagtgcct

cgagtcttgtcgcggcatctccaagaaggatgtgtttggcctgaggcgggaaatccccat

ゲノム全体はA,C,G,Tが均等に出現する

ランダム配列か？

大腸菌 Accession Number: NC_000913 Length of Sequence : 4639675 A Content : 1142228 (24.62%) T Content : 1140970 (24.59%) G Content : 1176923 (25.37%) C Content : 1179554 (25.42%) Others : 0 (0.00%) AT Content : 49.21% GC Content : 50.79% マイコプラズマ菌 (M. genitalium) Accession Number: NC_000908 Length of Sequence : 580074 A Content : 200543 (34.57%) T Content : 195695 (33.74%) G Content : 92312 (15.91%) C Content : 91524 (15.78%) Others : 0 (0.00%) AT Content : 68.31% GC Content : 31.69%

特定の部分での塩基の出現パターンの

法則は？

z

ゲノム全体

z

コード領域

z

非コード領域

z

シグナル配列部分

核酸配列、アミノ酸配列のシグナル

z

開始コドン、終止コドン

z

翻訳開始シグナル(SD, Kozak)

z

プロモータ配列

AUG

AGGAGG

AUUCCUCC

16S rRNA

fMet-tRNA

f

Methionine

Shine-Dalgarno

sequence

開始コドン

mRNA

リボソーム

16S rRNAの3‘末端はShine-Dalgarno配列と

対合する

ttacagagtacacaacatcc atg aaacgcatta [thr operon leader peptide:thrL]

aaggtaac

gagg

taacaacc atg cgagtgttga [aspartokinase I]

atggaagtt

aggag

tctgac atg gttaaagttt [homoserine kinase:thrB]

cacgagtactggaaaactaa atg aaactctaca [threonine synthase:thrC]

aatgataa

aaggag

taacct gtg aaaaagatgc [hypothetical protein:b0005]

atttcctgc

aagga

ctggat atg ctgattctta [hypothetical protein:yaaA]

gtttaaagagaaatactatc atg acggacaaat [transaldolase B:talB]

タンパク質のモチーフ(Helix-loop-Helix)

これもアミノ酸の偏りはあるけど曖昧…

1 50

YNP2_CAEEL_6-57 AKR..NARER TRVHTVNQAF .LVLKQHLPS .LR... ... Q9W7E6_51-103 KREMVNAKER LRIRNLNTMF .SRLKRMLPL MQ... ... O55208_60-112 RRRVANAKER ERIKNLNRGF .AKLKALVPF LP... ... ASH3_MOUSE_94-145 IRK.RNERER QRVKCVNEGY .ARLRRHLPE .D... ... O76488_77-138 ARR..NARER NRVKQVNDGF .NALRRHLPA .SVVAALS.. ... AST5_DROME_27-91 IRR..NARER NRVKQVNNGF .SQLRQHIPA .AVIADLSN. ... AST4_DROME_102-163 QRR..NARER NRVKQVNNSF .ARLRQHIPQ .SIITDLT.. ...

51 93

YNP2_CAEEL_6-57 ... ...Q FTKR.VSKLR ILNAAITYID TLL Q9W7E6_51-103 ... ... PDKK.PSKVD TLKAATEYIR LLL O55208_60-112 ... ... QSRK.PSKVD ILKGATEYIQ ILG ASH3_MOUSE_94-145 ... ... YLEKRLSKVE TLRAAIKYIS YLQ O76488_77-138 ... ....GGARRG SGKK.LSKVD TLRMVVEYIR YLQ AST5_DROME_27-91 ... ..GRRGIGPG ANKK.LSKVS TLKMAVEYIR RLQ AST4_DROME_102-163 ... ....KGGGRG PHKK.ISKVD TLRIAVEYIR RLQ

配列パターンのゆらぎ

z

結局ゲノム配列のどこを見ても、A,C,G,Tが均

等にランダムというところはあまりなさそうだ

z

シグナル配列もカッチリ決まっているわけではな

く、ゆらぎがある

z

これら塩基配列の出現の法則を何とかエレガン

トに表せないか？

確率モデ

塩基配列の出現を表す確率モデルは？

z

モデルの構築：各塩基が

それぞれ一定確率で出現

z

P(A) = P(C) = P(G) =

P(T) = ¼

z

種によって塩基の偏りが

ある

–

変異しやすい方向

–

tRNAの性質

–

二次構造およびその回避

–

その他たくさんの理由

z

P(A) = P(T) = 0.34, P(C)

= P(G) = 0.16ではどう

か？

(マイコプラズマ菌)

大腸菌 Accession Number: NC_000913 Length of Sequence : 4639675 A Content : 1142228 (24.62%) T Content : 1140970 (24.59%) G Content : 1176923 (25.37%) C Content : 1179554 (25.42%) Others : 0 (0.00%) AT Content : 49.21% GC Content : 50.79% マイコプラズマ菌 (M. genitalium) Accession Number: NC_000908 Length of Sequence : 580074 A Content : 200543 (34.57%) T Content : 195695 (33.74%) G Content : 92312 (15.91%) C Content : 91524 (15.78%) Others : 0 (0.00%) AT Content : 68.31% GC Content : 31.69%

ある塩基の出現率は前の塩基に影響される

１番目の塩基 ２番目の塩基 ２番目の塩基の出現率 a 0.296 c 0.225 g 0.208 A t 0.271 a 0.276 c 0.230 g 0.294 C t 0.200 a 0.227 c 0.326 g 0.230 G t 0.217 a 0.186 c 0.234 g 0.282 T t 0.298

z

終止コドンがコード領

域中にないから？

z

CpGがTpGに変異？

(高等生物)

z

前の塩基を考慮する

ような確率モデルが

必要

マルコフモデルとは？

z

状態＋状態間の遷移確率

z

前の状態に依存して確率的に次の状態が決ま

る

マルコフモデル

0.6 0.4 1.0 0.7 0.3 1.0

A

G

T

C

A

ACTAが観測される確率は…?

単純なマルコフモデル

A

G

C

T

P(A|G) P(G|A) P(C|T) P(T|C) P(C|A) P(A|C) P(T|G) P(G|T) P(A|A) P(C|C) P(T|T) P(G|G) P(C|G) P(G|C) P(A|T) P(T|A) b1 b2 P(b2|b1) a 0.296 c 0.225 g 0.208 A t 0.271 a 0.276 c 0.230 g 0.294 C t 0.200 a 0.227 c 0.326 g 0.230 G t 0.217 a 0.186 c 0.234 g 0.282 T t 0.298

塩基配列の観測確率 (1)

∏

= − = + −

→

=

→

×

×

→

×

→

×

1 1 1 1 1 3 2 2 1 1

)

(

)

(

)

(

)

(

)

(

)

(

i j j j j i i

k

k

t

k

P

k

k

t

k

k

t

k

k

t

k

P

L

P(k

₁

)を初期状態がk

₁

である確率、t(k→l)を状態kから

状態lへの遷移確率として、経路k

₁

,k

₂

,k

₃

,…,k

_i

を通過す

る確率は

0.6 0.4 1.0 0.3 0.7 1.0

A

G

T

C

A

塩基配列の観測確率 (2)

∏

= = − −

=

×

×

×

×

i j j j j i i

x

x

P

x

P

x

x

P

x

x

P

x

x

P

x

P

2 1 1 1 2 3 1 2 1

)

|

(

)

(

)

|

(

)

|

(

)

|

(

)

(

L

状態とそこから出力される塩基

の種類が1対1対応の場合、長

さiの塩基配列x

₁

,x

₂

,x

₃

,…,x

_i

が観

測される確率は、

A

G

C

T

P(A|G) P(G|A) P(C|T) P(T|C) P(C|A) P(A|C) P(T|G) P(G|T) P(A|A) P(C|C) P(T|T) P(G|G) P(C|G) P(G|C) P(A|T) P(T|A)

どこまで前の状態が影響するか？

z

０次のマルコフモデル

z

一次のマルコフモデル

z

n次のマルコフモデル

隠れマルコフモデルとは？

z

状態とその間をつなぐ遷移確率がある

隠れマルコフモデル

0.6 0.4 1.0 0.3 0.7 1.0 A:0.3 T:0.2 C:0.4 G:0.1 A:0.3 T:0.4 C:0.2 G:0.1 A:0.1 T:0.2 C:0.4 G:0.3 A:0.3 T:0.2 C:0.1 G:0.4 A:0.2 T:0.2 C:0.2 G:0.4

隠れマルコフモデルは以下の２つで定義される

z

t(k→l) … 状態kから状態lへの遷移確率

z

_e

_l

_{(b) … 状態lにおいて、記号bが出力される確率}

隠れマルコフモデルを使って

配列パターンをモデル化したい…

z

構造をどうする？

z

状態遷移確率、記号出力確率をどう決める？

z

与えられた配列の出力確率をどう求める？

z

与えられた配列の通った経路をどう求める？

隠れマルコフモデルの構造の一例

end

start

挿入

欠損

A:0.7 C:0.1 G:0.1 T:0.1 A:0.0 C:0.7 G:0.2 T:0.1 A:0.0 C:0.9 G:0.1 T:0.0 A:0.0 C:0.1 G:0.1 T:0.8

与えられた配列の出力確率は？

経路が分かっている場合

0.6 0.4 1.0 0.3 0.7 1.0 A:0.3 T:0.2 C:0.4 G:0.1 A:0.3 T:0.4 C:0.2 G:0.1 A:0.1 T:0.2 C:0.4 G:0.3 A:0.3 T:0.2 C:0.1 G:0.4 A:0.2 T:0.2 C:0.2 G:0.4

P(“ATCC”)

=0.3 x 0.4 x 0.4 x 0.7

x 0.4 x 1.0 x 0.1

与えられた配列の出力確率は？

経路が分かっていない場合

z

Forwardアルゴリズム

与えられた配列の通った経路は

どうやって求める？

z

通った経路は確率的にしか決まらない

z

最大の確率を与えるような経路はViterbiのアルゴリ

Viterbiのアルゴリズム

0

1

2

4

5

1.0 A:0.1 C:0.3 G:0.5 T:0.1 A:0.7 C:0.1 G:0.1 T:0.1

3

0.8 0.2 0.5 0.5 0.7 0.3 1.0 A:0.4 C:0.1 G:0.3 T:0.2 A:0.1 C:0.1 G:0.1 T:0.7 A:0.2 C:0.2 G:0.3 T:0.3

“ATGA”が観測されたとき、

どの経路を通った可能性が

最も高いか？

v(l, i)の定義

z

_{v(l, i) … 状態lにおいて、塩基配列x}

₁

_,x

₂

_,..x

_i

が観

測されているときの最大確率

0

1

2

4

5

1.0 A:0.1 C:0.3 G:0.5 T:0.1 A:0.7 T:0.1 C:0.1 G:0.1

3

0.8 0.2 0.5 0.5 0.7 0.3 1.0 A:0.4 C:0.1 G:0.3 T:0.2 A:0.1 C:0.1 G:0.1 T:0.7 A:0.2 C:0.2 G:0.3 T:0.3

今、

atgaが観測されたとして、

v(3,2)は状態3の時点でatが

観測されている最大確率

l

k

_max

t(k_max → l) e_l(x_i) e_kmax(x_i-1) 状態lにつながる直前の状態k v(l, i)

p_max= v(k_max, i-1)

k

₁

k

₂

k

₄

v(l,i)を与える直前の状態をk

_max

、そこまでの確

率をp

_max

とすると、

v(l,i)=p

_max

× t(k

_max

→l) ×e

_l

(x

_i

)

各々のkに対して、t(k→l)は一定

Viterbiのアルゴリズム

⎪

⎩

⎪

⎨

⎧

→

−

=

=

≠

=

=

))

(

)

1

,

(

(

max

)

(

0

and

0

if

1

0

and

0

if

0

)

,

(

l

k

t

i

k

v

b

e

l

i

l

i

i

l

v

k

i

l

最大確率と直前のノードの記録

v(l,i)

各状態lとi番目について、

(1)最大スコアv(l,i)と

(2)方角ptr(l,i)を記録しておく

l

k

₁

k

₂

k

₃

“atga”はどの経路を通った？

0

1

2

4

5

1.0 A:0.1 C:0.3 G:0.5 T:0.1 A:0.7 T:0.1 C:0.1 G:0.1

3

0.8 0.2 0.5 0.5 0.7 0.3 1.0 A:0.4 C:0.1 G:0.3 T:0.2 A:0.1 C:0.1 G:0.1 T:0.7 A:0.2 C:0.2 G:0.3 T:0.3

Viterbiのテーブルv(l,i)

状態l

配列のi番目

0

1

2

3

4

0

1

1

0.4

2

0.224

3

0.024

4

0.056

5

0.0392

Viterbiのテーブルptr(l, i)

配列のi番目

状態l

経路：５←４←２←１←０

0

1

2

3

4

0

1

0

2

1

3

1

4

2

5

4

演習問題

0

1

2

4

5

1.0 A:0.1 C:0.3 G:0.5 T:0.1 A:0.7 T:0.1 C:0.1 G:0.1

3

0.8 0.2 0.5 0.5 0.7 0.3 1.0 A:0.4 C:0.1 G:0.3 T:0.2 A:0.1 C:0.1 G:0.1 T:0.7 A:0.2 C:0.2 G:0.3 T:0.3

下記マルコフモデルで、

“agga”が出力されるときの

最大確率およびその経路を

求めよ。v(l,i), ptr(l,i)

を記録すること。

演習問題 解答 v(l, i)

0

1

2

3

4

0

1

0

1

0.4

0

2

0

0.032

3

0.024

0

4

0.0084

5

0.0059

演習問題 解答 ptr(l, i)

0

1

2

3

4

0

1

0

0

2

0

1

3

1

0

4

3

5

4

実装(C言語)

0

1

2

4

5

1.0 A:0.1 C:0.3 G:0.5 T:0.1 A:0.7 C:0.1 G:0.1 T:0.1

3

0.8 0.2 0.5 0.5 0.7 0.3 1.0 A:0.4 C:0.1 G:0.3 T:0.2 A:0.1 C:0.1 G:0.1 T:0.7 A:0.2 C:0.2 G:0.3 T:0.3 static double t[K][K] = { /* 0 1 2 3 4 5 */ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 ], /* 0 */ 0.0, 0.0, 0.8, 0.2, 0.0, 0.0 ], /* 1 */ 0.0, 0.0, 0.5, 0.0, 0.5, 0.0 ], /* 2 */ 0.0, 0.0, 0.0, 0.0, 0.7, 0.3 ], /* 3 */ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 ], /* 4 */ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ] /* 5 */ };

static double e[K][NUC_NUM] = { /* a c g t */ [ 0.00, 0.00, 0.00, 0.00 ], /* 0 */ [ 0.40, 0.10, 0.30, 0.20 ], /* 1 */ [ 0.10, 0.10, 0.10, 0.70 ], /* 2 */ [ 0.20, 0.20, 0.30, 0.30 ], /* 3 */ [ 0.10, 0.30, 0.50, 0.10 ], /* 4 */ [ 0.70, 0.10, 0.10, 0.10 ] /* 5 */ }; static char x[L]; strcpy(x, "atga");

塩基配列

状態遷移確率

記号出力確率

static int type[K]

= ( "TYPE_S", "TYPE_N", "TYPE_N", "TYPE_N", "TYPE_N", "TYPE_N" );

状態のタイプ

double viterbi(int l, int i){ int k, k_max; double p, p_max; switch(type[l]){ case TYPE_S: if(i <= 0){ return 0.0; } else { return LOG0; }

case TYPE_N:

if(i <= 0){ return LOG0; }

for(p_max = LOG0 - 1, k = 0;k < K;k ++){

if(tr[k][l] > 0.0)p = viterbi(k, i - 1) + log(tr[k][l]); else p = LOG0;

if(p > p_max){ p_max = p; k_max = k; } }

return p_max + log(e[l][ (int)cton(x[i - 1]) ]); }

}

実装(Perl)

0

1

2

4

5

1.0 A:0.1 C:0.3 G:0.5 T:0.1 A:0.7 C:0.1 G:0.1 T:0.1

3

0.8 0.2 0.5 0.5 0.7 0.3 1.0 A:0.4 C:0.1 G:0.3 T:0.2 A:0.1 C:0.1 G:0.1 T:0.7 A:0.2 C:0.2 G:0.3 T:0.3 @t = ( # 0 1 2 3 4 5 [ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 ], # 0 [ 0.0, 0.0, 0.8, 0.2, 0.0, 0.0 ], # 1 [ 0.0, 0.0, 0.5, 0.0, 0.5, 0.0 ], # 2 [ 0.0, 0.0, 0.0, 0.0, 0.7, 0.3 ], # 3 [ 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 ], # 4 [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ] # 5 ); @e = ( # a c g t [ 0.00, 0.00, 0.00, 0.00 ], # 0 [ 0.40, 0.10, 0.30, 0.20 ], # 1 [ 0.10, 0.10, 0.10, 0.70 ], # 2 [ 0.20, 0.20, 0.30, 0.30 ], # 3 [ 0.10, 0.30, 0.50, 0.10 ], # 4 [ 0.70, 0.10, 0.10, 0.10 ] # 5 ); @x = split("", "atga");

塩基配列

状態遷移確率

記号出力確率

@type = ( "TYPE_S", "TYPE_N", "TYPE_N", "TYPE_N", "TYPE_N", "TYPE_N" );

sub viterbi($$){ my($l, $i) = @_; my($k, $k_max); my($p, $p_max); if($type[$l] eq "TYPE_S"){ if($i <= 0){ return 0.0; } else { return $::LOG0; } }

elsif($type[$l] eq "TYPE_N"){ if($i <= 0){ return $::LOG0; }

for($p_max = $::LOG0 - 1, $k = 0;$k < $::K;$k ++){

if($tr[$k]->[$l] > 0.0){ $p = viterbi($k, $i - 1) + log($tr[$k]->[$l]); } else { $p = $::LOG0; }

if($p > $p_max){ $p_max = $p; $k_max = $k; } }

return $p_max + log($e[$l]->[ $cton{ $x[$i - 1] } ]); }

}