JC69の距離計算
1.1で、単純なp-距離ではなく、多重置換の補正が必要であることが述べられています。JC69を利用した多重置換の補正 を入れた距離推定を行います。
1. 距離計算（１）
JC69の瞬間置換行列から、ある塩基が単位時間当たりにそれ以外の塩基に変化する速度がλとなっており、変わりうる塩基は３種類あるので経過時間\(t\)に対する進化距離\(d\)は\(3λt\)となります。この進化距離を推定します。遷移確率行列\(P\)の1行めに注目します。この行の要素は、初期状態で塩基Tであったものは、時間\(t\)経過したのちに変化せずにTである確率（1行1列の要素\(P[1][1]\)）とC, A, Gに置換されている確率(１行目の２、３、４列の要素\(P[1][2],P[1][3],P[1][4]\)を表しています。
\(P[1][1]=p_{TT}(t)=\frac{1}{4}+\frac{3}{4}e^{-4λt} \)\(P[1][2]=p_{TC}(t)=\frac{1}{4}-\frac{1}{4}e^{-4λt} \)
\(P[1][2]=p_{TA}(t)=\frac{1}{4}-\frac{1}{4}e^{-4λt} \)
\(P[1][2]=p_{TG}(t)=\frac{1}{4}-\frac{1}{4}e^{-4λt} \)
これを全て足すと1になります。他の行も、それぞれ初期状態C(２行目), A（3行目), G(4行目)から時間\(t\)経過した後の状態に対する確率を表しており、各行の対角要素(=塩基が置換していない確率)は\(P[1][1]\)に一致しており、非対角要素はそれぞれ列に対応する塩基に置換される確率を表しており、上の\(P[1][2](=P[1][3]=P[1][4])\)に一致しています。従って、任意のサイトについて、時間\(t\)経過した後にそのサイトが他の塩基に置換されている確率\(p\)は次のようになります。
\(p=1-(\frac{1}{4}+\frac{3}{4}e^{-4λt})=3(\frac{1}{4}-\frac{1}{4}e^{-4λt}) = 3(\frac{1}{4}-\frac{1}{4}e^{-\frac{4}{3}d})\)
となります。この\(p\)を配列相違度\(\widehat{p}\)(\(=p\)-距離）に等しいと置いて\(d=3λt\)の推定値\(\widehat{d}\)を求めます。まず、
\(\widehat{p} = 3(\frac{1}{4}-\frac{1}{4}e^{-\frac{4}{3}d})\)
とおき、これを変形します。
\(1-\frac{4}{3}\widehat{p} = e^{-\frac{4}{3}d} \)
両辺の対数をとって変形すると
\(\widehat{d}=\frac{3}{4}\mathrm{log}(1-\frac{4}{3}\widehat{p})\)
となります。これが(1.6)式です。2. 距離計算（２）
1.4.1では、この距離の導出が最尤法で行われています。maximaを使って、この距離の最尤推定を行ってみましょう。今、長さ\(n\)の二本の塩基配列のアラインメントを考えます（gapは含まれていないものとします)。その中で差異のあるサイトの数が\(x\)であるとします。(1.3)より、任意のサイトで塩基が一致している確率\(p_{0}(t)\)と、差異のある確率\(3p_{1}(t)\)は、次の式で表されます。
\(p_{0}(t)=\frac{1}{4}+\frac{3}{4}\mathrm{e}^{-\frac{4}{3}d}\)
\(3p_{1}(t)=\frac{3}{4}-\frac{3}{4}\mathrm{e}^{-\frac{4}{3}d}\)
ここで、後の計算のため\(λt\)を\(d/3\)に変更しています。これを用いると、上記の配列アラインメントを観測する確率(=尤度）は、二項確率で表すことができます。
\(L(d:x)=\frac{n!}{x!(n-x)!}(\frac{3}{4}-\frac{3}{4}\mathrm{e}^{-4d/3})^{x}(\frac{1}{4}+\frac{3}{4}\mathrm{e}^{-4d/3})^{n-x}\)
これが(1.42)です。テキストでは、サイト•パターンを区別する形に拡張した(1.43)も与えられていますが、ここでは上の式で最尤法を行います。この確率が最大になるように\(d\)を推定するわけですが、その前に上の式を次のように変形します。まず、\(d\)に依存しない定数は計算に影響しないので、省きます。また、尤度の対数をとった対数尤度を最尤推定に用います。対数尤度関数は次のようになります。
\(\ell(d:x)=x\mathrm{log}(\frac{3}{4}-\frac{3}{4}\mathrm{e}^{-4d/3})+(n-x)\mathrm{log}(\frac{1}{4}+\frac{3}{4}\mathrm{e}^{-4d/3})\)
サイト•パターンを区別した(1.44)とは定数項だけ異なっており、尤度関数としては同じ式になります。この尤度関数をmaximaで定義しましょう。
lkf(d):=x*log((3/4)-(3/4)*exp(-4*d/3))+(n-x)*log((1/4)+(3/4)*exp(-4*d/3));
これを最大化する\(d\)を求めるには、この式を\(d\)で微分したものを0とおいて(\(\frac{\mathrm{d}\ell}{\mathrm{d}d}= 0\))、それを満たす\(d\)を求めます。微分には関数diffを用います。
diff(lkf(d),d);
すると次のように出力されます。
\(\frac{\%\mathrm{e}^{-\frac{4d}{3}}x}{\frac{3}{4}-\frac{3\%\mathrm{e}^{-\frac{4d}{3}}}{4}}-\frac{\%\mathrm{e}^{-\frac{4d}{3}}(n-x)}{\frac{3\%\mathrm{e}^{-\frac{4d}{3}}}{4} + \frac{1}{4}}\)
この式を0とおき、\(3exp(-4d/3)/4\)をaとおいてまとめなおすと
\((a +(1/4))x - ((3/4)-a)(n-x) = 0\)
となります。solveで\(a\)をもとめてみましょう。
partfrac(solve([(a +(1/4))*x - ((3/4)-a)*(n-x) = 0], a), n);
ここで、partfracは、部分分数展開のためのmaximaの関数です。すると
\([a = \frac{3}{4}-\frac{x}{n}]\)
これが\(a\)の値となります。\(a = 3\mathrm{exp}(-4d/3)/4\)とおいたので、\(d\)について解くために、この式に４/3をかけた後に、対数をとり、さらに-3/4をかけます。
(-3/4)*log((4/3)*(3/4-x/n));
すると、
\(-\frac{3\mathrm{log}\left[\frac{4\left(\frac{3}{4}-\frac{x}{n}\right)}{3}\right]}{4}\)
となり、(1.6)が得られました。
前ページ JC69の定常分布 | 次ページ K80