mcmctree计算分歧时间步骤2023-10-02

2023-10-01  本文已影响0人  土雕艺术家

根据mcmctree参数解读,mcmctree中需要调整的参数usedataclockRootAgergene_gammasigma2_gammaburninsampfreqnsample
真正需要计算使用的rgene_gammasigma2_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


关于rgene_gammasigma2_gamma,其实我也想过仔细研究一下,能不能精准计算一下?感觉基于Step-by-Step手册中直接使用时间单位替换速率作为计算有点粗糙,但是想要精准计算首先需要了解一些伽马的基础知识。

如果想了解伽马分布,请看下图
https://www.zhihu.com/question/34866983

image.png
image.png

看完上面的公式,再看手册中使用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


上一篇 下一篇

猜你喜欢

热点阅读