Ka/Ks (dN/dS) の計算 - 驚異のアニヲタ社会復帰の予備

MikuHatsune2013-06-21

ある塩基配列の変異パターンを解析しているのだが、それでKa/Ks (dN/dS)比というモデルを勉強した。
しかしこれには計算法がたくさんあるらしい。
サンプルデータは http://evolution.genetics.washington.edu/book/primates.dna をコピペして使う。

library(phangorn)
library(seqinr)
primates <- read.phyDat("http://evolution.genetics.washington.edu/book/primates.dna", format="phylip", type="DNA")
tree0 <- NJ(dist.ml(primates))
plot(tree0)

　
Miyata & Yasunaga. J Mol Evol. 1980 Sep;16(1):23-36.
Nei & Gojobori. Mol Biol Evol. 1986 Sep;3(5):418-26.
Sqdif Plotで両者計算できる。計算結果もコピペで取ってこれるが、複数の pairwise alignment をやったときには一括ダウンロードできるようにしてほしい。
　
Li. J Mol Evol. 1993 Jan;36(1):96-9.
seqinrの kaks 関数に実装されているアルゴリズム。
解説
 縮退 (degenerate)というコドンの特徴があるが、これを
$L_0$ ：すべてのありうる変化が非同義
$L_2$ ：変化のうち 1/3 が同義
( $L_3$ ：イソロイシンのコドン3番目の場合。 $L_2$ とみなしてしまう)
$L_4$ ：すべてのありうる変化が同義
$A_i$ ： $L_i;\hspace{3}i=\{0,2,4\}$ に対応するTransition (purine/pyrimidine が変わらない変異)の数
$B_i$ ： $L_i;\hspace{3}i=\{0,2,4\}$ に対応するTransversion (purine/pyrimidine が変わる変異)の数
と定義して、
$Ks=\frac{L_2A_2+L_4A_4}{L_2+L_4}+B_4$
$Ka=A_0+\frac{L_0B_0+L_2B_2}{L_0+L_2}$
を計算する。

# コドンのため3の倍数の塩基長になっていないとエラー
k0 <- kaks(as.alignment(length(primates), seq=apply(as.character(primates)[,1:231], 1, paste, collapse="")))

　
最尤法 Mol Biol Evol. 1998 Dec;15(12):1600-11.
$i$ 番目のコドンと $j$ 番目のコドンについて、 $q_{i,j}$ を
$q_{i,j}=0$ ：if i and j dier in more than one position
$q_{i,j}=\pi_j$ ：for synonymous transversion
$q_{i,j}=\pi_j\kappa$ ：for synonymous transition
$q_{i,j}=\pi_j\omega$ ：for nonsynonymous transversion
$q_{i,j}=\pi_j\omega\kappa$ ：for nonsynonymous transition
とする。ここで、 $\omega$ 、は　Ka/Ks, $\kappa$ は Transition と Transversion の比、 $\pi_j$ はコドン $j$ の平衡頻度である。 $\omega$ が1かうんたらかは前も言ったとおり。
phangornパッケージの pml 関数でできる。optim.pml で type を codon_n の何に指定するかで、 $\omega$ と $\kappa$ を1に固定するか自由に推定するか決められる。
サンプルはこちらを参考にしたが、サンプルがないわコピペ仕様になってないわで苦労した。

dat <- phyDat(as.character(primates), type="CODON")
fit0 <- pml(tree0, dat)
fit00 <- optim.pml(fit0, model="codon0", control=pml.control(trace=0))
fit1 <- optim.pml(fit0, model="codon1", control=pml.control(trace=0))
fit2 <- optim.pml(fit0, model="codon2", control=pml.control(trace=0))
fit3 <- optim.pml(fit0, model="codon3", control=pml.control(trace=0))
anova(fit00, fit2, fit3, fit1)

# 同じ結果
Likelihood Ratio Test Table
  Log lik. Df Df change Diff log lik. Pr(>|Chi|)    
1  -2905.6 25                                       
2  -2385.7 26         1       1039.65     <2e-16 ***
3  -2292.6 26         0        186.27     <2e-16 ***
4  -2291.4 27         1          2.39     0.1218    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

pairwise で比較して Ka/Ks を求めたいのだが、なぜかできない。
というわけでゴリ押ししようと思ったら、pairwise では optEdge=FALSE にしないと unroot がどうたらこうたら言われてエラーになる。
NJ法のところもエラーになってRがオシャカになるので、他の方法で距離行列を求める。

p0 <- as.character(primates)[, 1:231]
res <- diag(0, nrow(p0))
dimnames(res) <- list(rownames(p0), rownames(p0))
cmb <- combn(nrow(p0), 2)
for(i in seq(ncol(cmb))){
	tmp <- p0[cmb[, i], ]
	tree0 <- as.phylo(hclust(dist.ml(as.phyDat(tmp)))) # NJ だと壊れる
	dat <- phyDat(p0, type="CODON")
	fit0 <- pml(tree0, dat)
	# pairwise だと optEdge=FALSE にしないといけないっぽい
	pml1 <- optim.pml(fit0, model="codon1", optEdge=FALSE) # k and omega are not fixed to 1
	res[cmb[2, i], cmb[1, i]] <- pml1$dnds
}

系統樹を作成するときにもいろいろなアルゴリズムが使える。

# Maximau likelifood
# pairwise だと optEdge=FALSE にしないといけないっぽい
pml1 <- optim.pml(fit0, model="codon1", optEdge=FALSE) # k and omega are not fixed to 1
pml1$dnds

# Jukes-Cantor (starting tree from NJ)
fitJC <- pml(tree0, dat)

# optimize edge length parameter
fitJC <- optim.pml(fitJC, model="codon1", optEdge=FALSE)

# JC + Gamma + I - model
fitJC_GI <- update(fitJC, k=4, inv=0.2)

# optimize shape parameter + proportion of invariant sites
fitJC_GI <- optim.pml(fitJC_GI, optGamma=TRUE, optInv=TRUE, optEdge=FALSE)

# GTR + Gamma + I - model
fitGTR <- optim.pml(fitJC_GI, optBf=TRUE, optQ=TRUE, model="codon1", optEdge=FALSE)

# Pamilo-Bianchi-Li
k0 <- kaks(as.alignment(2, seq=apply(p0, 1, paste, collapse="")))
k0$ka / k0$ks