K80
1. K80のQ
JC69 は、どの塩基の間も同じ速度で置換が起きることを仮定していました。しかし、トランジションはトランスバージョンより生じやすいことから、その置換速度の違いを考慮したモデルがK80 (Kimura, 1980)です。K80の瞬間置換速度行列\(Q\)が、テキストの(1.8)式で与えられています。
\begin{equation} Q = \{ q_{i,j} \} = \left[ \begin{array} {rrrr} -(\alpha+2\beta) & \alpha & \beta & \beta \\ \alpha & -(\alpha+2\beta) & \beta & \beta \\ \beta & \beta & -(\alpha+2\beta) & \alpha \\ \beta & \beta & \alpha & -(\alpha+2\beta) \end{array} \right] \end{equation}JC69と同様に、この行列の各行、各列は、T,C,A,Gの４つの塩基に対応しています。 また、この\(Q\)と遷移確率行列の関係も、JC69と同様に次のようになります。
\begin{equation} P(t) = \mathrm{e}^{Qt} = \left[ \begin{array} {rrrr} p_{0}(t) & p_{1}(t) & p_{2}(t) & p_{2}(t) \\ p_{1}(t) & p_{0}(t) & p_{2}(t) & p_{2}(t) \\ p_{2}(t) & p_{2}(t) & p_{0}(t) & p_{1}(t) \\ p_{2}(t) & p_{2}(t) & p_{1}(t) & p_{0}(t) \end{array} \right] \end{equation}これがテキストの(1.9)です。また、この遷移確率行列の要素は、
\(p_{0}(t) = \frac{1}{4} + \frac{1}{4}\mathrm{e}^{-4\beta t}+\frac{1}{2}\mathrm{e}^{-2(\alpha+\beta)t} = \frac{1}{4} + \frac{1}{4}\mathrm{e}^{-4d/(\kappa+2)}+\frac{1}{2}\mathrm{e}^{-2d(\kappa+1)/(\kappa+2)}\)
\(p_{1}(t) = \frac{1}{4} + \frac{1}{4}\mathrm{e}^{-4\beta t}-\frac{1}{2}\mathrm{e}^{-2(\alpha+\beta)t} = \frac{1}{4} + \frac{1}{4}\mathrm{e}^{-4d/(\kappa+2)}-\frac{1}{2}\mathrm{e}^{-2d(\kappa+1)/(\kappa+2)}\)
\(p_{2}(t) = \frac{1}{4} - \frac{1}{4}\mathrm{e}^{-4\beta t} = \frac{1}{4} -\frac{1}{2}\mathrm{e}^{-4d(\kappa+1)/(\kappa+2)}\)
となります。これがテキストの(1.10)です。
2. maxima による導出
K80の\(Q\)から、遷移確率行列(1.9)やその要素(1.10)をmaximaで導きましょう。2-1. Qtの固有値、固有ベクトルの計算
まず、\(Qt\)をmaxima上で次のように設定します。
Qt:matrix([-(a+2*b),a,b,b],[a,-(a+2*b),b,b],[b,b,-(a+2*b),a],[b,b,a,-(a+2*b)]);
ここでは簡単のため、\(\alpha t\)を\(a\)で、\(\beta t\)を\(b\)で表しています。すると、
\begin{equation} \left[ \begin{array} {cccc} -2 b - a & a & b & b \\ a & -2 b - a & b & a \\ b & b & -2 b - a & a \\ b & b & a & -2 b - a \end{array} \right] \end{equation}が出力されます。この行列\(Qt\)の固有値問題を関数eigenvectorsを使って解きます。
ev:eigenvectors(Qt);
evにリストとして、固有値問題の解が代入されています。
ev[1];
と入力すると、
[[-4 b, -2 b - 2a, 0],[1, 2, 1]]
と出力されます。これは、Qtが縮退のない固有値, \(-4b\), 2重縮退した固有値 \(-2b-2a\) と、縮退のない固有値 0 を持つことを示しています。
ev[2];
と入力すると、固有値\(-4b\)に対する固有ベクトルと、\(-2b-2a\)に対する２つの固有ベクトルが得られます。
[[[1,1,-1,-1]],[[1,-1,0,0],[0,0,1,-1]]]
また、
ev[2][1][1];
ev[2][2][1];
ev[2][2][2];
として、個々の固有ベクトルを取り出すことができます。この３つのベクトルが線形独立であることは視察からも明らかですが、maximaで確認しましょう。
rank(transpose(matrix(ev[2][1][1], ev[2][2][1], ev[2][2][2])));
transposeは、ここでは必ずしも必要ありません。 この処理の結果は、３と出力され、これら3つの固有ベクトルは1次独立あることがわかります。また、
ev[2][1][1].ev[2][2][1];
ev[2][1][2].ev[2][2][2];
は、いずれも\(0\)になり、異なる固有値に属す固有ベクトルは直交することがわかります。
今回も固有値0に対する固有ベクトルが出力されていません。JC69で説明したのとは異なる方法で、固有値0に対する固有ベクトルをもとめてみましょう。JC69もK80も各行の要素の総和は0になっています。そのため、\(Q\)と(1,1,1,1)の積は(0,0,0,0)になり、(1,1,1,1)が0の固有ベクトルになっていることがわかります。(1,1,1,1)と他の固有値に対する固有ベクトルとの内積をとってみましょう。
ev[2][1][1].[1,1,1,1];
ev[2][2][1].[1,1,1,1];
ev[2][2][1].[1,1,1,1];
は、いずれも0になり、(1,1,1,1)は、他の固有値に対する固有ベクトルとは直交していることがわかります。
2-2. Qtの対角化
JC69と同様に、上で求めた固有ベクトルを列とした行列を作成します。最初に、固有値-4*xに対応する固有ベクトルを並べ、次に固有値0に対する固有ベクトルを並べます。
U:transpose(matrix(ev[2][1][1],ev[2][2][1],ev[2][2][2],[1,1,1,1]));
また、Uの逆行列UIをinvertを使って次のように求めておきます。
UI:invert(U);
ここで、
ratsimp(UI.Qt.U);
を計算すると、固有値を対角要素とした次の対角行列が出力されます。
\begin{equation} \left[ \begin{array} {rrrr} -4 b & 0 & 0 & 0 \\ 0 & -2 b - 2 a & 0 & 0 \\ 0 & 0 & -2 b - 2 a & 0 \\ 0 & 0 & 0 & 0 \end{array} \right] \end{equation}2-3. 遷移確率行列を求める
まず、対角要素に\(exp(-4b), exp(-2b-2a), exp(-2b-2a), 1(=exp(0))\)を持つ対角行列\(M1\)を作成します。
M1:matrix([exp(-4*b),0,0,0],[0,exp(-2*b-2*a),0,0],[0,0,exp(-2*b-2*a),0],[0,0,0,1]);
ここで、M1;とウィンドウに入力すると、\(M1\)が次のようになっていることが確認できます。
\begin{equation} \left[ \begin{array} {rrrr} \%\mathrm{e}^{-4 b} & 0 & 0 & 0 \\ 0 & \%\mathrm{e}^{-2 b - 2 a} & 0 & 0 \\ 0 & 0 & \%\mathrm{e}^{-2 b - 2 a} & 0 \\ 0 & 0 & 0 & 1 \end{array} \right] \end{equation}次に、JC69のページの1の説明に従い、遷移確率行列\(P\)を次のようにして求めます。
P:U.M1.UI;
すると\(P\)として次の行列がウィンドウ上に表示されます。
\begin{equation} \left[ \begin{array} {cccc} \frac{\%\mathrm{e}^{-4 b}}{4} + \frac{\%\mathrm{e}^{-2 b-2 a}}{2} + \frac{1}{4} & \frac{\%\mathrm{e}^{-4 b}}{4}-\frac{\%\mathrm{e}^{-2 b-2 a}}{2} + \frac{1}{4}& \frac{1}{4}-\frac{\%\mathrm{e}^{-4 b}}{4} & \frac{1}{4}-\frac{\%\mathrm{e}^{-4 b}}{4} \\ \frac{\%\mathrm{e}^{-4 b}}{4}-\frac{\%\mathrm{e}^{-2 b-2 a}}{2} + \frac{1}{4} &\frac{\%\mathrm{e}^{-4 b}}{4} + \frac{\%\mathrm{e}^{-2 b-2 a}}{2} + \frac{1}{4} & \frac{1}{4}-\frac{\%\mathrm{e}^{-4 b}}{4} & \frac{1}{4}-\frac{\%\mathrm{e}^{-4 b}}{4} \\ \frac{1}{4}-\frac{\%\mathrm{e}^{-4 b}}{4} & \frac{1}{4}-\frac{\%\mathrm{e}^{-4 b}}{4} & \frac{\%\mathrm{e}^{-4 b}}{4} + \frac{\%\mathrm{e}^{-2 b-2 a}}{2} + \frac{1}{4} & \frac{\%\mathrm{e}^{-4 b}}{4}-\frac{\%\mathrm{e}^{-2 b-2 a}}{2} + \frac{1}{4} \\ \frac{1}{4}-\frac{\%\mathrm{e}^{-4 b}}{4} & \frac{1}{4}-\frac{\%\mathrm{e}^{-4 b}}{4} & \frac{\%\mathrm{e}^{-4 b}}{4}-\frac{\%\mathrm{e}^{-2 b-2 a}}{2} + \frac{1}{4} & \frac{\%\mathrm{e}^{-4 b}}{4} + \frac{\%\mathrm{e}^{-2 b-2 a}}{2} + \frac{1}{4} \end{array} \right] \end{equation}aとbは表記を簡単するために\(\alpha t\)と\(\beta t\)を置き換えたものでしたから、得られた行列の要素は式(1.10)に対応することがわかります。
また、各行の要素は、その行に対応する塩基から、列に対応する塩基への遷移確率を表しています。例えば1行２列は、塩基Tから塩基Cへの遷移確率になります。そこで各行の遷移確率の和をとると、行に対応する塩基から、可能な全ての塩基への遷移確率の和をとることになるので、１になります。実際に、
ratsimp(P[1,1]+P[1,2]+P[1,3]+P[1,4]);
ratsimp(P[2,1]+P[2,2]+P[2,3]+P[2,4]);
ratsimp(P[3,1]+P[3,2]+P[3,3]+P[3,4]);
ratsimp(P[4,1]+P[4,2]+P[4,3]+P[4,4]);
を実行すると、いずれも1になります。
前ページ JC69の距離計算 | 次ページ K80の定常分布