1.构建常规的有根树

利用单拷贝基因的CDS，构建获得 newick格式的进化树，构建有根树

2.化石时间/生成phy文件

2.1 化石时间

获得fossil calibration time：
通过网站http://www.timetree.org/ ，该网站根据多篇文献支持提供两两物种间的分化时间，并给出置信度范围，单位是Mya(million years ago)。另外一个可以查询分化时间的网站是https://fossilcalibrations.org/ ，可以互相参考。

在已经构建好的树中添加化石时间，主要形式是 '>0.22<6.75'，'>'后接分化时间的最小值，'<'后接最大值，注意单位是100Mya

7 1
((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);

(
	( 
		 (
			   (
					(chimpanzee, bonobo) 
					,human
				)'>.06<.08'
				,gorilla
		),   
		(orangutan, sumatran)
	)'>.12<.16' 
	,gibbon                 
);

(chimpanzee, bonobo)) '>.06<.08', 表示human与这俩的分化时间

2.2 获取phy文件

perl -lane 'BEGIN{print "7 67364"};if(/>(.*)/){printf $1}else{print"    $_"}'  Single_gene_seq.aln > input.phy
7 67364为7个物种,比对后是65146个碱基,根据实际来更改

3.mcmctree运行

3.1 配置文件

cp /share/nas6/xul/project/chloroplast/GP-20200702-2068_henankejidaxue_nongxueyuan_2samples_zhiwu_yelvti/analysis1/28samples/ref/gene/phytree2/mcmctree.ctl ./

根据实际情况修改

* infile
seqfile = input.phy # 生成的phy文件
treefile = phytree.nwk # 修改的tree文件

* outfile
outfile = tmp1.out
mcmcfile = mcmc.txt 

seqtype = 0	#数据类型；0，表示核酸数据；1，表示密码子比对数据；2，表示氨基酸数据。

ndata = 1	#设置输入的多序列比对的数据个数

usedata = 1   #设置是否利用多序列比对的数据：0，表示不使用多序列比对数据，则不会进行likelihood计算，虽然能得到mcmc树且计算速度飞快，但是其分歧时间结果是有问题的；1，表示使用多序列比对数据进行likelihood计算，正常进行MCMC，是一般使用的参数; 2，进行正常的approximation likelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的in.BV文件。该文件是使用usedata = 3参数生成的out.BV文件重命名而来的。

clock = 2   #设置分子钟方法：1，global clock方法，表示所有分支进化速率一致； 2，independent rates方法，各分支的进化速率独立且进化速率的对数log(r)符合正态分布;  3，correlated rates方法，和方法2类似，但是log(r)的方差和时间t相关。

RootAge = '<2'	#设置root节点的分歧时间，一般设置一个最大值。

model = 4        #设置碱基替换模型：0，JC69；1，K80；2，F81；3，F84；4，HKY85。

noisy = 3

alpha = 0.5      #核酸序列中不同位点，其进化速率不一致，其变异速率服从GAMMA分布。一般设置GAMMA分布的alpha值为0.5。

ncatG = 5     #设置离散型GAMMA分布的categories值。

kappa_gamma = 6 2     #设置kappa（转换/颠换比率）的GAMMA分布参数。

alpha_gamma = 1 1      #设置GAMMA形状参数alpha的GAMMA分布参数.

rgene_gamma = 2 20 1   #设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数：alpha=2、beta=20、初始平均替换率为每100million年（取决于输入有根树文件中的时间单位）1个替换。若时间单位由100Myr变换为1Myr，则要设置成"2 2000 1"。

sigma2_gamma = 1 10 1   #设置所有位点进化速率取对数后方差（sigma的平方）的Dirichlet-GAMMA分布参数：alpha=1、beta=10、初始方差值为1。

print = 1       #设置打印mcmc的取样信息：

burnin = 20000     #将前20000次迭代burnin后，再进行取样（即打印出该次迭代计算的结果信息，各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等）。

sampfreq = 2      #每10次迭代则取样一次。

nsample = 100000    #当取样次数达到该次数时，则取样结束（程序也将运行结束）。

finetune = 1: .1 .1 .1 .1 .1 .1   #冒号前的值设置是否自动进行finetune，一般设置成1，然后程序自动进行优化分析；

3.2 运行

/home/pub/biosoft/paml/paml-4.10.6/bin/mcmctree mcmctree.ctl
/share/nas6/zhouxy/biosoft/paml/4.9/bin/mcmctree mcmctree.ctl

比较耗时间的步骤主要在于取样的百分比进度：

第一列：取样的百分比进度。
第2~6列：参数的接受比例。一般，其值应该在30%左右。20~40%是很好的结果，15~70%是可以接受的范围。若这些值在开始时变动较大，则表示burnin数设置太小。
第7~x列：各内部节点的平均分歧时间，第7列则是root节点的平均分歧时间。若有y个物种，则总共有y-1个内部节点，从第7列开始后的y-1列对应这些内部节点。
倒数第3~x列：r_left值。若输入3各多序列比对结果，则有3列。x列的前一列是中划线 - 。
倒数第1~2列：likelihood值和时间消耗。

屏幕信息最后，给出各个内部节点的分歧时间(t)、平均变异速率(mu)、变异速率方差(sigma2)和r_left的Posterior信息：
均值(mean)、95%双侧置信区间(95% Equal-tail CI)和95% HPD置信区间(95% HPD CI)等信息。
此外，倒数第二列给出了各个内部节点的Posterior mean time信息，可以用于收敛分析。

在当前目录下，生成几个主要结果：
FigTree.tre    生成含有分歧时间的超度量树文件

mcmc.txt       MCMC取样信息，包含各内部节点分歧时间、平均进化速率、sigma2值等信息，可以在Tracer软件中打开。
通过查看各参数的ESS值，若ESS值大于200，则从一定程度上表示MCMC过程能收敛，结果可靠。

out.txt        包含由较多信息的结果文件。例如，各碱基频率、节点命名信息、各节点分歧时间、进化速率和进化树等。

4.生成进化树图片

假设有根树为phytree.nwk,存在tmp1.out 文件

1.
sed -i '1d' phytree.nwk #删除phytree.nwk文件的第一行
cp phytree.nwk input.tree #替换tree里的名字
touch desc #生成一个名为desc的空文件

2.
perl /home/cp_basic/program/draw_time_tree.pl  input.tree desc tmp1.out #测试环境中的程序位置，此命令会生成divtree.newick
perl /share/nas6/zhouxy/modules/MCMCtree/00.bin/03.divergence_time/bin/draw_time_tree.pl  input.tree desc tmp1.out  #生产环境中的程序位置
替换名字

sed -i '2d' divtree.newick #删除divtree.newick的第二行

3.
perl /home/cp_basic/program/draw_tree.pl divtree.newick desc -cali input.tree > divtree.svg #新服务器的程序位置
perl  /share/nas6/liuyj/00.Reseq/Population_evolution/00.bin/Modules/08.MCMCtree/00.bin/03.divergence_time/bin/draw_tree.pl divtree.newick desc -cali input.tree > divtree.svg #日常用的程序位置

perl /share/nas6/xul/program/nuclear/tree/draw_time_tree.pl -i tmp1.out -cali input.tree -o divtree.svg #有点问题

生成的svg需要去改一下参数
修改最前面出现的两个width,改成1200

svg2xxx -t png -dpi 600 divtree.svg && svg2xxx -t pdf divtree.svg