准备文件：

1，化石标定的物种关系树 （参考已发表的文献） 改名为input.tree
2，mcmctree.ctl文件 记得修改里面的树文件，序列文件路径
3，序列文件(orthofinder得到的物种数,phylip格式)

#准备文件：
# 1，化石标定的物种关系树 （参考已发表的文献） 改名为input.tree
#2，mcmctree.ctl文件 记得修改里面的树文件，序列文件路径
#3，序列文件(orthofinder得到的物种数,phylip格式)
#软件安装
conda install paml

#1，把fasta格式序列改成菲利普格式
cat SpeciesTreeAlignment.fa |tr '\n' '\t' |sed 's/>/\n/g' |sed 's/\t/        /' |sed 's/\t//g' |awk 'NF > 0' > supergene.phy.tmp

#2在前面加上一部分内容，得到最终菲利普文件
awk '{print " "NR" "length($2)}' supergene.phy.tmp|tail -n 1 | cat - supergene.phy.tmp > supergene.phy

####改成菲利普格式方法2
  trimal -in SpeciesTreeAlignment.fa -out supergene.phy -phylip_paml

#对树进行修改
sed 's/:[^,)(]\+//g' SpeciesTree_rooted.txt |sed 's/)1/)/g' > input.tree
并加上进化时间，时间点的选择。最远的，中间的，最近的各一个
在上面加上物种数量，1棵树

sed -i 's/\.pro//g' input.tree 
sed -i 's/\.pro//g' supergene.phy

##step 1,估算位点替换率
#将mcmctree.ctl 配置文件中seqtype 设置为2；将usedata 设置为3  运行mcmctree.ctl
mcmctree mcmctree.ctl

#将生成的tmp0001.ctl 拷贝为codeml.ctl 添加 clock = 1: 修改getSE为0
cp tmp0001.ctl codeml.ctl
echo "clock = 1" >> codeml.ctl
把getSE = 2 改为getSE = 0 或 sed -i 's/\(getSE.*=\s*\d/getSE = 0/' codeml.ctl

#将生成的tmp0001.trees 文件进行修改，定根，添加化石时间（'@中间值'）
#运行codeml,
codeml codeml.ctl
#查看tmp0001.out结果中的替换速率，替换速率在“Substitution rate is per time unit后
tail -n 20 tmp0001.out
#知道替换速率后，修改rgene_gamma =2 x 中x的值，上述都是为了修改x。所以修改成功后重新运行



# step2,使用近似然法计算分化时间,即重新开始
#调整mcmctree.ctl中的rgene_gamma 第二个参数b,使得a/b约等于之前得到的替换率；将usedata设置为3 运行mcmctree.ctl
mcmctree mcmctree.ctl > run.log
#将生成的out.BV 文件改名为in.BV
mv out.BV in.BV

# 将mcmctree.ctl中usedata 设置为2，运行
mcmctree mcmctree.ctl > run1.log

#最终结果仍为figtree