根据mcmctree参数解读,mcmctree中需要调整的参数usedata
、clock
、RootAge
、rgene_gamma
、sigma2_gamma
、burnin
、sampfreq
、nsample
真正需要计算使用的rgene_gamma
、sigma2_gamma
:
1. Rough estimation of the substitution rate
树文件:请注意,根节点需要一个简单的化石校准,比如“@4.5”(450Ma ago)。如果没有这个标记计算出来的结果就没有Substitution rate is per time unit
,也就没办法计算后续。
序列文件:fas 或者 phy
使用baseml/codeml获得Substitution rate is per time unit
,通过它可以计算rgene_gamma
2. Estimation of the gradient and Hessian
在mcmctree.ctl
修改设置usedata = 3
,这告诉MCMCTREE,我们还不会进行分歧时间估计。我们只想计算gradient (g)
、Hessian (H)
、the maximum likelihood estimates (MLE) of the branch lengths
. 其生成的结果是out.BV
,将 out.BV
改名为in.BV
,作为下一步的输入文件。
3.Estimation of divergence times with the approximate likelihood method
- 设置
usedata = 2
,使用上步生成in.BV
-
clock
使用2/3都可以,区别不大,具体描述参考mcmctree参数解读。推测种群的分歧时间估计可能clock=1
比较好,分化时间较近;如果系统发育分析,涉及类群和时间尺度较大的,用clock=2/3
也许可能更好。 -
RootAge
设置一个最古老祖先的时间,这一步中不设置会报错,昆虫纲一般可以无脑<4.0
(400 Mya) -
rgene_gamma
&sigma2_gamma
设置内部参数rgene_gamma
(根据第一步的结果)、sigma2_gamma=1 4.5
(昆虫纲经验数值,可以根据情况再调整。REF:mcmctree参数解读)
burnin
、sampfreq
、nsample
遗弃代数、抽样频率、最后结果代数
一般遗弃代数设置为总运行代数的20%,抽样频率我一般选用默认的5,觉得有随机抽样频率就行,感觉理论上再大一点好。最后运行总代数,手册说的是至少2k,实际运行速度还是挺快的,可以设置几万代跑跑看,一般结果都跟现有文献差距不大。
关于rgene_gamma
、sigma2_gamma
,其实我也想过仔细研究一下,能不能精准计算一下?感觉基于Step-by-Step
手册中直接使用时间单位替换速率作为计算有点粗糙,但是想要精准计算首先需要了解一些伽马的基础知识。
如果想了解伽马分布,请看下图
https://www.zhihu.com/question/34866983
看完上面的公式,再看手册中使用Substitution rate
=m(mean of the gamma distribution)
=s(the standard deviation of the gamma distribution)
然后算算出rgene_gamma
,再根据文献选择昆虫纲的经验值sigma2_gamma = 1 4.5
,是不是合理多了?毕竟是放在PAML官网首页的推荐手册,我选择相信A Step-by-Step Tutorial
- A Step-by-Step Tutorial: Divergence Time Estimation with Approximate Likelihood Calculation Using MCMCTREE in PAML <9 June 2011>
- User Guide PAML: Phylogenetic Analysis by Maximum Likelihood <Version 4.9j (February 2020)>