分子進化解析をしているときに、PAMLを使って自然淘汰の検出を行うことがある。特定の枝について（大抵の場合は、注目している系統や種、分類群を指定することになる）どの塩基に正の自然淘汰の影響があるのかを検出したい時に、branch-site test of positive selectionをすることが多い。テストをする際に、PAML (Phylogenetic Analysis by Maximum Likelihood) のcodemlというプログラムを用いる (1)。解析にはnull modelとalternative modelを記述するためのコントロールファイルを準備する必要があるのだが、これをどう準備すればいいのかをよく忘れる。その度にマニュアルやYangの本を引っ張り出して確認して思い出す。今回行ったので、ここに備忘録としてまとめておきたいと思う。
内容については、textや論文、PAMLのドキュメントなどを参照しているが、間違いを含んでいる可能性がある。
PAML
PAML (Phylogenetic Analysis by Maximum Likelihood)は、最尤法を使って分子進化の解析をすることができるプログラム集で、Ziheng Yangによって作成されている。進化系の論文では、分子進化速度を推定したり、モデルに基づいて自然淘汰の検出を行ったりする際によく用いられている。いくつかのプログラムからなるが、DNA配列の解析によく用いられるのは、その中でもcodemlというプログラムだと思う。Windows、 Unix系OS、MacOSなどで使うことができる。これまでの経験では、たまにOSによって準備したファイルの不備に対して返してくれるエラーメッセージがあったりなかったりといった細かな違いがあった。自分は主にMacOSX版を使っている。
帰無仮説となるモデルと対立仮説となるモデルでの解析を行ったのち、likelihood ratio testをしてpositive selectionの有無を検定することができる。例えば正の自然淘汰が枝や、サイトに働いているというモデルを、帰無仮説のモデルと比較する。ある枝に対して自然淘汰が働いているとする場合はBranch model、特定のサイトに働いているとする場合はsite modelによる検定となる。
一方で、branch modelやsite modelは正の自然淘汰を検出するためにはやや保守的な仮定を置くことになるかもしれない。branch modelの場合はある特定の枝で、遺伝子の塩基配列全体にわたってdN/dS比が上昇していることを仮定することになるが、ある特定のサイトのみが自然淘汰の対象となっている場合（こちらの方が多くの場合自然に感じられる。そのほかのサイトはpurifying selectionの対象になっているだろう）、dN/dSを遺伝子配列全体に渡って推定すると、正の自然淘汰を受けているサイトが埋没してしまいうまく検出できないかもしれない。同様に、site modelでは、系統樹のどの枝でも、正の自然淘汰の対象になっていると仮定する。これは、host-pathogenの軍拡競争的な淘汰が働いている場合は当てはまるだろう。しかし、そのほかの大部分の遺伝子で、系統樹のうちのいくつかの種では自然淘汰が働いていて、そのほかの種では特に働いていないということは普通にあるだろうが、そういった場合はsite modelでは保守的すぎることになるだろう。
branch-site testは、Yang et al., 2005、Zhang et al., 2005で示されたモデルで、Yang and Nielsen (2002)のモデルの改良版だ。alternative modelとなるmodel Aと、branch-siteモデルは、dN/dSがサイトごとにも、枝によってもばらついていることを仮定する。そして、特定の枝（foregroundの枝）の特定のサイトのみ自然淘汰の影響下にあると仮定してそれを検出しようとする。
プログラムのドキュメントには、PAMLでは何が
できないか (What PAML programs cannot do)、についても書いてあり、1) 配列のアラインメント、2) 遺伝子予測、3) 大規模データからのtree探索、用のプログラム
ではない旨が書かれている。
codeml
protein coding regionを対象とした塩基配列の解析である場合、プログラムのcodemlを使うことになる。codemlはcodon substitution modelsから来ているらしい。必要になるのは、1) 配列のphylipファイル、2) 系統関係を表すnewick形式のファイル、そして3) パラメータなどを記したコントロールファイルの三つとなる。これらが同じディレクトリにありコントロールファイルの名前がcodeml.ctlとなっていれば、シンプルに、
> codeml
とコマンドを打てばプログラムがスタートする。もし、コントロールファイルの名前を変更していたら（例えば、codeml.MA1.ctlなど）、引数に指定して
> codeml codeml.MA1.ctl
とすれば良い。phylipファイルやnewickファイルはコントロールファイルの中に指定する。あとは、アウトプットファイルの名前も指定する必要がある。
よくあるのが、配列の中にストップコドンがある場合。この場合、ランが始まったのちに警告メッセージが出るようになっている（ver.4.9で確認)。この場合はエンターキーで続行となる。ちなみに、ソースコードを編集すると、警告メッセージなしで先に進めるようにできる。検索するといくつか結果が出てくる。
codeml.ctlファイルの中身は、以下のようなパラメータの指定で構成されている。
seqfile = stewart.aa * sequence data filename
treefile = stewart.trees * tree structure file name
outfile = mlc * main result file name
noisy = 9 * 0,1,2,3,9: how much rubbish on the screen
verbose = 1 * 0: concise; 1: detailed, 2: too much
runmode = 0 * 0: user tree; 1: semi-automatic; 2: automatic
* 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise
seqtype = 2 * 1:codons; 2:AAs; 3:codons-->AAs
CodonFreq = 2 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table
* ndata = 10
clock = 0 * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
aaDist = 0 * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a
aaRatefile = dat/jones.dat * only used for aa seqs with model=empirical(_F)
* dayhoff.dat, jones.dat, wag.dat, mtmam.dat, or your own
model = 2
* models for codons:
* 0:one, 1:b, 2:2 or more dN/dS ratios for branches
* models for AAs or codon-translated AAs:
* 0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F
* 6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)
NSsites = 0 * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
* 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
* 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
* 13:3normal>0
icode = 0 * 0:universal code; 1:mammalian mt; 2-10:see below
Mgene = 0
* codon: 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff
* AA: 0:rates, 1:separate
fix_kappa = 0 * 1: kappa fixed, 0: kappa to be estimated
kappa = 2 * initial or fixed kappa
fix_omega = 0 * 1: omega or omega_1 fixed, 0: estimate
omega = .4 * initial or fixed omega, for codons or codon-based AAs
fix_alpha = 1 * 0: estimate gamma shape parameter; 1: fix it at alpha
alpha = 0. * initial or fixed alpha, 0:infinity (constant rate)
Malpha = 0 * different alphas for genes
ncatG = 8 * # of categories in dG of NSsites models
getSE = 0 * 0: don't want them, 1: want S.E.s of estimates
RateAncestor = 1 * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)
Small_Diff = .5e-6
cleandata = 1 * remove sites with ambiguity data (1:yes, 0:no)?
* fix_blength = 1 * 0: ignore, -1: random, 1: initial, 2: fixed, 3: proportional
method = 0 * Optimization method 0: simultaneous; 1: one branch a time
* Genetic codes: 0:universal, 1:mammalian mt., 2:yeast mt., 3:mold mt.,
* 4: invertebrate mt., 5: ciliate nuclear, 6: echinoderm mt.,
* 7: euplotid mt., 8: alternative yeast nu. 9: ascidian mt.,
* 10: blepharisma nu.
* These codes correspond to transl_table 1 to 11 of GENEBANK.
これは、元のファイルに添付されていたものである。
seqfile, treefile, outfileのそれぞれのところにシークエンスファイル、系統関係のファイル、結果ファイルの名前を記入する。パスを指定することもできる。
seqtypeは、DNA塩基配列の場合は1を指定する。モデルによって頻繁に指定を変更するところはmodelとNSsitesの項目、あとfix_omega, omegaの項目など。
phylip形式は、以下のようなファイル形式になっている。
6 1527
Sp1 acaacaaatggaacattgaagatcccag...
Sp2 tcaacgaatggatcattgaagatcccag...
Sp3 acaacaaatggaacatagaagctcccag...
Sp4 acaacaaatggaacattgaagatcccag...
Sp5 acaacaaatggaacattgaagatcccaa...
Sp6 acaccaactgcaaccttcaacatgccaa...
この中で、最初の６が配列数、1527が配列長を示す。Species1から6まであるとして、配列名の後にいくつか空白を開けて配列を指定する。
newick形式は、以下のような形式で系統関係を記述する。
(((Sp1, Sp2), Sp3) , (Sp4, Sp1), Sp6);
newickとphylipのファイルの配列名が一致していないとエラーになる。
Branch model
branch modelはdN/dSが系統樹の各枝ごとにばらついているという仮定をとる。このため、特定の枝（特定の系統だったり、種だったり）に自然淘汰が働いている場合にそれを検出しようとするのに適している。mode=1と指定すると、枝ごとにdN/dSを推定してくれる。これがいわゆるfree-ratios modelとなる。
free-ratios modelは枝の数だけdN/dSを推定するのでパラメータがとても多くなる。model=2とすると、枝ごとにいくつかのdN/dSが指定できるようになる。どのように指定するかというと、newick形式のファイルの方で、
(((Sp1 #1, Sp2 #1) #1, Sp3) , (Sp4, Sp1), Sp6);
というように、#1というラベルでSp1と2に続く枝に一つのdN/dSの値を想定し、そのほかの枝に別の値を想定することができる。ラベルを増やして
(((Sp1 #1, Sp2 #1 ) #1, Sp3 #2) #2 , (Sp4, Sp1), Sp6);
といったように指定することもできる。この場合、#1と#2とそれ以外の枝でそれぞれdN/dSの値が異なることになる。
Site model
Site modelでは、dN/dSの値がサイトごとにばらつくと想定する。この場合、model=0を指定して、NSsitesをモデルにより違う値に設定する。自然淘汰の検出をsite modelを使ってやるときには、特に2つの比較をすることが多い。
M1a (neutral) vs M2a (Positive selection):
M1aは中立を仮定したモデルで、帰無仮説のモデルになる。パラメータとしてNSsites=1, model=0を指定する。この場合のfree parameterの数は2になる。
M2aは自然淘汰を仮定したモデルで、alternative modelになる。パラメータとしてNSsites=2, model=0を指定する。free parameterの数は4になる。
このペアの間でLikelihood Ratio Test (LRT)をする場合は、d.f.は2になる（4-2=2）。
M7 (beta) vs M8 (beta&omega)
M7はdN/dSの分布にベータ分布を仮定している。例えばbeta(0.2, 1)の分布だと、０のところが一番多くて、dN/dSが大きくなるにつれてなだらかに少なくなるような分布になる。パラメータとして、NSsites=7, model=0を指定する。 free parameterは2となる。
M8は対立仮説のモデルとなる。ベータ分布＋dN/dS＞１となるようなサイトがあることを仮定する。NSsites=8, model=0, free parameterは４となる。この場合もd.f.は2となる。ベータ分布beta(
p,
q)に沿った分布になる。
各モデルのパラメータは以下の通りになる(Yang (2006) p.274 Table 8.2より)。
| Model |
#P |
Parameters |
| M1a (neutral) |
2 |
P0 (P1 = 1- P0) |
|
|
W0 < 1, W1 = 1 |
| M2a (selection) |
4 |
P0, P1 (P2 = 1 - P0
- P1) |
|
|
W0 < 1, W1 = 1, W2
> 1 |
|
|
|
| M7 ( beta) |
2 |
p, q |
| M8 (beta & W) |
4 |
P0 (P1 = 1 - P0) |
| |
|
p, q, Ws > 1 |
ここで、#Pはパラメータ数。モデル間の比較では、パラメータ数の引き算がdegree of freedomになる。
P0, P1, P2などが、それぞれオメガ(
W)のクラスを割り当てるサイトの割合になる。M7,8の
p,
qは、ベータ分布のパラメータになる。例えば、beta(0.2, 1)の分布は下の図の黒線で示した分布になる。PAMLでsite testを行なっているKusumi et al., (2015) Gene Genet. System.の推定では、
p=0.97924,
q = 2.89520となっていた。これをプロットすると、下の赤線のようになる。Kusumi et al. (2015)では、このような分布に従う割合が
P0=0.97924、
Ws=11.31830の値をとる
positive siteの割合が
P1=0.02076と推定されている。
Branch-site model
Branch-Site modelでは、dN/dSの値がサイト間、枝間の両方でばらつくという仮定をとる。目的は、特定の枝について、正の自然淘汰を受けている少数のサイトを検出することであり、このときの特定の枝をforeground branchesという。Yang et al., (2005), Zhang et al., (2005)に示されたテストの仕方が推奨されている。ここではnull modelとしてMA1、alternative modelとしてMAを使う。
Model A1:
帰無仮説となるMA1はmodel=2, NSsites=2, fix_omega=1, omega=1を指定する。
Model A:
対立仮説となるMAは、model=2, NSsites=2, fix_omega=0を指定する。
branch-siteモデルのパラメータをまとめると以下の通りになる（Yang (2006) p.280 Table 8.4より）
| Site
class |
Proportion |
Background W |
Foreground W |
| 0 |
P0 |
0 < W0 <
1 |
0 < W0 <
1 |
| 1 |
P1 |
W1 = 1 |
W1 = 1 |
| 2a |
(1 - P0 - P1)P0/(P0 + P1) |
0 < W0 <
1 |
W2 > 1 |
| 2b |
(1 - P0 - P1)P1/(P0 + P1) |
W1 = 1 |
W2 > 1 |
ここで、positive selectionを受けていると仮定するのが、Foregroundで指定される枝である。ファイルでは、newickファイルに#1などのマークで指定してやる必要がある。
branch-site testでは、やはり二つのモデル間のパラメータ数の差をとって、d.f.は１となる。
プログラムのチュートリアルには、LRTのやり方が解説されているが、p値を計算する際にカイ二乗（df=1)とpoint mass 0 の50:50が帰無分布となるが、それの代わりにカイ二乗（df=1）によって計算する方法を勧めている。
二つのモデルで解析が終わったら、結果ファイルにlog likelihoodが表示される。これらから、Likelihood ratio testのstatisticを計算する。すなわち、
2Δl = 2 x (l_MA - l_MA1)
ここで、l_MAはmodel Aのlog likelihood、l_MA1はmodel A1のlog likelihoodとなる。 この統計量2Δl を用いてdf=1でカイ二乗検定してp値を求めることができる。Rを使う場合は、
1 - pchisq(2Δl, df=1)
でp値を計算することができる。 例えば、2Δl＝19.88の場合、1 - pchisq(19.88, 1) = 8.245848E-06となる。
References
1. Yang, Z. (2007) PAML 4: aprogram package for phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586-1591.
2. Yang, Z., Wong, W.S.W., and Nielsen, R. (2005) Bayes empirical Bayes inference of amino acid sites under positive selection. Mol. Biol. Evol. 22, 1107-1118.
3. Zhang, J., Nielsen, R., and Yang, Z. (2005) Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level. Mol. Biol. Evol. 22, 2472-2479.
4. Yang, Z. (2006) Computational Molecular Evolution. Oxford Univ. Press
5. Kusumi, J., Tsumura, Y., and Tachida, H. (2015) Evolutionary rate variation in two conifer species, Taxodium distichum (L.) Rich. var. distichum (baldcypress) and Cryptomeria japonica (Thunb. ex L.F.) D. Don (Sugi, Japanese cedar). Gene Genet. Syst. 90, 305-315
