..微笑着看着杯中的花茶一片片撑开.. ..透明的花瓣里水破开的声音很轻微..
分类： 信息化
2016-07-12 16:53:24
By 斑斑 <23920620>
basemlg使用同baseml相同的控制文件baseml.ctl。由于不允许使用树搜索和假定的分子钟，因此要选择runmode = 0和clock = 0。basemlg可用的模型有JC69，F81，K80，F84和HKY85，并且对位点的比率总是假定一个gamma分布。变量ncatG，given_rho，rho，nhomo在这里无效。因为迭代过程中计算S.E.s，所以参数估计的S.E.s总是会输出，因此getSE变量也没有用。
由于basemlg需要密集计算，因此对于数据分析推荐baseml中的离散gamma模型。如果你选择使用basemlg，不需要首先运行baseml，然后在运行basemlg。这可以让baseml将初始值收集到一个名为in.basemlg的文件中，来供basemlg使用。注意，basemlg执行的只是baseml中模型的一个子集。
这里有一系列密码子替代模型。详细的讨论参见Yang和Bielawski（2000），Yang（2002）和Yang（2006：第8章）。
核酸模型通常使用Goldman和Yang（1994）模型的简化版本并制定密码子i到密码子j的替代速率为
如果两个密码子在一个以上的位置不同。
同义颠换
同义转换
非同义颠换
非同义转换
（Yang等，1998）。密码子j的均衡频率（πj）可以被认为是一个自由参数，但也可以根据密码子三个位置上的核酸频率计算得到（控制变量 CodonFreq ）。这个模型下，ω = dN/dS，也就是非同义/同义替代速率。对应这个碱基模型，需要在控制文件codeml.ctl中设定model = 0，NSites = 0。
ω比率用于度量蛋白的正选择作用。简而言之，ω值 < 1，= 1，> 1表示负的纯净选择，中性进化和正选择。然而所有位点的平均比率和所有的种系几乎从不 > 1，因为正选择不可能在漫长的时间中作用于所有的位点。因此，真正要探究的只是一些种系和一些位点所受的正选择影响。
分枝模型允许系统发育树分枝中使用多个ω比率，它用于检测个别种系承受的正选择作用（Yang 1998；Yang和Nielsen 1998）。它们用变量model指定。model = 1对应自由率模型，它假定每个分枝都有一个独立的ω比率。这个模型的参数非常丰富，不推荐使用。model = 2允许你有几个ω比率。你必须指定有多少个比率和树文件中用分枝标签来标定哪个分枝对应哪个比率。参见第四章“树文件格式”中的“枝或节点标签”。溶菌酶例子的数据在文件夹examples/lysozyme/里。查看说明文档。
表2. 位点模型的参数
|
模型 |
NSsites |
#p |
参数 |
|
M0 （one ratio） |
0 |
1 |
ω |
|
M1a （neutral） |
1 |
2 |
p0 （p1 = 1 – p0）， ?ω0 < 1，ω1 = 1 |
|
M2a （selection） |
2 |
4 |
p0， p1 （p2 = 1 – p0 – p1）， ?ω0 < 1，ω1 = 1，ω2 > 1 |
|
M3 （discrete） |
3 |
5 |
p0， p1 （p2 = 1 – p0 – p1） ?ω0，ω1，ω2 |
|
M7 （beta） |
7 |
2 |
p， q |
|
M8 （beta&ω） |
8 |
4 |
p0 （p1 = 1 – p0）， p， q，ωs > 1 |
注意 #p是ω分布中自由参数的数量。括号中的参数不是自由的也不能被计算：例如，在M1a中，p1 = 1 – p0 它就不是自由参数。在比对M1a和M2a及M7和M8的似然率测试中，df = 2。位点模型通过使用NSsites指定。
位点模型允许位点中使用多个ω比率（蛋白中就是密码子和氨基酸）（Nielsen和Yang 1998；Yang等，2000b）。一些这样的模型通过设定codeml中的变量NSsites来运行（并且model = 0）。你可以通过为NSsites指定几个值来一次运行几个NSsites模型。例如，NSsites = 0 1 2 7 8在一次运行中对相同的数据计算5个模型。位点模型可用于真实数据的分析和计算机模拟研究的评估（Anisimova等，2001；Anisimova等，2002；Anisimova等，2003；Wong等，2004）。两对模型在进行正选择的两个似然率测试时似乎非常有用。第一对比对M1a（接近中性）和M2a（正选择），第二对比对M7（beta）和M8（beta和ω）。M1a（接近中性）和M2a（正选择）对模型M1（中性）和M2（选择）（Nielsen和Yang 1998）稍作了修改。两组测试中，要使用df = 2。M1a-M2a的比对似乎比M7-M8的比对更稳定。参见下面的表格。旧模型M1和M2设定ω0 = 1和ω1 = 1，而且不能对0 < ω < 1之间的位点进行计算。Wong等人（2004）和Yang等人（2005）的工作中描述的新模型M1a和M2a在设定ω1=1时可以估算出数据中的0 < ω0 < 1。版本3.14以后用于识别正选择位点的BEB程序也可以执行了。
第三组测试比对零假设的M8a（beta和ωs = 1）和M8（Swanson等，2003; Wong等，2004）。M8a用NSsite = 8，fix_omega = 1，omega = 1指定。零分布是50：50的由0和x1^2的混合体（self和Liang 1987）。5%的临界值为2.71，1%的是5.41（相对x1^2，5%的为3.84，1%的为6.63）。注意，50：50的xj^2和xk^2的混合体的p值只是两个分布的p值的平均数。在M8a-M8的比对中，我们从x1^2得到p值并且用它的一半来获得混合分布的p值。你也可以使用x1^2（Wong等，2004）。
我们建议M0-M3比对作为对位点中变量ω的测试而不是用于正选择测试。然而，对所有位点只用一个ω的模型在每一个功能蛋白中可能都会有错误，所以没有测验的点。
朴素经验贝叶斯（NEB）（Nielsen和Yang 1998；Yang等，2000b）与贝叶斯经验贝叶斯（BEB）（Yang等，2005）可用来计算位点类型的后验概率，如果似然率测试显著也可以用于识别正选择位点。NEB使用MLEs参数（如比例和ω速率）但不能计算它们的抽样误差，而BEB可以用先验贝叶斯处理抽样误差。BEB只能在M2a和M8下运行。我们建议你忽略NEB的输出，仅使用BEB的结果。
BEB输出结果格式如下：
Prob（w>1） mean w
135 K 0.983* 4.615 +- 1.329
解释：4.615是w后验分布的近似均值，1.329是w后验分布的方差的平方根。如果后验概率大于95%程序将输出*，如果大于99%则输出**。
Suzuki和Gojobori （1999）检测正选择位点的方法。使用的术语中，SG99方法检测每个位点的ω比率是否大于1或小于1。使用简约发重构祖先序列，那么对于每个位点，计算同义和非同义的差异数目以及同义和非同义的位点数。并测试位点上的ω与1有显著不同。祖先序列重构中的错误被忽略。Suzuki写了一个叫做AdaptSite的程序来执行这个测试。这只是一个直观的测试，缺乏说服力。
Codeml中，这种测试仅作为codeml或baseml中使用ML来重构祖先序列时的副产品来执行。使用RateAncestor = 1。与codeml相对的baseml的选择以及每个程序中替代模型的选择仅影响祖先序列的重构。接下来的步骤也一样，参照Suzuki和Gojobri（1999）。Codeml中我们使用M0（NSsites = 0和model = 0）。如果你愿意，你可以尝试其它模型，如NSsites = 2或8。典型的祖先重构模型几乎没什么差异。对于baseml，你需要在序列数据文件的第一行使用“GC”来指明序列是编码蛋白的。在控制文件中使用icode（= 0 用于普通密码子；= 1 用于脊椎动物线粒体密码子）来指定遗传密码子，如在codeml中。“多基因”模型与M0很接近：model = 4 Mgene = 4（参见Yang 1996a和“用于分段模型的结合分析”一节）
枝-位点模型允许蛋白中的位点和树中的枝的ω都可以变化并检测个别种系（称为前景枝）中影响很少位点的正选择。最初Yang和Nielsen执行模型A（model = 2 NSsites = 2）和模型B（model = 2 NSsites = 3）。模拟中（Zhang 2004）的测试并不顺利，因此构建了两个测试对模型A做了改动（表3）（Yang等，2005；Zhang等，2005）。测试2也是正选择的枝-位点测试，推荐使用这个测试。这个比较是将修改后的模型A与相对应的ω2设为1的空模型（fix_omega =1，omega = 1）对照的。零分布是由0和x1^2按50：50混合而成的，5%的临界值为2.71，%1的临界值为5.41。我们推荐替换x1^2（用临界值3.84和5.99）来防止模型假设的违例。
类似的都要对修改后的枝-位点模型A（不是模型B）执行针对位点类型计算后验概率的NEB和BEB方法。你需要结合BEB来使用模型A并忽略NEB的输出。
表3. 枝-位点模型的参数
|
位点类型 |
旧模型A （np = 3） |
新模型A （np = 4） |
|||
|
|
比例 |
背景 |
前景 |
背景 |
前景 |
|
0 |
p0 |
?ω0 = 0 |
?ω0 = 0 |
0 <ω0 < 1 |
0 <ω0 < 1 |
|
1 |
p1 |
?ω1 = 1 |
?ω1 = 1 |
?ω1 = 1 |
?ω1 = 1 |
|
2a |
（1 – p0 – p1） p0/（p0 + p1） |
?ω0 = 0 |
?ω2 >=1 |
0 <ω0 < 1 |
?ω2 >=1 |
|
2b |
（1 – p0 – p1） p1/（p0 + p1） |
?ω1 = 1 |
?ω2 >=1 |
?ω1 = 1 |
?ω2 >=1 |
注意 枝位点模型A用model = 2 Nssites = 2指定。枝位点测试的正选择是个可选的模型。空模型设定ω2 = 1。似然率测试 df = 1（见文本）。
进化枝模型通过model = 3 Nssites = 2指定，而进化枝模型D通过model = 3 Nssites = 3指定并用ncatG来指定位点类型的数量（Bielawski和Yang 2004；也可以参见Forsberg 和Christiansen， 2003）进化枝模型C做了些类似枝位点模型A所做的一些变化。新的模型C用 0 <ω0 < 1替换了ω0 = 0，并在ω分布中设置了5个参数：p0，p1，ω0，ω2和ω3。新的模型C可以同新的M1a（两个参数，d.f. ≈ 3）（接近中性）比较。
表4. 进化枝模型中的参数 （新旧对比）
|
|
旧模型C （np = 4） |
新模型C （np = 5） |
|||
|
位点类型 |
比例 |
进化枝1 |
进化枝2 |
进化枝1 |
进化枝2 |
|
0 |
p0 |
?ω0 = 0 |
?ω0 = 0 |
0 <ω0 < 1 |
0 <ω0 < 1 |
|
1 |
p1 |
?ω1 = 1 |
?ω1 = 1 |
?ω1 = 1 |
?ω1 = 1 |
|
2 |
p2 = 1 – p0 – p1 |
?ω2 |
?ω3 |
?ω2 |
?ω3 |
进化枝模型D可以使用ncatG = 3 或 2。当指定枝位点模型A或B及进化枝模型C，并且模型中的类别数量确定时选项ncatG的变化是被忽略的。
BEB程序可以对模型C执行但不能用于模型D。你需要使用模型C与BEB程序结合。忽略NEB的输出。
模型C和模型D做了一个扩展允许针对两个以上的进化枝或分支类型。分支类型通过树文件中的标签指定。如果你有四个分支类型（标签#0，#1，#2，#3），进化枝模型如下。这里的ω2，ω3，ω4，ω5是（0，∞）范围中最佳的独立参数。
进化枝模型C下的BEB计算开销非常大（每个额外的分枝类型要增加10翻的计算量），因此只能用于分枝类型少的模型。
表5. 含两个以上分枝类型的进化枝模型C或D的参数
|
位点类型 |
比例 |
进化枝1 |
进化枝2 |
进化枝3 |
进化枝4 |
|
0 |
p0 |
0 <ω0 < 1 |
0 <ω0 < 1 |
0 <ω0 < 1 |
0 <ω0 < 1 |
|
1 |
p1 |
?ω1 |
?ω1 |
?ω1 |
?ω1 |
|
2 |
p2 = 1 – p0 – p1 |
?ω2 |
?ω3 |
?ω4 |
?ω5 |
进化枝模型C中，ω1是固定的，而在进化枝模型D中，ω1是作为自由参数估算的。
突变-选择模型（Yang和Nielsen，2008）使用下面的控制变量执行。
CodonFreq = 7 * 0：1/61 each， 1：F1X4， 2：F3X4， 3：codon table
* 4：F1x4MG， 5：F3x4MG， 6：FMutSel0， 7：FMutSel
estFreq = 1
CodonFreq = 6 指定FMutSel0，7指定FmutSel。如果estFreq = 1，frequency/fitness参数通过ML从数据中估算，如果estFreq = 0，则用数据中的观察频率计算。模型的详细内容参见Yang和Nielsen （2008）。examples/mtCDNAape/文件夹中的说明文档说明了如何重复文献中表1的结果。
我对经验和机械的氨基酸替代模型做了区分（Yang 等，1998； 2006：第2章）。Codeml中的经验模型包括Dayhoff （Dayhoff 等1978），JTT （Jones 等1992）， WAG （Whelan和Goldman 2001）， mtMam （Yang 等1998）， mtREV （Adachi和Hasegawa 1996a）等。它们从真实数据中按普通时间可逆模型估算替代率参数。
机理模型是考虑到氨基酸间的突变距离构成的，突变距离通过编码密码子的位置决定，对自然选择的突变的过滤是在蛋白水平上操作的（Yang等 1998）。Aaml程序实现了很少的这类模型，通过变量aaDist指定。
由于基于密码子的分析和基于氨基酸的分析使用不同的模型，一些控制变量也就有了不同的意义，针对密码子和氨基酸序列使用不同的控制文件是个不错的主意。codeml默认的控制文件是codeml.ctl，如下：
seqfile = stewart.aa * sequence data file name
outfile = mlc * main result file name
treefile = stewart.trees * tree structure file name
noisy = 9 * 0，1，2，3，9： how much rubbish on the screen
verbose = 0 * 1： detailed output， 0： concise output
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：TipDate
aaDist = 0 * 0：equal， +：geometric; -：linear， 1-6：G1974，Miyata，c，p，v，a
* 7：AAClasses
aaRatefile = wag.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， 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-11：see below
Mgene = 0 * 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 = 3 * # of categories in dG of NSsites models
fix_rho = 1 * 0： estimate rho; 1： fix it at rho
rho = 0. * initial or fixed rho， 0：no correlation
getSE = 0 * 0： don't want them， 1： want S.E.s of estimates
RateAncestor = 0 * （0，1，2）： rates （alpha>0） or ancestral states （1 or 2）
Small_Diff = .5e-6
* cleandata = 0 * remove sites with ambiguity data （1：yes， 0：no）?
* fix_blength = 0 * 0： ignore， -1： random， 1： initial， 2： fixed
method = 0 * 0： simultaneous; 1： one branch at a time
变量seqfile， outfile， treefile， noisy， Mgene， fix_alpha， alpha， Malpha， fix_rho， rho， clock， getSE， RateAncestor， Small_Diff， cleandata， ndata， fix_blength和method的使用与前面baseml.ctl中描述的一样。Seqtype在数据中指定序列的类型；seqtype = 1表示密码子序列（程序为codeml）；2表示氨基酸序列（程序为aaml）；3表示被翻译成蛋白的密码子序列。
CodonFreq 在密码子替代模型中指定均衡的密码子频率。这些频率可以假设为相等的（每个标准遗传密码都是1/61，CodonFreq = 0），或从平均核酸频率计算得到（CodonFreq = 1），或来自三个密码子位置的平均核酸频率（CodonFreq = 2），又或使用自由参数（CodonFreq = 3）。这些密码子频率模型中参数的数值为0，3，9和60（通用密码子），分别对应CodonFreq = 0，1，2和3。
Clock 参见baseml控制文件。
aaDist 指定是否假定氨基酸距离是相等的（= 0）或使用Grantham矩阵（= 1）（Yang等，1998）。论文中用于分析的线粒体数据示例在软件包中的example/mtdna文件夹。aaDist = 7（AAClasses），可以使用密码子和氨基酸序列来运行，允许你有一些氨基酸替代类型并让程序估算它们的不同速率。Yang等，（1998）使用了这个模型。替代类型的数量和氨基酸对的变化属于那种类型是在叫做OmegaAA.dat的文件中指定的。你可以对“保守”和“变化”的氨基酸替代情况使用模型来配合不同的ω速率。examples/mtCDNA中的例子是对应这个模型的。查看文件夹内的说明文件。
Runmode = -2 执行ML来估算成对比较的蛋白编码序列中的ds 和dn（seqtype = 1）。程序将估算出的ds和dn结果输出到2ML.ds和2ML.dn。由于许多使用者似乎对种系中的dn/ds比率很感兴趣，那么他们可能也对通过计算两个比率得到的枝长来表示的树型的检验感兴趣，尽管这个分析有点特别。如果你的物种名不超过10个字符，你可以使用输出的距离矩阵作为Phylip程序的输入，这不需要任何改动。否则你需要编辑文件将名字裁短。对于氨基酸序列（seqtype = 2），选项runmode = -2使程序在替代模型下通过数值迭代来计算ML距离，要么在所有位点使用一个比率的模型下（alpha = 0），要么在位点的比率为gamma模型下（alpha ≠ 0）。在后者中，使用连续的gamma并且变量ncatG被忽略。（runmode = 0，使用离散gamma。）
Model指定枝模型（Yang 1998；Yang 和 Nielsen 1998）。Model = 0表示对所有分枝使用一个比率；1表示每个分枝一个比率（自由比率模型）；2表示任意数目的比率（如2-ratios或3-ratios模型）。当model = 2时，你必须将树中的分枝用符号#或$进行分组。参见前面关于指定枝位点模型的章节。Model = 2，变量fix_omega设定最后的ω比率（ωk-1如果总计有k个比率）为文件中指定的omega的值。这个选项用于测试作为前景枝的比率是否显著于其它的值。参见examples/lysozyme文件夹来重复Yang（1998）中的结果。
NSsites 指定位点模型，NSsites = m相当于Yang（2000b）中的Mm模型。变量ncatG用于指定一些模型下ω分布中的类别数量。在Yang等（2000b）的工作中，3对应M3（离散），5对应M4（频率），10对应连续分布（M5：gamma，M6：2gamma，M7：beta，M8：beta&w，M9：beta&gamma，M10： beta&gamma+1， M11：beta&normal>1， and M12：0&2normal>1， M13：3normal > 0）。例如M8有11个位点类别（10个来自beta分布加上1个额外的用于正选择的类别，同时ωs≥1）。参见前面密码子替代模型章节中关于Paml3.14中M1和M2变化的部分。
你可以在一个分枝中运行几个NSsites模型。如下，在Yang等（2000b）中，使用默认的ncatG并使用
NSsites = 0 1 2 3 7 8
这些位点的后验概率以及这些位点预期的ω值列在了rst文件中，它可能用于识别正选择位点。
参见examples/hivNSsites/和examples/lysine/下使用位点模型分析的例子。
枝位点模型A（参见前面的密码子替代模型章节）通过使用变量model和NSsites指定。
Model A： model = 2， NSsites = 2， fix_omega = 0
对枝位点正选择测试的模型是可选的。空模型也是枝位点模型，但是设定了ω2 = 1，并指定
Model A1： model = 2， NSsites = 2， fix_omega = 1， omega = 1
这里有些值得注意的是对与两个旧的测试我们并不推荐。旧的枝位点模型B（Yang 和 Nielsen 2002）通过下面参数指定。
Model B： model = 2， NSsites = 3
旧测试B的空模型是含2种位点类型的NSsites模型3（离散）：
model = 0， NSsites = 3， ncatG = 2
使用d.f. = 2。Yang和Nielesn（2002）的大规模的溶解酶数据集分析在examples/lysozyme文件夹中。
同样适用于Yang等（2005）中描述的枝位点测试1和Zhang等（2005）使用的修改后的枝位点模型A作为被择假设，而零假设是d.f. = 2的新位点模型M1a（接近中性）。当前景枝在宽松的选择约束或正选择压力时这个测试是显著的。建议使用测试2来替代它，测试2也称为正选择的枝位点测试。
Bielawski和Yang (2004)的进化枝模型C和D通过下面方式指定：
Model C： model = 3， NSsites = 2
Model D： model = 3， NSsites = 3 ncatG = 2
详细内容参见文献。类似的模型C也进行了修改并且BEB程序只能用于模型C。参见前面。
Icode 指定遗传密码。十一个遗传密码表通过使用icode = 0到10来应用，它们分别对应于GeneBank中的transl_table 1到11。0为通用密码子；1为哺乳动物线粒体密码子，3为霉菌线粒体；4为无脊椎线粒体；5为纤毛虫细胞核密码子；6为棘皮动物线粒体；7为euplotid线粒体；8为替代酵母细胞核；9为海鞘线粒体；10为赭纤虫细胞核。还有一种额外的密码子，叫做Yang氏规则密码子，通过icode = 11来指定。这套密码子中有16种氨基酸，密码子第一位和第二位的偏差为非同义替代，第三位的偏差为同义替代。也就是说所有的密码子都是四重简并。目前没有任何生物体使用这套密码子。
Mgene 同序列文件中的选项G结合使用，按表6中概括的来指定分割模型（Yang 和 Swanson2002）。文章中使用的溶菌酶数据包含在examples/文件夹中。分析将溶解酶中的氨基酸分为包埋的和暴露的两部分（“基因”），并在分析中对各部分的异质性做了解释。你可以阅读说明文件并重复表6中Yang和Swanson（2002）的结果。
表6. 密码子替代的分割模型的设置
|
序列文件 |
控制文件 |
跨越基因的参数 |
|
无G |
Mgene = 0 |
每个都相等 |
|
选项 G |
Mgene = 0 |
（k，ω）和π相同，但cs不同。（成比例的枝长） |
|
选项 G |
Mgene = 2 |
（k，ω）相同，但πs和cs不同。 |
|
选项 G |
Mgene = 3 |
π相同但（k，ω）和cs不同 |
|
选项 G |
Mgene = 4 |
（k，ω），πs和cs都不同 |
|
选项 G |
Mgene = 1 |
单独分析 |
Fix_alpha，alpha：对于密码子模型，fix_alpha和alpha为位点指定gamma比率模型，模型中密码子中针对位点变化的相关比率服从gamma分布，但是ω比率在全部位点中保持不变。这是对核酸和氨基酸gamma模型的一个懒惰的扩展。我不喜欢这个模型并且建议你用NSsites模型来替代它（使用NSsites变量来指定，并且设定fix_alpha = 1，alpha = 0）。如果你对这个程序同时指定NSsites和alpha，它将中止运行。
RateAncestor：参见baseml控制文件中的描述。
密码子序列（seqtype = 1）的输出：每个序列中的密码子频率被计算并列在遗传密码表中，以及它们在各个物种中的总和。每个表格包含六个或更少的物种。对于多基因数据（序列文件中指定了选项G），每个基因中的密码子频率也会被列出。密码子位置上的核酸分布也会列出。Nei和Gojobori (1986)的方法用来计算每个同义替代位点（ds）上发生了同义替代的数量以及每个非同义替代位点（dn）上发生了非同义替代的数量及它们的比值（dn/ds）。这些用于构建对似然性分析中枝长的最初估算，但是它们不能似然性估算自身。
祖先重构（RateAncestor）的结果在rst文件中。在计算位点中变量dn/ds比率的模型里（NSsites模型），位点类型的后验概率和正选择位点也都在rst文件中。
Model指定氨基酸替代模型：0为泊松模型，它假设所有的氨基酸替代是均等的速度（Bishop 和 Friday 1987）；1为比例模型，该模型中一个氨基酸改变的速度是与氨基酸的频率成比例的。Model = 2指定一类经验模型，先验的氨基酸替代速率矩阵通过aaRatefile文件指定。文件包含在软件包中用于Dayhoff（dayhoff.dat），JTT，WAG（wag.dat），mtMAM（mtmam.dat），mtREV24（mtREV24.dat）等的经验模型。
如果你想指定自己的替代速率矩阵，可以参看这些文件，里面注释了文件的结构。氨基酸替代模型的其它选择被忽略。简而言之，变量model，aaDist，CodonFreq，NSsites和icode用于指定密码子模型（seqtype = 1），而model，alpha和aaRatefile用于氨基酸序列。
Mgene 同序列文件中的选项G结合使用，按表6中概括的来指定分割模型（Yang 和 Swanson2002）。文章中使用的溶菌酶数据包含在examples/文件夹中。分析将溶解酶中的氨基酸分为包埋的和暴露的两部分（“基因”），并在分析中对各部分的异质性做了解释。你可以阅读说明文件并重复表6中Yang和Swanson（2002）的结果。（此处与上面完全相同，原文中如此，可能是作者编辑时出现了错误。）
表7. 氨基酸替代的分割模型的设置
|
序列文件 |
控制文件 |
跨越基因的参数 |
|
无G |
Mgene = 0 |
每个都相等 |
|
选项G |
Mgene = 0 |
π相同，但cs不同。（成比例的枝长） |
|
选项G |
Mgene = 2 |
πs和cs都不同 |
|
选项G |
Mgene = 1 |
单独分析 |
Runmode 同baseml.ctl中的工作方式一致。指定runmode = 2将使程序在成对比较中计算ML距离。你可以改变控制文件中的变量aaRatefile，model和alpha。
如果你进行成对的ML比较（runmode = -2）并且数据包含模糊字符或者对齐空位如果cleandata = 1，程序将从所有序列中移除这类的位点。这也称为“全部删除”。如果cleandata = 0，将成对的移除对齐空位和模糊字符。{{这似乎不是真的。如果runmode = -2,程序当前移除所有模糊的位点。需要检查。注意ziheng 31/08/04}}注意，在种系发生的多序列的似然性分析中，如果cleandata = 0，对齐空位被当做模糊字符对待，如果cleandata = 1，对齐空位和模糊字符都会被删除。注意移除对齐空位和将它们作为模糊字符对待都会降低序列分歧。数据中的模糊字符（cleandata = 0）使似然率计算变慢。
氨基酸序列（seqtype = 2）输出：输出文件不需解释，它与核酸和密码子分析的结果非常相似。氨基酸替代的经验模型（通过dayhoff.dat，jones.dat，wag.dat，mtmam.dato或mtREV24.dat指定）在替代速率矩阵中不包括任何参数。当RateAncestor = 1时，祖先重构的结果在rst文件中。用位点变化速率模型计算位点替代速率的结果在rates文件中。
这个程序生成一个简单的菜单，如下：
(1) Get random UNROOTED trees?
(2) Get random ROOTED trees?
(3) List all UNROOTED trees into file trees?
(4) List all ROOTED trees into file trees?
(5) Simulate nucleotide data sets (use MCbase.dat)?
(6) Simulate codon data sets (use MCcodon.dat)?
(7) Simulate amino acid data sets (use MCaa.dat)?
(8) Calculate identical bi-partitions between trees?
(9) Calculate clade support values (read 2 treefiles)?
(0) Quit?
选项 1，2，3，4。这个程序可以用于生成一个随机树，无根或有根都可以，也可以有枝长或无枝长。当然，你只能用于少量的物种否则你的硬盘将被无用的树塞满。选项8用于从一个树文件中读取多个数并且计算两部分的距离，要么是第一个和剩余的树之间，要么是每对之间。
选项9（进化枝支持值）可用于概括自举和贝叶斯分析。它读入两个树文件。第一个树文件只能包含一棵树，如最大似然树。你在文件的前端需要有物种的数量和树的数量（只能是1棵树的）。第二个树文件可以包含树的集合，比如从1000次自举假样本中估算出的1000个最大似然树。这个选项会为树文件里的ML树中的每个进化枝计算自举支持值，也就是说，第二个树文件中的比例树包含第一个文件中树的进化枝或节点。第二个树文件中的第一行并不必须含有树和物种的数量。如果你运行MrBayes，你可以将最大似然树或最大经验树移动到第一个文件中，第二个树文件可以由MrBayes生成.t文件来替代，不需要任何改动。现在物种只是用树文件中的数字来表示。你可以选择这个选项9来运行evolver。程序将向你询问两个输入文件的名字。另一种方法是绕开这个简单的菜单，而在命令行中输入选项和两个文件的名字：evolver 9
选项 5，6，7（模拟）evolver程序模拟带有用户指定树拓扑结构和枝长的核酸，密码子和氨基酸序列。用户在控制文件中指定替代模型和参数。见下面。程序在一个文件中以PAML（输出mc.paml）或PAUP（mc.paup）格式生成多个数据集如果你选择PAUP格式，程序将寻找名为paupstart（程序将其拷贝到数据文件的前面），paupblock（程序将其拷贝到每个模拟数据集的末尾）和paupend（程序将其合并在文件的末尾）文件。
这可以是PAUP程序在一次运行中分析全部数据集。模拟的参数在三个文件中指定：MCbase.dat， MCcodon.dat和MCaa.dat分别对应核酸，密码子和氨基酸序列。运行默认的参数并在屏幕输出中查看。然后查看appropriate .dat 文件。作为例子，MCbase.dat文件在下面列出。注意文件的第一块作为evolver的输入，其它的为注释。树的长度是种系发生中沿着所有分枝的每位点替代数的期望值，计算为枝长的总和。当我用模拟来评估序列分歧的效率又要保持树型的固定时，引入了这个变量。Evolver将测量树来计算指定树长的枝长。如果你对树长使用-1，程序将使用树中给定的枝长而不再重新计算。还要注意碱基频率必须按固定的顺序。MCaa.dat和MCcodon.dat中的氨基酸与密码子频率也一样。
0 * 0: paml format (mc.paml); 1:paup format (mc.nex)
367891 * random number seed (odd number)
5 1000000 1 * <# seqs> <# nucleotide sites> <# replicates>
-1 *
(((A :0.1, B :0.2) :0.12, C :0.3) :0.123, D :0.4, E :0.5) ;
3 * model: 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85, 5:T92, 6:TN93, 7:REV
5 * kappa or rate parameters in model
0 0 * <#categories for discrete gamma>
0.1 0.2 0.3 0.4 * base frequencies
T C A G
Evolver的模拟选项（5，6，7）可以使用命令行格式避开简单的菜单。因此下面是一些运行evolver可能的方式。
evolver
evolver 5 MyMCbaseFile
evolver 6 MyMCcodonFile
evolver 7 MyMCaaFile
密码子替代模型使用选项6来假设系统发生树中所有分枝和基因上的所有位点都具有相同的ω比率。这就是所谓的M0模型（一比率）。在位点中使用带有ω变量的位点模型来模拟，或者在分枝中使用带有不同ω的分枝模型，又或在位点和分枝中使用变化的ω的枝位点模型，请阅读paml/Technical/Simulation/Codon/文件夹下的CodonSimulation.txt文件。
第一个变量控制使用的序列数据格式：0表示paml/phylip格式，1位点模式记数和2nexus格式。位点模式记数可以被baseml或codeml读入，如果你有带有很多位点的大规模数据（>106），这是很有用的。（参见前面的序列数据格式部分）
Evolver也可以对每个模拟数据使用带有随机枝长的随机树进行模拟。你必须改变源代码并重新编译。打开evolver.c然后找到函数Simulate中的fixtree = 1并将它改成0。重新编译并将程序命名为evolverRandomTree，如。
cc -o evolverRandomTree -O2 evolver.c tools.c –lm evolver 5 MCbaseRandomTree.dat
控制文件如MCbase.dat也需要改变。软件包中含有一个名字为MCbaseRandomTree.dat的例子文件。“树长”和树拓扑结构的一行被出生率λ，死亡率μ，样本比率ρ和树高所替代（从根到顶端的每个位点替代数的期望值）。树用带有物种样本的生-灭过程生成（Yang和Rannala 1997）。时钟是假设的。你必须改变源码来释放时钟。如果你对文件格式选择0（paml），随机树将输出在文件ancestral.txt中；它可以用Rod Page的TreeView查看。如果文件格式选择2（nexus），程序将把树输出在序列文件的树的一块内容中。
Evolver程序还有几个选项用来列出几个少量物种的所有的树，并从分枝进化，含有物种样本的生-灭过程模型中生成随机树（Yang和Rannala 1997）。它还有一个选项来计算树拓扑结构间的分歧距离。
Evolver中使用的蒙特卡罗模拟算法。你可以在“模型与分析”一节中读到更相信的内容。参见第九章中Yang（2006）。这有一些简短的注释。Evolver通过沿着树的“进化”序列模拟数据集。首先，靠模型或数据文件指定使用均衡的核酸，氨基酸或密码子频率对根生成一个序列（MCbase.dat，MCcodon.dat和MCaa.dat）。每个位点按照枝长、替代模型的参数等沿着树枝进化。当序列中的位点按照相同的过程进化，转移概率矩阵对每个枝的所有位点只计算一次。对所谓的位点-分类模型（如gamma，NSsites密码子模型），不同位点可能拥有不同的转移概率矩阵。枝起始处的字符和结尾处的字符是从通过源字符得到的转移概率所指定的多项式分布中抽样得来的。通过模拟生成祖先节点的序列并输出在ancestral.txt文件中。
一些人想在根处指定序列而不是让程序随机生成序列。这可以通过在RootSeq.txt中放入一条序列来实现。这个序列不能包含模糊字符或空位或终止密码子。在几乎所有的模拟中，确定根序列是完全错误的，因此你需要避免犯这种错误。如果你想通过模拟来反映你的特殊基因，你可能要用一个模型从那个基因来估算参数并使用参数估计来模拟数据集。
程序yn00执行Yang和Nielsen（2000）的方法来估算两个序列间的同义和非同义替代速率（ds和dn）。Nei和Gojobori（1986）的方法也包含在内。程序中的这个专门方法用来解释颠换/转换比的偏差和密码子使用偏差，它是假设密码子频率模型为F3×4时的ML法的近似值，用来解释颠换/转换比率速率。我们推荐你尽可能使用ML法（控制文件中runmode= -2，CodonFreq = 2），即使对成对比较的序列也是。
seqfile = abglobin.nuc * sequence data file name
outfile = yn * main result file
verbose = 0 * 1: detailed output (list sequences), 0: concise output
icode = 0 * 0:universal code; 1:mammalian mt; 2-10:see below
weighting = 0 * weighting pathways between codons (0/1)?
commonf3x4 = 0 * use one set of codon freqs for all pairs (0/1)?
上面展示了一个yn00的控制文件，指定了序列文件的名称（seqfile），输出文件的名称（outfile）和遗传密码子（icode）。序列中包含对齐空位和模糊核酸的位点将从全部的序列中移除。变量weighting决定当在密码子之间计算时使用相等的权重或不等的权重。两种方法对相异的序列是不同的，不等的权重将计算的较慢。从数据文件中对所有序列估算颠换/转换比速率并用于随后的成对比较。Commonf3×4指定密码子频率是从数据文件中的每对序列间估算还是所有的序列中估算。除了主要的结果文件，程序还生成三个距离矩阵：同义替代率的2YN.ds，非同义替代率的2YN.dn，结合密码子比率的2YN.t（t是作为每密码子核酸替代率的数目来估算的）。这些都是下对角线矩阵并且可被一些诸如Felesenstein 开发的PHYLIP中的NEIGHBOR之类的程序直接读取。
Mcmctree程序可能是第一个贝叶斯系统发育程序（Yang和Rannala1997；Rannala和Yang1996），但它是运行的比较慢，自从MrBayes（Huelsenbeck和Ronquist 2001）发布后就退役了。
从PAML3.15（2005）版本起，mcmctree执行Yang和Rannala（2006）的MCMC算法以及Rannala和Yang（2007）的算法使用多重化石校正的有根树来估算物种分歧时间。这与Jeff Thorne和Hiro Kishino的multidivtime程序类似。Yang和Rannala（2006）和Yang（2006，7.4）中讨论了两个程序间的不同。也可以参见后面。
请参考任何与贝叶斯计算相关的书籍，如，Yang（2006）MCMC算法基础中的第五章。
这个程序还有些需要注意的地方：
在开始运行程序前将窗口从80列设置为100列。（在windows XP/Vista中，右键命令提示符装口标题栏并改变属性-布局-窗口大小-宽度）
树文件中的树必须是有根的二叉树：每个内部节点必须只有两个子节点。你不能用MCMCTREE来计算一个多分枝的树的分歧时间。相对的你需要使用分支的ML树或NJ树或传统的形态学的树。注意，一个二叉树有一个调整的机会，而多叉树则没有。
树不能有枝长。如，((a:0.1, b:0.2):0.12, c:0.3) '>0.8<1.0'; does not work, while ((a, b), c) '>0.8<1.0'就很好。
在不严格的时钟模型下（clock = 2 或 3）并且根部没有校准，宽松的上限值（最大年代的约束）必须通过控制文件指定（RootAge）。（如果clock = 1就没有必要使用RootAge，如果程序坚持要这个参数，你必须指定。我会尝试来订正这个问题。）
Time unit需要选择满足节点年代在0.01 – 10之间的时间单位。如果分歧时间在100 – 1000MY，那么100MY就是一个时间单位。基于你选择的时间尺度的时间的先验值和速率以及化石校准都要指定好。例如，如果一个时间单位是10MY，那么
rgene_gamma = 2 20 * gamma prior for overall rates for genes
sigma2_gamma = 10 100 * gamma prior for sigma^2 (for clock=2 or 3)
就表示每10MY每位点的总平均替代速率是2/20 = 0.1或者说每年每位点的替代率为10-8，这对哺乳类线粒体基因是合理的。这里的gamma先验值形状参数为α= 2，相当的弥散。如果你改变时间单位，你需要保持形状参数的固定，并改变尺度参数β为一个恰当的平均值。换句话说，使用100MY做为时间单位，上面的就变成了
rgene_gamma = 2 2 * gamma prior for overall rates for genes
sigma2_gamma = 10 100 * gamma prior for sigma^2 (for clock=2 or 3)
不要改变预设的σ2，因为σ2是log比率的方差：当你用一个常数来调节速度时，这个对数比值的方差不能改变。树文件中的化石校准需要相应的改变。理想的状况就是当改变时间尺度时，生物学的结果没有改变。我们知道，模型的两个部分（速度的log正态分布和时间的生-灭模型）在时间的尺度上不是不变的。虽然如此，Mathieu Groussin的测试表明时间尺度的选择对推导出的时间和估算的速率影响是非常小的。
你要运行同一个分析至少两次来证实不同运行间的结果是非常相似的（尽管不是完全一致），这一点很重要。关键是你需要确认接受的比例既不是太高也不是太低。关于变量的微调见后面。
运行没有序列数据的程序是首先要检查预先指定的分歧时间的举止和CIs。理论上，全部时间的联合先验分布需要用户自己指定。尽管如此，指定这样一个复杂的高维分布几乎是不可能的。取而代之的是程序使用校正分布和根部约束生成的联合先验值和生-灭过程生成的联合先验值。程序使用这个先验值进行年代测定分析。你必须检查它确保它是明确无误的，用你对物种的知识及有关的化石记录来判断。如果需要，你可能必须改变化石标尺以便先验值看起来更合理。
程序现在对MCMC样本做一个简单的概括，计算均值，中值和对参数95%的CIs。如果你想更复杂的总结，如1-D和2-D密度测试，你可以在mcmctree运行末尾通过键入ds mcmc.out运行一个小程序ds。
你可以使用强制的界限来编译程序，设定末尾的概率为10-300来取代0.025。从mcmctree.c文件的起始处复制带有-DHardBounds标志的命令到命令行来编译一个名为mcmctreeHB的可执行程序。
化石校准信息以分歧时间（或物种树的节点年代）的统计分布的格式来通过树文件指定。表8中进行了概括。这里的化石表示包括地址事件等任何形式的外部校准数据。对于一个明确的分析，树中必须包含一个下界和一个上界，即使他们可能不在同一个节点。Gamma，斜正态，斜t分布都可以充当两个边界。因此这样的一个标度足够标定这棵树来确保一个明确的分析。
表8 校准分布
|
校准 |
#p |
规格 |
密度 |
|
L（下游或最低边界） |
4 |
'>0.06' 或'L(0.06)' 或 'L(0.06, 0.2)' 或 'L(0.06, 0.1, 0.5)' 或 'L(0.06, 0.1, 0.5, 0.025)' |
L(tL, p, c, pL) 用偏移量p，规模参数c和左侧的尾部概率pL指定最低年代边界tL。默认值为 p = 0.1, c = 1，和pL = 0.025, so >0.06 或 L(0.06) 表示L(0.06, 0.1, 1, 0.025)，和L(0.06, 0.2)表示L(0.06, 0.2, 1, 0.025)。如果你想让最低边界可靠，用pL = 1e-300，但不要使用pL = 0。 |
|
U （上游或最高边界） |
2 |
'<0.08' 或 'U(0.08)' 或 'U(0.08, 0.025)' |
YR06中的Eq. 16和fig. 2b U(tU, pR)用右尾部概率pR指定最高年代边界tU。默认值为pU = 0.025, so <0.08 或 U(0.08) 表示U(0.08, 0.025). 例如U(0.08, 0.1)表示8MY这个最高边界有10%的可能性被违反（例如真实的年代早于8MY） |
|
B（下游和上游或最低和最高边界） |
4 |
'>0.06<0.08' 或 'B(0.06, 0.08)' 或 'B(0.06, 0.08, 0.025, 0.025)' |
YR06中的Eq. 17 & fig. 2c。 B(tL, tU, pL, pU)用左尾概率pL和右尾概率pU 指定了一对儿边界，以便真实的年代在tL和tU之间。默认值为pL = pU = 0.025. |
|
G （gamma） |
2 |
'G(alpha, beta)' |
YR06中的Eq. 18 & fig. 2d。 |
|
SN （斜正态） |
3 |
'SN(location, scale, shape)' |
Eq. 2和下面的图。 |
|
ST（斜t） |
4 |
'ST(location, scale, shape, df)' |
Eq. 4和下面的图。 |
|
S2N（斜2正态） |
7 |
'SN2(p1, loc1, scale1, shape1, loc2, scale2, shape2)' |
p1: 1 – p1 两个斜正态的混合。 |
注意 #p为用户提供的分布中的参数值。YR06（Yang和Rannala2006）中的图2就是Yang（2006）中的图7.11。
（1）下界（最小年代）通过‘>0.06’或者‘L（0.06）’指定，表示节点年代至少在6百万年。这里我们假定单位是1亿年。Paml4.2中对最小年代的计算有所改变。一种基于柯西分布的重尾密度分布（Inoue等2010）替代了Yang和Rannala（2006）中图2或Yang（2006）图7.11a中的错误的扁平密度（flat density）分布。带有位置参数tL（1+p）和范围参数tLC的柯西分布缩短了tL，并在tL的左侧增加了2.5%的密度使其具有柔性。结果的分布为tL（1+p），2.5%到97.5%限制在tL和Tl[1+p+c tan(π(1/2-0.025A/0.975))]之间。如果最小边界tL基于较好的化石数据，真实的分析时间可能接近最小边界，以便使用一个小的p值和c值。值得注意的是c比p对后面的时间估算有更强的影响。程序默认的值为p = 0.1 和 c = 1。然而建议你通过对化石数据仔细评估后，对不同的下界使用不同的p值和c值。下面是这种分布的一些平面图。下界设定为tL = 1，但是单位可能是任意的，如1亿年或者十亿年。对p的每个值（0.1和0.5），四条曲线对应c = 0.2，0.5，1，2（从最上面的峰到最下面的峰）。2.5%界限为1，而当p = 0.1时，97.5%界限各自对c的值为4.93，12.12，24，43，49.20，当p = 0.5时，分别为4.32，9.77，20.65，44.43（注意，这些值在Inoue等2010中计算的不正确）。
（2）上界（最大年代）通过如‘0.08’或‘U（0.08）’指定，表示节点年代最大为8百万年。
（3）上下界通过‘>0.06<0.08’或‘B(0.06,0.08)’对同一个节点指定，表示这个节点的年代在6百万至8百万年之间。注意，在所有的三种标定（L,U,B）中，边界都是柔性的，年代有2.5%的概率偏离边界（参见Yang和Rannala2006中图2；或Yang2006中的图7.11）
（4）伽马分布‘G（188，2690）’指定了形状参数 α = 188，范围参数 β = 2690的伽马分布。
