教程 | 柚子的祖先啥时候出来的？物种分歧推断
介绍
还是接上之前Orthofinder分析后，很多老师和同学都会选择做时间分歧分析。
Orthofinder下篇
Orthofinder上篇
"时间分歧" 简单来说，就是确定物种分化的大概时间区域。
分歧时间区域计算原理：基于已知的时间节点（已被报道的分支时间节点/具有化石时间）和分子钟计算，推算系统发育树的其他节点发生的时间。
简单补充分子钟原理：任意两个物种，自从分化成两个物种后，他们之间的遗传差异的进化速率应该与分化的时间保持相对稳定。举个例子：人与青蛙的分化点肯定比人与猴子的早，那么人与猴子的序列差异理论上会比人与青蛙的序列差异要小。
目前，采用估算物种分歧时间的程序常见有：mcmctree、r8s..
这次推送主要分享mcmctree的使用方法~
运行
PAML安装
mcmctree主要是在PAML软件包。
conda安装~
两个方法下载的paml的软件包都是4.9版本的~
主要用到paml软件包的mcmc，其中学习手册的介绍：Bayesian estimation of species divergence times incorporating uncertainties in fossil calibrations (mcmctree).
paml的学习手册：http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf
准备文件
我们在上次推文已经完成了orthofinder的分析，所以从其分析结果获mcmctree运行所需的文件。
主要需要两个文件：
1.单拷贝直系同源基因数据集合的基因ID
2. 物种的系统发育树文件
1. 准备序列比对文件
2. 准备树文件
在进行分歧时间预测，我们至少知道一个节点的大概分歧时间
TIMETREE数据库：http://timetree.org/ 是一个记录多个植物分歧进化时间的公共数据库。而且提供相关文献，我们可以进一步进行确定
这边我们暂定葡萄与拟南芥/柑橘的分歧点约为大于 100
参考orthofinder的系统进化关系
稍微解释一下如何编写这个物种树，借鉴Orthofinder的构建好的树结构，补充了葡萄分歧时间大于100MYB。
也可以直接使用 'L(1.0)'、'>1.0'
如果是在某个时间段可以：'B(1.0,1.17)'、'>1.0<1.17'
如果是小于某个时间节点可以：'U(1.17)'、'<1.17'
而括号后的两个参数分为 p=0.1, c=0.2, 由于物种数量比较小而且距离还相对较远，增加两个内置函数更严格限制已知的分歧点。
更详细请参考PAML User Guide的 第9节mcmctree的Fossil calibration 部分。
mcmctree运行
简单说一下，最主要填写的项：
seqfile = all.singlecopy.phy
treefile = treefile
outfile = out *输出前缀
ndata = 30
建议将随机采取的比对序列情况设置大一点；因为我这里六个物种单拷贝基因数量太多了，运行起来耗时多，所以我仅采用30个情况进行时间分歧计算。
usedata
设置是否利用多序列比对的数据：0，表示不使用多序列比对数据，则不会进行likelihood计算；1，表示使用多序列比对数据使用Felsenstein（1981）的似然计算，也是MCMC常使用的参数，不过计算速度较慢; 2，进行approximation likelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的in.BV文件。该文件是使用usedata = 3参数生成的out.BV文件重命名而来的。
clock = 3
model = 4
0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 说明书列举了四种情况，其中JC69是适合比较简单的关系，HKY85是相对比较远源的关系（更复杂的模型，配合alpha=0.5使用）。
alpha = 0.5
首先设置usedate=3
获得 out.BV文件
usedate=2
查看FigTree.tre
在Figtree可以可视化，后期导出SVP。
至于如何更美化这个分歧树... 就各自修行。
mcmctree.ctl配置
