06进化树流程-分歧时间

一.构建常规的有根树

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

perl /share/nas6/pub/pipline/genome-assembly-seq/chloroplast-genome-seq/v1.2/phytree/raxml_phytree.pl -i gene_seq.trim.aln -aln -o time_phytree -r Danio_rerio_NC_002333 &

python3 /share/nas1/yuj/script/chloroplast/phytree/simplify_newick.py  time_phytree/sample.genome.nwk   phytree.nwk

二.化石时间/生成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 /share/nas6/xul/program/mt2/phytree/gene_tree/src/fasta2line.pl -i gene_seq.trim.aln -o gene_seq.trim.aln.line # 转换成一行序列
ir.py -l -i gene_seq.trim.aln.line # 获取长度

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

三.mcmctree运行

3.1 配置文件

cp /share/nas2/yuj/project/2024/mitochondrion/GP-20240614-8555_20241021/QA2_20241119/ref_tre/gene/mcmctree.ctl ./ # 一般不用改了

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        包含由较多信息的结果文件。例如，各碱基频率、节点命名信息、各节点分歧时间、进化速率和进化树等。

3.3 FAQ

1.
Error: you should specify # seqs in the tree file..

进化树少第一行，个数+棵树

四.生成进化树图片

假设有根树为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 #日常用的程序位置

生成的svg需要去改一下参数
修改最前面出现的两个width,改成1200
svg2xxx -t png -dpi 600 divtree.svg && svg2xxx -t pdf divtree.svg

或者，使用的校正节点很多时，用下面的程序直接生成
perl /share/nas6/xul/program/nuclear/tree/draw_time_tree.pl -i tmp1.out -cali input.tree -o divtree.svg

perl /share/nas1/yuj/program/nuclear/tree/draw_time_tree.pl -i tmp1.out -cali input.tree -o divtree.svg
# 修改字体为arial
# 去除 95％ CI
# 分别绘制物种名和登录号
# 加粗测序物种

五.参数

seed：这设置了程序所使用的随机种子。当设置为一个负整数时（比如-1，就像上面给出的例子），程序将使用计算机的当前时间来设置种子，所以每次你运行MCMCTree时，它将从一个不同的种子开始，MCMC的结果看起来也会不同。如果你需要可重复的结果，你可以将种子设置为一个奇数或一个偶数。你应该使用这个选项，并至少运行两次程序，以确定运行的结果非常相似。

sefile和trefile：这分别是序列比对文件和树文件的名称。

mcmfile：关于分析的参数的MCMC运行报告。

oufile：一旦程序完成，它将把结果的摘要写到这个文件中。

ndata：比对文件中分区的数量。在我们的例子中，有三个分区，分别对应于线粒体蛋白质中的三个密码子位置。

usedata：当设置为1时，似然函数以正常方式计算，MCMC分析正常进行。当usedata=0时，不计算似然函数（它被设置为1），所以只计算先验。当usedata=2和=3时，将进行近似似然计算和分支长度的ML估计。它们可用于分析核苷酸、氨基酸和密码子序列；分别使用核苷酸、氨基酸和密码子替换模型。它规定了近似值计算，输入（梯度和Hessian矩阵等）在ﬁle中。
     * 0：no data不使用，不会进行likelihood估算，会快速得到mcmc树，但分歧时间不可用; 
     * 1：seq like，使用多序列比对数据进行likelihood估算，正常进行mcmc; usedata=1时model无法选择；
     * 2：normal 进行正常的approximation likelihood分析，不读取多序列比对数据，直接读取当前目录的in.BV文件，in.BV是由usedata = 3时生成的out.BV重命名得来；
     * 3：程序利用多序列比对数据调用baseml/codeml命令对数据进行分析，生成out.BV文件。

seqtype：变量，为核苷酸序列设置0，为密码子序列设置1，为氨基酸序列设置2。

clock：要使用的时钟模型。这里我们将使用独立速率模型（clock=2），其中速率遵循对数正态分布（也就是说，速率的对数是正态分布）。

RootAge：如果树文件中没有提供，则使用根的校准。这里我们使用"<1.0"即所有类人猿最近的共同祖先最大限制为100Myr。

model, alpha and ncatG：要使用的替换模型。在这个例子中，我们将使用JC69（它的计算非常快）。因为α=0，所以我们在本例中不使用速率变化的gamma模型。

cleandata=0：意味着在似然计算中，比对gap和模糊字符将被视为缺失数据。= 1意味着在分析之前，任何至少有一个序列有比对gap或模糊字符的位点都将被删除。
这个变量用于usedata=1和3，如果usedata=2，则没有影响。

BDparas：控制出生-死亡过程的参数。出生-死亡过程用于构建树中没有化石标定的节点的时间先验。这里我们使用默认的1 1 0.1，它产生了统一的节点年龄先验。

kappa_gamma和alpha_gamma：设置kappa（转换/颠换比率）的GAMMA分布参数；设置GAMMA形状参数alpha的GAMMA分布参数。

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

sigma2_gamma：设置所有位点进化速率取对数后方差（sigma的平方）的Dirichlet-GAMMA分布参数：alpha=1、beta=10、初始方差值为1。当clock参数值为1时，表示使用全局的进化速率，各分枝的进化速率没有差异，即方差为0，该参数无效；当clock参数值为2时，若修改了时间单位，该参数值不需要改变；当clock参数值为3时，若修改了时间单位，该参数值需要改变。

print：如果设置为1，MCMC的输出和结果的摘要将被写入硬盘（MCMC被写入mcmc.out文件，摘要被写入上述设置的outﬁle文件）。如果设置为2，每个基因座（分区）的分支率将被附加到输出中。请注意，对于宽松的时钟模型（时钟=2或3），输出的速率相当大，可能需要更多的时间将它们全部写入输出文件（打印=2）。因此，根据你使用的模型，你可以决定 "print "是设置为1还是2。如果设置为0，结果将只被打印到屏幕上。你不希望这样。

burnin、sampfreq和nsample：在我们的例子中，程序将放弃前2000次迭代作为burn-in，然后它将每10次迭代采样一次，直到收集到20000个样本。总的来说，MCMC将运行2,000+10×20,000=202000次迭代。通常情况下，你应该收集10,000到20,000个样本以获得良好的统计总结。大样本量（例如100,000）往往会浪费大量的硬盘空间，提供很少的统计改进。程序也可能需要很长的时间来总结结果。如果你需要增加MCMC的长度（以提高收敛性），增加sampfreq，但保持nsample在一个合理的值。
运行程序

种子值（seed）

seed = -1

作用：设置程序使用的随机数种子。
说明：
- 值为负数时（如 -1），程序会根据计算机当前时间生成种子，每次运行结果可能不同。
- 值为正数或固定负数时，可以确保结果可重复。
建议：至少运行两次程序，确保结果一致。

输入文件

seqfile = mtCDNApri123.txt
treefile = mtCDNApri.trees

seqfile：序列比对文件，包含核苷酸或蛋白质序列。
treefile：系统发育树文件，包含树结构和化石校准信息。

输出文件

mcmcfile = mcmc.txt
outfile = out.txt

mcmcfile：记录 MCMC 运行参数。
outfile：保存运行结果的摘要。

数据分区数量

ndata = 3

作用：指定序列文件中分区数量。
说明：本例中分为 3 个部分，对应密码子的第 1、第 2 和第 3 位点。

序列类型

seqtype = 0 * 0: 核苷酸；1: 密码子；2: 氨基酸

作用：设置输入序列的类型。
- 0：核苷酸序列。
- 1：密码子序列。
- 2：氨基酸序列。
本例设置：seqtype = 0。

数据使用方式

usedata = 1 * 0: 无数据；1: 序列；2: 近似计算；3: 输出.BV（输入.BV）

作用：指定程序如何使用数据。
- 0：不使用数据，仅计算先验分布。
- 1：使用完整序列，进行标准计算。
- 2：使用近似方法计算。
- 3：使用梯度和海森矩阵文件。
本例设置：usedata = 1。

分子钟模型

clock = 2 * 1: 全局时钟；2: 独立时钟；3: 相关速率

作用：选择分子钟模型。
- 1：全局时钟模型，假定所有谱系的进化速率相同。
- 2：独立速率模型，假定谱系速率独立。
- 3：相关速率模型，假定谱系速率相关。
本例设置：clock = 2。

根节点年龄校准

RootAge = '<1.0' * 对根年龄的安全限制

作用：设置树根校准范围。
说明：如果树文件中未指定校准点，则此参数提供限制值。本例中，设置根年龄最大值为 1.0（单位：亿年）。

替代模型

model = 0 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85

作用：指定替代模型。
选项：
- 0：JC69 模型。
- 1：K80 模型。
- 2：F81 模型。
- 3：F84 模型。
- 4：HKY85 模型。
本例设置：model = 0。

Gamma 分布参数

alpha = 0 * alpha for gamma rates at sites
ncatG = 5 * No. categories in discrete gamma

alpha：设置 Gamma 分布的形状参数。
- 0：不使用 Gamma 模型。
ncatG：离散 Gamma 分布类别数量。

数据清理

cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)?

作用：是否移除不确定或模糊的数据。
- 0：保留模糊数据。
- 1：移除含缺失值的数据。

出生-死亡过程参数

BDparas = 1 1 0.1 * birth, death, sampling

作用：设置出生-死亡模型参数。
- 第一个值：出生率。
- 第二个值：死亡率。
- 第三个值：采样比例。

Gamma 先验分布

kappa_gamma = 6 2 * gamma prior for kappa
alpha_gamma = 1 1 * gamma prior for alpha

kappa_gamma：转换/颠换比率（κ）的 Gamma 先验分布。
alpha_gamma：Gamma 模型形状参数（α）的先验分布。

基因速率和速率漂移

rgene_gamma = 2 20 1 * gammaDir prior for rate for genes
sigma2_gamma = 1 10 1 * gammaDir prior for sigma^2

rgene_gamma：基因速率的 Gamma-Dirichlet 分布。
sigma2_gamma：速率漂移参数（σ²）的 Gamma-Dirichlet 分布。

自动微调步长

finetune = 1: .1 .1 .1 .1 .1 .1 * auto (0 or 1)

作用：控制 MCMC 步长是否自动调整。
本例设置：1（启用自动调整）。

输出控制

print = 1 * 0: no mcmc sample; 1: everything except branch rates; 2: everything

作用：控制输出详细程度。
- 0：仅屏幕输出。
- 1：输出除分支速率外的所有结果。
- 2：输出所有结果。

MCMC 参数

burnin = 2000
sampfreq = 10
nsample = 20000

burnin：烧入阶段迭代次数，前 2000 次丢弃。
sampfreq：每隔多少次迭代保存一次样本。
nsample：总共保存的样本数量。

六.方法说明

使用Paml v4.9中的MCMCTree方法，采用松弛分子钟方法估算了分歧时间（Yang, 2007）。
校准使用了从时间树数据库（http://timet.ree.org/）中获得的8个分歧时间点。第一个校准点对应于Mene maculata + Xiphias gladius的分歧（84.3–57.4百万年）(http://timet.ree.org/)。……最后一个校准点对应于Epinephelus akaara + Plectropomus leopardus的分歧（62.3-22.3百万年）（http://timet.ree.org/）。
为了确保收敛，我们检查了来自两个独立运行的链，在Tracer v.1.7（Rambaut et al., 2018）中评估有效样本量（ESS）值是否大于200，以确保每个参数从后验分布中适当采样。
生成的时间树通过Figtree v.1.4.4（http://tree.bio.ed.ac.uk/software/figtree/）查看和编辑。


## 文章来源 https://www.frontiersin.org/journals/genetics/articles/10.3389/fgene.2024.1414074/full
Divergence time estimation
Divergence time was computed across each dataset (UCE, BUSCO, and AHE) at 50% species occupancy. We estimated the divergence time utilizing a relaxed molecular clock method using MCMCTree in Paml v4.9 (Yang, 2007). Calibration was performed using three divergence time points obtained from the timetree database (http://timet.ree.org/) and literature. The first calibration point corresponds to the divergence of Dacus + Zeugodacus (86.3–59.3 Mya) (Krosch et al., 2012). The second calibration point represents the most recent common ancestor of subgenus Tetradacus (30.9–12.4 Mya) (Krosch et al., 2012). The final calibration point is the divergence of Drosophila (38–62 Mya) (http://timet.ree.org/). To ensure convergence, chains from two independent runs were checked in Tracer 1.7 (Rambaut et al., 2018) to assess the effective sample size (ESS) values above 200, indicating appropriate sampling from the posterior distribution of each parameter. The resulting time trees were viewed and edited using Figtree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).