最近在做硕士毕业课题的收尾工作，正好做到了物种分歧时间的部分，虽然之前也做过，但是还是复习一遍。
1.安装Paml
conda install -c biodonda paml
2.配置带有时间校准点的系统发育树
首先说明Mcmctree中时间的单位是100Mya为一个单位，例如60Mya就是0.6，160Mya就是1.6。配置时间标准店的方法为在你想要添加的节点位置的后方添加时间范围，例如‘>.3302<.6386’，代表的就是33.02Mya-63.86Mya，这里的范围一般指的是95％置信区间。以下是一个时间校准点的系统发育树的例子。
76 1
(Asparagus_officinalis,((((Liriope_muscari,Liriope_spicata),(Ophiopogon_chingii,(Ophiopogon_jaburan,Ophiopogon_japonicus))),Peliosanthes_macrostegia),(((((Disporopsis_aspersa,Disporopsis_pernyi),Disporopsis_fuscopicta),Disporopsis_longifolia),((Polygonatum_cirrhifolium,Polygonatum_stenophyllum),Polygonatum_sp.)),(((((Maianthemum_bicolor,Maianthemum_japonicum),Maianthemum_henryi),Maianthemum_bifolium),Dracaena_cambodiana),(Nolina_atopocarpa,(((Convallaria_majalis,Convallaria_pseudomajalis),Convallaria_keiskei),(Speirantha_gardenii,((Reineckea_carnea,((Rohdea_aurantiaca,((Rohdea_chinensis,Rohdea_delavayi),Rohdea_japonica)),Rohdea_wattii)‘>.0251<.0786’),((Tupistra_grandistigma,Tupistra_muricata),((Aspidistra_zhangii,Aspidistra_erecta),((((Aspidistra_connata,Aspidistra_punctatoides),Aspidistra_subrotata),(Aspidistra_alternativa,Aspidistra_punctata)),(((Aspidistra_sichuanensis,(Aspidistra_zongbayi,Aspidistra_leshanensis)),Aspidistra_fimbriata),(Aspidistra_longipedunculata,(((((((((((Aspidistra_minutiflora,Aspidistra_retusa),Aspidistra_elatior),Aspidistra_hainanensis),(Aspidistra_sutepensis,Aspidistra_longifolia)),(Aspidistra_caespitosa,Aspidistra_omeiensis)),Aspidistra_yingjiangensis),(Aspidistra_arnautovii,Aspidistra_linearifolia)),((Aspidistra_cavicola,(Aspidistra_formosa,(((Aspidistra_obconica,(Aspidistra_bicolor,Aspidistra_nankunshanensis)),Aspidistra_dolichanthera),Aspidistra_papillata))),Aspidistra_quadripartita)),Aspidistra_revoluta),(Aspidistra_fungilliformis,Aspidistra_claviformis)),(Aspidistra_clausa,((Aspidistra_luodianensis,(((Aspidistra_guangxiensis,Aspidistra_obliquipeltata),Aspidistra_radiata),Aspidistra_longgangensis)),((Aspidistra_superba,Aspidistra_ovatifolia),(Aspidistra_daxinensis,Aspidistra_patentiloba)))))))))))))‘>.0962<.1985’))))‘>.1166<.3106’)‘>.3302<.6386’;
第一行的76代表这颗树有多少个物种 1代表有几颗树
PS：一般的常规建树软件如IQ-TREE和RAxML构建的系统发育树都带有BP值枝长，在使用Mcmctree的时候输入的树文件只能带有时间标点信息，所以要删除这些信息，如果物种少可以进行手动删除，如果物种较多可以使用newick_utils这个软件进行删除。
#下载
conda install newick_utils
#删除枝长和BP值等结点信息
nw_topology -I Plastom_gene_ML.contree > Plastom_gene_ML_un_leaf_infomations.contree
#定跟
nw_reroot Plastom_gene_ML_un_leaf_infomations.contree Asparagus_officinalis > Plastom_gene_ML_un_leaf_infomations_reroot.contree
3.配置mcmctree.ctl文件
Mcmctree不同于BEAST不能进行参数配置界面的可视化，所以需要通过修改mcmctree.ctl进行配置。
seed = -2
seqfile = Plastom_gene_mafft.phy *用于估算分歧时间的序列，就是建树所用到的比对后的序列，必须要是phylip格式的
treefile = Plastom_gene_ML_un_leaf_infomations_reroot.contree *带有时间校准点的树文件
outfile = out *输出结果的文件夹
ndata = 1 *输入多序列比对的序列个数
seqtype = 0 * 0: nucleotides; 1:codons; 2:AAs *多序列比对的数据类型：1，核酸；2，密码子；3，氨基酸
usedata = 2 * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV *是否使用我们指定的多序列比对序列去进行分歧时间的估算，0，不使用，虽然也会得到结果，但是结果不准确，一般不使用；1，正常使用多序列进行分析，一般情况下设置为1；2，进行正常的分析，但是要使用到in.BV文件，是当该参数设置为3的时候生成的文件。
clock = 1 * 1: global clock; 2: independent rates; 3: correlated rates *使用的分子种模型1，所有的物种/分支有着相同的进化速率；2，各分支的物种有着独立的进化速率，并且符合正态分布；3，和2类似但是对进化速率的评估方式不同
RootAge = <1.0 * safe constraint on root age, used if no fossil for root. *对进化树的根进行时间设点
model = 0 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 *碱基替代模型，可以通过软件来进行估算
alpha = 0.5 * alpha for gamma rates at sites *一般设置为0.5
ncatG = 5 * No. categories in discrete gamma *一般为5
cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)? *是否要去除不确定字符，例如一些不确定碱基N
BDparas = 1 1 0.1 * birth, death, sampling *一般不需要设置，但是当带有时间标定的树中的时间单位发生改变的时候要进行对应的修改，例如100Mya设置成1Mya的时候，要设置成.01 .01 0.1
kappa_gamma = 6 2 * gamma prior for kappa *kappa的GAMMA分布参数
alpha_gamma = 1 1 * gamma prior for alpha *alpha的GAMMA分布参数
rgene_gamma = 2 20 1 * gamma prior for overall rates for genes *设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数：alpha=2、beta=20、初始平均替换率为每100million年（取决于输入有根树文件中的时间单位）1个替换。若时间单位由100Myr变换为1Myr，则要设置成"2 2000 1"。总体上的平均进化速率为：2 / 20 = 0.1 个替换 / 每100Myr，即每个位点每年的替换数为 1e-9
sigma2_gamma = 1 10 1 * gamma prior for sigma^2 (for clock=2 or 3) **设置所有位点进化速率取对数后方差（sigma的平方）的Dirichlet-GAMMA分布参数：alpha=1、beta=10、初始方差值为1。当clock参数值为1时，表示使用全局的进化速率，各分枝的进化速率没有差异，即方差为0，该参数无效；当clock参数值为2时，若修改了时间单位，该参数值不需要改变；当clock参数值为3时，若修改了时间单位，该参数值需要改变。
finetune =1: .1 .1 .1 .1 .1 .1 * auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr *冒号前的值设置是否自动进行finetune，一般设置成1，然程序自动进行优化分析；冒号后面设置各个参数的步进值：times, musigma2, rates, mixing, paras, FossilErr。由于有了自动设置，该参数不像以前版本那么重要了，可能以后会取消该参数。
print = 1 *是否打印取样信息，0，不打印；1，打印除了分支进化速率的所有信息；2，打印全部信息
burnin = 20000 *将前20000次的取样信息进行舍弃，随后在进行取样并统计
sampfreq = 100 *取样频率
nsample = 10000 *取样次数
*** Note: Make your window wider (100 columns) before running the program.
接下来就是运行Mcmctree进行分歧时间的估算
mcmctree mcmctree.ctl
最后带有分歧时间的系统发育树就会输出在你设置的输出文件下面，名为Figtree.tre，可以通过Figtree进行可视化。最后，可以通过将mcmc.txt（取样信息）输入进Tracer软件进行ESS统计，查看是否收敛，如果所有参数的ESS＞200，则一般认为结果已经收敛，物种的分期时间准确。以下是Tracer的下载网址。
(https://beast.community/tracer)
