JC69

JC69 (Jukes and Cantor, 1969)は、最も単純な塩基置換モデルです。JC69モデルは、次の瞬間置換速度行列(instantaneous substitution rate matrix)で表現されます。テキストの(1.1)式です。

この行列の各行、各列は、T,C,A,Gの４つの塩基に対応しています。表1.1や図1.2には、塩基置換の制約を変えた他のモデルが紹介されています。

1. 遷移確率行列と\(Q\)の関係

一方、\(p_{i,j}(t)\)を、塩基iが時間t経過した時に塩基jに置換される確率とします。これを遷移確率とよびます。

詳細は1.5節で触れられていますが、遷移確率行列を\(P(t)\)={\(p_{i,j}(t)\)}とすると \begin{equation} \frac{\mathrm{d}P(t)}{\mathrm{d}t}=P(t)Q \end{equation} となります。これが(1.52)式で、この解が(1.2)式 \begin{equation} P(t)=\mathrm{e}^{Qt} \end{equation} となります。行列の指数関数は次のようになります。 \begin{equation} \mathrm{e}^{Qt}=I + Qt + \frac{1}{2!}(Qt)^2 + \frac{1}{3!}(Qt)^3 + ... \end{equation} \(Q\)の固有値問題を解くことができ、次のように対角化できたとします。 \begin{equation} Q = UΛU^{-1} \end{equation} \(Λ\)は、\(Q\)の固有値{\(v_{1}, v_{2}, v_{3}, v_{4}\)}を対角成分とする対角行列です。各固有値に対応する固有ベクトルが\(u_{1}, u_{2}, u_{3}, u_{4}\)とすると、\(U\) = {\(u_{1}, u_{2}, u_{3}, u_{4}\)}であり、その逆行列が\(U^{-1}\)です。これを行列の指数関数の式に代入します。すると、 \begin{equation} \mathrm{e}^{Qt} = UU^{-1} + UΛU^{-1}t + \frac{1}{2!}UΛU^{-1}UΛU^{-1}t^{2} + \frac{1}{3!}UΛU^{-1}UΛU^{-1}UΛU^{-1}t^{3} + ... \end{equation} \begin{equation} = UU^{-1} + UΛU^{-1}t + \frac{1}{2!}UΛ^{2}U^{-1}t^{2} + \frac{1}{3!}UΛ^{3}U^{-1}t^{3} + ... \end{equation} \begin{equation} =U(I + Λt + \frac{1}{2!}Λ^{2}t^{2} + \frac{1}{3!}Λ^{3}t^{3} + ...)U^{-1} \end{equation} 対角行列のべき乗は、対角要素のべき乗を要素とした対角行列となることから、上式は次のようになります。 \begin{equation} \mathrm{e}^{Qt}=U \left[ \begin{array} {rrrr} \mathrm{e}^{v_{1}t} & 0 & 0 & 0 \\ 0 & \mathrm{e}^{v_{2}t} & 0 & 0 \\ 0 & 0 & \mathrm{e}^{v_{3}t} & 0 \\ 0 & 0 & 0 & \mathrm{e}^{v_{4}t} \end{array} \right] U^{-1} \end{equation} つまり、\(Q\)の固有値と固有ベクトルが得られれば、遷移確率行列\(P(t)\)は上式の形で簡単に求めることができ、それが(1.3)式で表されるものになるわけです。 \begin{equation} P(t) = \mathrm{e}^{Qt} = \left[ \begin{array} {rrrr} p_{0}(t) & p_{1}(t) & p_{1}(t) & p_{1}(t) \\ p_{1}(t) & p_{0}(t) & p_{1}(t) & p_{1}(t) \\ p_{1}(t) & p_{1}(t) & p_{0}(t) & p_{1}(t) \\ p_{1}(t) & p_{1}(t) & p_{1}(t) & p_{0}(t) \end{array} \right] \end{equation} ここで、\(p_{0}(t) = \frac{1}{4} + \frac{3}{4}\mathrm{e}^{-4λt}, p_{1}(t) = \frac{1}{4} - \frac{1}{4}\mathrm{e}^{-4λt}\)となります。 この計算をMaximaで確認しましょう。

2. maxima による導出

アプリケーションフォルダ内のwxMaximaのアイコンをダブルクリックすると、wxMaximaのウィンドウがたちあがります。 このウィンドウ上で、操作をおこないます。

2-1. Qtの固有値、固有ベクトルの計算

まず、行列\(Qt\)をmatrixを使って次のように設定します。

Qt:matrix([-3*x,x,x,x],[x,-3*x,x,x],[x,x,-3*x,x],[x,x,x,-3*x]);

ここでは簡単のため、\(λt\)を\(x\)で表しています。[と]で囲まれた部分で行列の行を指定しています。Qtの後ろのコロン(:)は代入を意味します。maximaでは、代入にコロンが使用されます。 また、コマンドの行末にセミコロン(;)をつける必要があります。 入力が終わったら、Shiftキーを押したままで、Enterキーを押すとコマンドが実行され、次の行列がウィンドウ上に表示されます。

この行列\(Qt\)の固有値問題を関数eigenvectorxを使って解きます。

ev:eigenvectors(Qt);

と入力して、Shiftキーを押したまま、Enterキーを押します。eigenvectorsは（）の中で指定された行列の固有値問題を解く関数で、固有値と固有ベクトルが得られます。ここでは、コロンを使ってその結果をevという変数に代入しています。evは、複数のデータを組にしたリストと呼ばれる構造を持っています。（注：古いバージョンのmaximaのeigenvectorsの出力のリスト構造は少し違っているようです。）

ev[1];

と入力すると、

[[-4 x, 0],[3, 1]]

と出力されます。これは、Qtが３重縮退した固有値 -4 x と、縮退のない固有値 0　を持つことを示しています。

ev[2];

と入力すると、固有値-4 x に対する3つの固有ベクトルが得られます。

[[[1,0,0,-1],[0,1,0,-1],[0,0,1,-1]]]

ev[2][1][1], ev[2][1][2], ev[2][1][3]として、個々の固有ベクトルを取り出すことができます。この３つのベクトルが線形独立であることが明らかですが、maximaで確認しましょう。

rank(transpose(matrix(ev[2][1][1], ev[2][1][2], ev[2][1][3])));

transposeは、引数となる行列を転置する関数です。この操作は、ここでは必ずしも必要ありませんが、transposeは後で使用します。 rankは、引数となる行列の階数を出力してくれます。今回の場合、３と出力され、固有値の縮退度と、固有空間の次元が一致していることがわかります。

以前のバージョンのmaximaでは、固有値0に対する固有ベクトル\([1,1,1,1]\)を出力してくれていたのですが、現在のバージョンでは出てこないようなので計算してみましょう。今, 固有値0に対する固有ベクトルが\([a, b, c, d]\)であるとします。これとQtの積が,\([0,0,0,0]\)となるように、\([a, b, c, d]\)を決めます。まず、このベクトルをAとして表します。

A:[a,b,c,d];

次にQtとAの積をとります。

Qt.A;

ここで、QtとAの間の'.'は、ベクトルの内積、行列とベクトル、また行列同士の積の計算を意味します。次の縦ベクトルが出力されます。

となりますが、見難いのでratsimpを用いて簡約化します。

ratsimp(Qt.A);

すると、上と同じものですが、次のように簡約化されて表示されます。

solve([d+c+b-3*a=0,d+c-3*b+a=0,d-3*c+b+a=0,-3*d+c+b+a=0],[a,b,c,d]);

rank(matrix([1,1,1,-3],[1,1,-3,1],[1,-3,1,1],[-3,1,1,1]));

solve([1+c+b-3*a=0,1+c-3*b+a=0,1-3*c+b+a=0],[a,b,c]);

[1,1,1,1].ev[2][1][1];
[1,1,1,1].ev[2][1][2];
[1,1,1,1].ev[2][1][3];

U:transpose(matrix(ev[2][1][1],ev[2][1][2],ev[2][1][3],[1,1,1,1]))

UI:invert(U);

UI.Qt.U;

M1:diagmatrix(4,exp(-4*x));

M1[4][4]:1;

P:U.M1.UI;