マルコフ過程についての概論


塩基置換に関する数理モデルは,基本的に確率論の中の「確率過程(stochastic process)」,とりわけ「マルコフ過程(Markov process)」の理論に全面的に依拠しています.今回の講義では,確率過程に関する内容に深入りすることはできないのですが,分子進化モデルに直接関係する部分だけでも触れておく方がいいでしょう.

分子系統学での確率過程に関する数学的論議のほとんどは,下記の本:
Cox, D.R. and H.D. Miller (1977)
"The Theory of Stochastic Processes".
Chapman and Hall, London, x+398 pp.
を引用しているので,本書をふまえてざっと説明しましょう.

時間tにおいて確率変数X(t)が離散的な状態i(i=1,2,...,m)のいずれかをとるケースを考える.塩基配列の場合m=4となり,X(t)はある時刻tで特定のサイトがいずれかの状態(塩基)をとるとみなされる.

ここで,変数X(t)の状態遷移に関する確率p[ij](t)を導入する:
p[ij](t)=Prob[X(t)=j|X(0)=i]
※時刻0のときに状態X(0)=iである変数が,
時間tの経過後に状態X(t)=jとなる確率
微小時間Δtにおける状態変化を線形モデルによって近似する(仮定):
 p[ii](Δt)=1+q[ii]・Δt+ο(Δt)
p[ij](Δt)=q[ij]・Δt+ο(Δt) (ただしi≠j)
この線形近似における定数係数q[ij] は瞬間的な状態遷移速度を表す.ο(Δt)はΔtの高次項である.

ここで時刻t+Δtにおける確率p[ik](t+Δt)を考える.時刻0,t,t+Δtの3時点で排反的に可能な状態1)と2)を図示すると:
   0   t   t+Δt
----------------------------------------
1) i → k → k
2) i → j(≠k)→ k
1)の経路での確率は:
p[ik](t)・p[kk](Δt)=p[ik](t)・{1+q[kk]・Δt+ο(Δt)}
となる.

2)の経路での確率は,すべてのj(≠k)に関する総和を取ると:
Σ[j≠k]p[ij](t)・p[jk](Δt)=Σ[j≠k]p[ij](t)・{q[jk]・Δt+ο(Δt)}
となる.

したがって:

p[ik](t+Δt)
=p[ik](t)・{1+q[kk]・Δt+ο(Δt)}
+Σ[j≠k]p[ij](t)・{q[jk]・Δt+ο(Δt)}
=p[ik](t)・{1+q[kk]・Δt}
+Σ[j≠k]p[ij](t)・{q[jk]・Δt}+ο(Δt)

という式を得る.Δtにおける平均変化率:

{p[ik](t+Δt)−p[ik](t)}/Δt
=p[ik](t)・q[kk]+Σ[j≠k]p[ij](t)・q[jk]+ο(Δt)
=Σ[j]p[ij](t)・q[jk]+ο(Δt)

のΔt→0の極限を取ることで(このときο(Δt)→0):

d{p[ik](t)}/dt=Σ[j]p[ij](t)・q[jk]

という1階の線形微分方程式が導かれる.ここでm×m型の正方行列PとQ:

P(t)=(p[11](t) ... p[1m](t))
    (p[21](t) ... p[2m](t))
    (      ...      )
    (p[m1](t) ... p[mm](t))

Q=(q[11] ... q[1m])
  (q[21] ... q[2m])
  (    ...    )
  (q[m1] ... q[mm])

を定義すると,上の微分方程式はもっと簡潔に下記のように表せる:

dP(t)/dt=P(t)・Q
∴一般解 P(t)=C・exp(Qt) C:定数行列

ここで: p[ii](0)=1
p[ij](0)=0 (ただしi≠j)
であるから,

初期値P(0)=I(単位行列)
∴一般解 P(t)=exp(Qt)

右辺はテイラー展開による行列の無限級数として計算される:

P(t)=exp(Qt)=Σ[r=0〜∞]Q^r・t^r/r!

行列Qの冪乗(Q^r)は,Qを固有値と固有ベクトルを用いて対角化することにより求まる.

Q=A(λ1        )A^[-1]
   (  λ2      )
   (     ...   )
   (       λm )

λ1,λ2,...,λm :Qの固有値を対角要素として並べた行列.
         非対角要素はすべて0である.
A:Qの固有ベクトルを並べた行列.A^[-1]はその逆行列

このとき両辺をr乗すると:

Q^r=A(λ1^r         )A^[-1]
    (   λ2^r      )
    (       ...   )
    (        λm^r )

となるから,一般解の式に代入すると,下記のように変形される:

P(t)=Σ[r=0〜∞]Q^r・t^r/r!
   =Σ[r=0〜∞]A(λ1^r         )A^[-1]・t^r/r!
           (   λ2^r      )
           (       ...   )
           (        λm^r )
   =A(Σ[r=0〜∞]λ1^r ・t^r/r!        )A^[-1]
     (   Σ[r=0〜∞]λ2^r ・t^r/r!     )
     (        ...              )
     (        Σ[r=0〜∞]λm^r ・t^r/r!)
   =A(exp(λ1・t)        )A^[-1]
     (   exp(λ2・t)     )
     (       ...      )
     (        exp(λm・t))

最後の式を展開すると,行列P(t)の各要素は,指数関数 exp(λi・t)(i=1,2,...,m)のある一次結合として表現できることがわかる.

----------------------------
時間逆転性をもつマルコフ過程
----------------------------

離散的状態に関するマルコフ連鎖:

dP(t)/dt=P(t)・Q

の一般解:

P(t)=A(exp(λ1・t)        )A^[-1]
     (   exp(λ2・t)     )
     (       ...      )
     (        exp(λm・t))

が得られたとき,t→∞のときのP(t)の各要素の極限値を考える.この極限値が有限確定である場合,すなわち:

lim[t→∞]p[ij](t)=π[ij]

となるとき,極限値π[ij]をこのマルコフ過程の「定常平衡値(stationary equilibrium)」と呼ぶ.π[ij]を並べた行列をΠと表すならば:

lim[t→∞]P(t)=Π

と表記できる.ここで,極限では,最初の状態iがどの状態であっても,最後の状態jをとる確率はすべて互いに等しくなる.すなわち:

lim[t→∞]p[ij](t)=lim[t→∞]Prob[X(t)=j]=πj

であるから,定常平衡値は初期状態iとは無関係になる:

∀i,π[ij]=πj

よって:

Π=(π[11] .... π[1m])=(π1 .... πm)
  (π[21] .... π[2m]) (π1 .... πm)
  (    ....    ) (  ....  )
  (π[m1] .... π[mm]) (π1 .... πm)

となる.したがって,このマルコフ過程における安定平衡値はπi(i=1,2,...,m)によって与えられる.塩基配列の場合,m=4なので,安定平衡値はt→∞においてA,C,G,Tそれぞれの存在頻度πA,πC,πG,πT を意味する.現在,分子系統学で用いられている分子進化モデルのほとんどすべては,マルコフ過程が安定平衡に達していることを前提にしている.

確率過程が進行する時間軸の上での2点t1,t2(t1<t2)を考える.いま,次の条件付き確率を定義する:

Pji=Prob[X(t2)=i|X(t1)=j]

Pjiは「t1で状態jをとるという条件のもとで,将来のt2で状態iをとる確率」である.ここで時間軸を逆にさかのぼる条件付き確率を定義する:

P'ij=Prob[X(t1)=j|X(t2)=i]

このP'ijは「t2で状態iをとるという条件のもとで,過去のt1で状態jをとる確率」である.

条件付き確率の定義により:

P'ij=Prob[X(t1)=j&X(t2)=i]/Prob[X(t2)=i]
=Prob[X(t2)=i|X(t1)=j]・Prob[X(t1)=j]/Prob[X(t2)=i]
=Pji・Prob[X(t1)=j]/Prob[X(t2)=i]

となる.この式は,順方向と逆方向の確率がそれぞれの時点での状態の存在確率を係数として関係づけられることを意味する.

いまt1→∞(したがってt2→∞)とすると:

lim[t1→∞]Prob[X(t1)=j]=πj
lim[t2→∞]Prob[X(t2)=i]=πi

となるので,この極限においては:

P'ij=Pji・πj/πi
∴πi・P'ij=πj・Pji

が成立する.左辺は逆方向(j←i)の遷移確率を,右辺は順方向(j→i)の遷移確率を表わす.両者を等置するこの条件式を「時間逆転性(time-reversibility)」と呼ぶ.時間逆転性とは,安定平衡に達した極限において,順方向と逆方向の遷移確率が互いに等しいという性質である.分子進化モデルにおいて,時間逆転性の仮定は,系統の方向性を前提としない解析をする上で必要になる.

--------------------------------
確率過程としての塩基置換のモデル
--------------------------------

塩基置換を定常マルコフ過程とみなしたとき,時刻tでの塩基A,C,G,Tの極限頻度をそれぞれπ(A),π(C),π(G),π(T)とする.微小時間Δtにおける塩基置換確率の行列Qは下記のように表記される:

t+Δt
A C G T
| ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄
A| q[1,1] q[1,2] q[1,3] q[1,4]
t C| q[2,1] q[2,2] q[2,3] q[2,4]←置換確率行列「Q」
G| q[3,1] q[3,2] q[3,3] q[3,4]
T| q[4,1] q[4,2] q[4,3] q[4,4]

ここで時間逆転性の条件を満たす塩基置換過程を考えると,Δtでの線形近似により:

πi・P'ij=πj・Pji(時間逆転性)
∴πi・q[i,j]Δt=πj・q[j,i]Δt
∴πi・q[i,j]=πj・q[j,i]

となる.この条件式を満足する置換確率行列Qの成分として:

対角成分 t+Δt
A C G T
| ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄
A| μaπ(C) μbπ(G) μcπ(T)
t C|μaπ(A) μdπ(G) μeπ(T)
G|μbπ(A) μdπ(C) μfπ(T)
T|μcπ(A) μeπ(C) μfπ(G)

非対角成分
A→A −μ(aπ(C)+bπ(G)+cπ(T))
C→C −μ(aπ(A)+dπ(G)+eπ(T))
G→G −μ(bπ(A)+dπ(C)+fπ(T))
T→T −μ(cπ(A)+eπ(C)+fπ(G))

で指定される確率過程を「一般時間逆転モデル」(GTR:general time-reversible model)と呼ぶ.ここに,μは全塩基にわたる平均瞬間置換速度であり,各塩基対間の置換速度の相対パラメーターとの積:

t+Δt
A C G T
| ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄
A| μa μb μc
t C| μa μd μe ←置換速度行列「R」
G| μb μd μf
T| μc μe μf

により塩基置換速度を与える.

なお,極限塩基頻度の行列Πを:

| ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄
| π(A) 0 0 0
| 0 π(C) 0 0 ←塩基頻度行列「Π」
| 0 0 π(G) 0
| 0 0 0 π(T)

と書くと,Q=RΠとなる.

GTRには塩基頻度に関するパラメータが4つ(π(A),π(C),π(G),π(T))および置換確率に関するパラメータが6つ(a〜f)含まれている.これらのパラメータ間にさらに強い制約を置くことで,より特定的な塩基置換モデルを指定することができる.たとえば:

1)JC: Jukes & Cantor (1969) 1-parameter model
π(A)=π(C)=π(G)=π(T)=0.25(塩基は等頻度)
a=b=c=d=e=f=1

2)K2P: Kimura (1981) 2-parameter model
π(A)=π(C)=π(G)=π(T)=0.25(塩基は等頻度)
a=c=d=f=1, b=e=κ(transition:transversion比)

3)HKY85: Hasegawa et al. (1985) model
a=c=d=f=1, b=e=κ(transition:transversion比)
※このモデルでは塩基が等頻度という仮定を置かない

のようなモデルがGTRから導出される.

これらのモデルの間には,パラメータに関する制約条件の厳しさを反映する階層的な順序関係がある:GTR→HYY85→K2P→JC.GTRが最も緩いモデルであるのに対し,JCは最も制約の厳しいモデルである.塩基置換モデルには上記以外にもさまざまなモデルがあるが,現在でいずれも GTR モデルを出発点とする階層性の枠内のいずれかの位置を占めている.

どのような分子進化モデルを採用するかは,分子系統解析(最尤法と距離法)にとって大きな問題である.最適な確率モデルをどのように選択するかはかなりやっかいな問題である.一般に,データをうまく説明でき(最尤推定値)しかもパラメータの少ないモデルが望ましいとされる(AIC:赤池情報量基準).塩基置換モデルの場合についても,観察された塩基置換パターンをうまく説明でき,しかも自由パラメータ数ができるだけ少ない単純なモデルを決定する必要がある.これは最節約原理(オッカムの剃刀)の適用にほかならない.

--------------------------------------------------------
Jukes-Cantorモデル(JC)のもとでの塩基置換確率と進化距離
--------------------------------------------------------

微小時間Δtにおける瞬間的な塩基置換速度を記述する「GTRモデル」:

非対角成分 t+Δt
A C G T
| ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄
A| μaπ(C) μbπ(G) μcπ(T)
t C|μaπ(A) μdπ(G) μeπ(T)
G|μbπ(A) μdπ(C) μfπ(T)
T|μcπ(A) μeπ(C) μfπ(G)

対角成分
A→A −μ(aπ(C)+bπ(G)+cπ(T))
C→C −μ(aπ(A)+dπ(G)+eπ(T))
G→G −μ(bπ(A)+dπ(C)+fπ(T))
T→T −μ(cπ(A)+eπ(C)+fπ(G))

π(A),π(C),π(G),π(T)――A,C,G,Tの極限頻度
μ ――平均瞬間置換速度
a〜f ――置換速度の相対パラメータ

において,パラメータ間の更なる制約条件:

1)π(A)=π(C)=π(G)=π(T)=0.25(塩基は等頻度)
2)a=b=c=d=e=f=1

を置くことにより Jukes & Cantor (1969) の1パラメータモデル(JC)を導出することができる:

A C G T
| ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄ ̄
A|−3α α α α
C| α −3α α α
G| α α −3α α
T| α α α −3α ただしα=μ/4

この JC のもとで,時刻t〜t+Δtにおける塩基置換をたどってみる.いま時刻0から時刻tまでの間に塩基iがjに置換する確率を P[ij](t) と書く.たとえば,A→Cの置換確率 P[AC](t) を求めるとしよう.このとき:

t t+Δt
-----------------------------------
A A
C C
G G
T T

より,時刻t〜t+Δtの間で,次のような関係式が成り立つ:

P[AC](t+Δt) =P[AC](t)(1−αΔt)+(1−P[AC](t) )・(α/3)Δt

左辺は「時刻t+ΔtにおいてAからCへの置換が生じている確率」を意味する.右辺第1項は「時刻tですでにCへの置換が生じていて(P[AC](t)),続くΔtでは置換が生じない確率(1−αΔt)」を,右辺第2項は「時刻tではまでCへの置換が生じておらず(1−P[AC](t) ),続くΔtにAまたはGまたはTからCへの置換が生じる確率(3×(1/3)×(α/3)Δt)」をそれぞれ表す.右辺第2項では,AまたはGまたはTである確率(それぞれ1/3)ごとに,置換先がCである確率((α/3)Δt)を乗じた積を排反事象として合計したものである.

上の式を変形すると:

{P[AC](t+Δt) −P[AC](t)}/Δt=α/3−(4α/3)P[AC](t)

Δt→0の極限をとると,線形微分方程式:

dP[AC](t)}/dt=α/3−(4α/3)P[AC](t)

となる.変数分離すると:

d{P[AC](t)−1/4}/dt=−4α/3{P[AC](t)−1/4}

となるので,

log|P[AC](t)−1/4|=−4α/3・t+c (c:定数)
∴ P[AC](t)=1/4+c'exp[−4α/3・t] (c':定数)

任意の塩基iとjに対して,初期値 P[ii](0)=1;P[ij](0)=0を代入して定数c'を決定し,α=μ/4を代入すると,下記のような解が求まる::

P[ii](t)=1/4+3/4・exp[−μ/3t]
P[ij](t)=1/4−1/4・exp[−μ/3t] (i≠j)

2本の塩基配列で対応するあるサイトの塩基が時間t後に異なる塩基をもつ確率D(t)は:

D(t)=(1/4)Σ[i]Σ[j≠i]P[ij](t)
=3/4−3/4・exp[−μ/3t]

ことばで言えば,「そのサイトがもともと塩基iであったとして,時間tの経過によって別の塩基jに置換される確率をjに関して集計して,iに関する期待値(各塩基の存在確率1/4)として表わした値」がD(t)である.

JCのもとでの進化距離を「αt」(=∫αdt)と定義すると,α=μ/4より上の式から:

αt=−3/4log{1−4/3・D(t)}

となる.これが JC のもとでの進化距離を表す式である.

--------
参考文献
--------
マルコフ過程(一般には確率過程)そのものについての参考文献は日本語のものも含めてたくさんあります。しかし、その大部分は確率論・統計学そのものの「数学本」、あるいはシステム制御に関わる「工学本」であって、生物系の読者にはつらいものがあります。

過去を振り返って、確率過程と生物学が密接に関わってきたのは、数理生態学(e.g.個体群動態モデル)と集団遺伝学(e.g.遺伝子頻度変化)です。いずれも日本語の教科書が多数あるでしょう。

分子系統学に直接関わる確率過程論を論じているのはそう多くありません。前に紹介した

Swofford, D.L. et al. (1996), Phylogenetic inference. Pp.407-514 in: D. M.Hillis et al. (eds.), Molecular Systematics, Second Edition. Sinauer, Sunderland.

長谷川政美・岸野洋久(1996)『分子系統学』岩波書店, 東京, x+257pp.

の二つは,現時点でも最良のテキストです。しかし、確率過程そのものに関しては説明を「はしょっている」ように思われます。遠からず刊行される予定の PAUP* のマニュアルには最尤系統推定に関連する確率過程の解説が載るものと期待されます.

系統発生に関係する確率過程について多少とも厳密に追求しているのは、イェール大学の下記博士論文です:

Miklos Cusros (2001) Reconstructing Phylogenies in Markov Models of Sequence Evolution.
http://www.iro.umontreal.ca/~csuros/papers/thesis.pdf

pdfでダウンロードできるこの学位論文の前半部分に、系統マルコフ過程の理論が詳述されています。

それ以外には、オンライン公開されている講義資料も参考になります。たとえば、John Huelsenbeck & Mitsunori Ogiwara のバイオインフォマティクスの講義資料の中には、生物系統学と確率過程論にまたがる内容を含んでいます。時間逆転性を論じた下記:

Time-reversibility
http://www.cs.rochester.edu/users/faculty/ogihara/cs120/slides/main03.pdf

などをご覧ください。

最近の本では:

【書名】Phylogenetics
【著者】Charles Semple and Mike Steel
【刊行】6 February 2003
【出版】Oxford University Press, Oxford
【叢書】Oxford Lecture Series in Mathematics and Its Applications 24
【頁数】xiv+239 pp.
【定価】£45.00 (hardcover)
【ISBN】0-19-850942-1
【備考】目次はここ↓
http://cse.niaes.affrc.go.jp/minaka/files/phylogenetics.html

の第8章が系統発生のマルコフ過程理論に当てられています.