In [1]:
import numpy as np
import matplotlib.pyplot as plt
塩基置換のマルコフモデル¶
| $$P(t)$$ | $$Q$$ |
|---|---|
| 離散 | 連続（に限りなく近づけた極限） |
| 遷移確率行列 | 置換速度行列（遷移強度行列） |
| 時間 $t$ 後に塩基 $i$ が $j$ になっている確率を $P_{ij}(t)=p_{ij}(t)$ とする。 | 間隔 $(t,t+\Delta t)$ に起きた遷移 $i\rightarrow j$ の確率が $q_{ij}\Delta t$ |
| 定義より、$t=0$ の遷移確率行列 $P(0)$ は単位行列 $I$ | $Q$ は積の形 $Qt$ のみで現れ、積のみが意味を持つ。（＝相対的な置換速度しか表せない。 |
チャップマン-コルモゴロフの定理¶
以下が成立することは、定義から明らか。
$$ \begin{aligned} &p_{ij}(s+t) = \sum_{k}p_{ik}(s)p_{kj}(t)\\ &\Longleftrightarrow P(s+t) = P(s)P(t) \end{aligned} $$$P(t)$ と $Q$ の関係¶
$P(t)$ は微分可能であり、上式を用いて導関数を計算することができる。
$$\begin{aligned} P(t+\Delta t) &= P(t)P(\Delta t) \\ &= P(t)P(0)(I + Q\Delta t)\\ &= P(t) + P(t)Q\Delta t\\ \frac{P(t+\Delta t) - P(t)}{(t+\Delta t) - t}&= P(t)Q\\ \underset{t\rightarrow0}{\Longrightarrow} \frac{dP(t)}{dt} &= P(t)Q \end{aligned}$$ここで、この微分方程式は簡単に解くことができる。
- 一般の $n\times n$ 行列 $A$ に対し、$e^{tA} = \sum_{k=0}^{\infty}\frac{t^k}{k!}A^k$
- 上の式で $t$ に対して微分すると、$\frac{d}{dt}e^{tA} = Ae^{tA}$
- 先の式と照らし合わせると、以下を得る。 $$P(t) = e^{Qt}$$
対角化と $P(t)$ の解¶
行列 $Q$ が 異なる固有値を $n(=4)$ 個持つ（対角化できる）とすると、正則行列 $U$ を用いて以下のように対角化できる。
$$\Lambda = U^{-1}QU\Longleftrightarrow Q=U\Lambda U^{-1}$$これを先の式に代入すると、
$$ \begin{aligned} P(t) &=\exp\left(Qt\right)\\ &= \sum_{k=0}^{\infty}\frac{t^k}{k}Q^k\\ &= \sum_{k=0}^{\infty}\frac{t^k}{k}\left(U\Lambda U^{-1}\right)^k\\ &= U\left(\sum_{k=0}^{\infty}\frac{t^k}{k}\Lambda^k\right)U^{-1}\\ &= U\exp\left(\Lambda t\right) U^{-1}\\ &= U \text{diag}\left\{\exp(\lambda_1t), \exp(\lambda_2t),\cdots,\exp(\lambda_nt)\right\}U^{-1} \end{aligned} $$よって、$\Sigma$ が消えて簡単な形で表せるようになった。
各塩基の分布 $\pi(t)$¶
状態 $k$ の初期確率分布が $\pi(0)$ で与えられる時、時刻 $t$ の分布は次式によって計算できる。 $$\pi(t) = \pi(0)P(t)$$
定常分布の時、
$$ \begin{aligned} \pi &= \pi P(t)\\ &= \pi P(0)(I + Qt)\\ &= \pi + \pi Qt\\ \therefore \pi Qt &= 0 \end{aligned} $$最後の式は、$t$ の値に依存せずに成り立つので、$\pi Q=0$
提案されているモデル¶
基本的に $Q$ を定義して、以下の式から $P(t)$ を求める、というのが一般的。
$$ \begin{cases} \begin{aligned} P(t) &=\exp\left(Qt\right) \\ &= U \text{diag}\left\{\exp(\lambda_1t), \exp(\lambda_2t),\cdots,\exp(\lambda_nt)\right\}U^{-1} \end{aligned}\\ Q=U\Lambda U^{-1} \end{cases} $$JC69モデル¶
JC69 (Jukes and Cantor, 1969)は、最も単純な塩基モデルであり、以下の置換速度行列で表される。
$$ Q = \{ q_{i,j} \} = \left[ \begin{array}{rrrr} -3\mu & \mu & \mu & \mu \\ \mu & -3\mu & \mu & \mu \\ \mu & \mu & -3\mu & \mu \\ \mu & \mu & \mu & -3\mu \end{array} \right] $$ここで、$Q$ の固有方程式を解くと、
$$ \begin{aligned} \det\left(Q-\lambda I\right) &= \mu\left(x^4 - 6x^2 + 8x -3\right), \qquad \left(x = (-3-\lambda/\mu)\right)\\ &= (x-1)(x^3 + x^2 - 5x + 3)\\ &= (x-1)(x-1)(x^2 + 2x -3)\\ &= (x-1)(x-1)(x-1)(x+3) \end{aligned} $$したがって、以下の対角化の結果（の一つ）が得られる。
$$ Q = \left( \begin{array}{rrrr} -1 & -1 & -1 & 1 \\ 0 & 0 & 1 & 1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array} \right) \left( \begin{array}{rrrr} -4\mu & 0 & 0 & 0 \\ 0 & -4\mu & 0 & 0 \\ 0 & 0 & -4\mu & 0 \\ 0 & 0 & 0 & 0 \end{array} \right) \left(\frac{1}{4} \left( \begin{array}{rrrr} -1 & -1 & -1 & 3 \\ -1 & -1 & 3 & -1 \\ -1 & 3 & -1 & -1 \\ 1 & 1 & 1 & 1 \end{array} \right) \right) $$In [2]:
class JC69():
def __init__(self):
self.substitution_rate_matrix = np.array([
[-3, 1, 1, 1],
[ 1,-3, 1, 1],
[ 1, 1,-3, 1],
[ 1, 1, 1,-3],
])
self.eigvals, self.U = np.linalg.eig(self.substitution_rate_matrix)
self.U_inv = np.linalg.inv(self.U)
def Q(self, lamda):
""" Substitution Rate Matrix """
return lamda*self.substitution_rate_matrix
def P(self, t, lamda):
Lambda = np.diag(np.exp(t*lamda*self.eigvals))
return np.dot(self.U, np.dot(Lambda, self.U_inv))
def p0(self, t, lamda):
return self.P(t, lamda)[0][0]
def p1(self, t, lamda):
return self.P(t, lamda)[0][1]
In [3]:
time = np.linspace(0, 2, 10000)
In [4]:
model = JC69()
In [5]:
plt.plot(np.vectorize(model.p0)(time, lamda=1), color="red", label="$p_0$")
plt.plot(np.vectorize(model.p1)(time, lamda=1), color="blue", label="$p_1$")
plt.axhline(1/4, linestyle="--", color="black")
plt.legend(), plt.grid(), plt.xticks([0,2500,5000,7500,10000], [0,0.5,1,1.5,2])
plt.show()
