image.png
1. 多序列比对文件转phylip
trimal -in input.fa -out out.phy -phylip_paml
树文件，在orthofinder下面由物种树直接拷贝，删除枝长。
sed -i 's/:[^,[]]\+//g' rooted.txt ^表示非
在开头加序列条数和树的个数，以及化石点，如果根部没有时间，必须在ctl文件中指定一个根的范围。
例如：
4 1
(((1,2),3)‘>2<3’,4)'<4.5'; ##给的时间范围
或者
4 1
(((1,2),3)‘@0.2’,4); ##给定时间点，根部没有时间点，则需要做mcmctree.ctl文件中的RootAge处给出范围。
3.修改ctl文件
seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs
usedata = 3 * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
RootAge = * safe constraint on root age, used if no fossil for root. 如果根部没有化石点则加上范围，比如<5.4表示根要小于5.4个mya
4. 运行第一次mcmc，生成out.BV，利用codeml对生成的temp.ctl其进行修改重新计算。
cp temp.ctl codeml.ctl (核酸则为baseml)
第一添加‘clock = 1 ’ ##生成gama分布的先验概率
第二修改 getSE = 2 为getSE = 0
第三修改model = 0 为 model = 2
第四修改model需要的模型文件aaRatefile = 为aaRatefile = wag.dat #这个文件在paml软件目录下
codeml codeml.ctl
##查看计算的替换速率
grep -A 2 'Substitution' temp0001.out
Substitution rate is per time unit
0.303717 ##这个就是替换速率
5.修改mcmctree.ctl中的rgene_gama参数以及usedata参数，再次运行mcmctree生成out.BV
usedata = 3
rgene_gama = 2 2 修改改为 regene = 2 6 #2除以6=0.3
mcmctree mcmctree.ctl
6.把out.BV修改为in.BV，修改mcmctree.ctl中的usedata = 2 ，再次运行mcmctree
sed 's/usedata = 3 /usedata = 2/g' mcmctree.ctl
mv out.BV in.BV
mcmctree mcmctree.ctl
##记得这一步计算两次，如果两次结果差别很小，且ESS>200，则认为结果收敛。
mv out.BV in.BV
mcmctree mcmctree.ctl
7.检查log文件
5% 0.32 0.30 0.33 0.30 0.26 4.931 0.547 4.894 3.951 3.534 0.161 - 0.106 0.146 -334.9 5:04
10% 0.31 0.30 0.32 0.29 0.26 4.954 0.545 4.917 3.971 3.556 0.162 - 0.104 0.143 -347.6 5:27
15% 0.31 0.31 0.32 0.29 0.26 4.949 0.544 4.912 3.959 3.543 0.161 - 0.100 0.143 -317.3 5:50
20% 0.31 0.30 0.31 0.30 0.26 4.943 0.545 4.907 3.959 3.543 0.161 - 0.101 0.142 -334.4 6:14
25% 0.31 0.30 0.31 0.29 0.26 4.952 0.546 4.917 3.972 3.555 0.161 - 0.100 0.142 -344.4 6:37
30% 0.31 0.30 0.32 0.29 0.26 4.954 0.547 4.918 3.974 3.557 0.161 - 0.101 0.142 -315.9 7:00
35% 0.31 0.30 0.32 0.29 0.26 4.954 0.547 4.918 3.974 3.558 0.161 - 0.099 0.142 -317.6 7:23
40% 0.31 0.30 0.31 0.29 0.25 4.952 0.547 4.917 3.973 3.558 0.161 - 0.100 0.141 -368.3 7:46
45% 0.31 0.30 0.31 0.29 0.26 4.956 0.547 4.921 3.974 3.559 0.161 - 0.101 0.141 -316.2 8:09
50% 0.31 0.30 0.31 0.29 0.26 4.960 0.548 4.924 3.977 3.562 0.161 - 0.101 0.141 -327.2 8:32
55% 0.31 0.30 0.31 0.29 0.25 4.958 0.548 4.922 3.976 3.561 0.161 - 0.101 0.141 -311.9 8:55
60% 0.31 0.31 0.31 0.29 0.25 4.959 0.548 4.923 3.975 3.561 0.161 - 0.101 0.141 -330.9 9:18
65% 0.31 0.30 0.31 0.29 0.25 4.957 0.548 4.922 3.974 3.559 0.161 - 0.100 0.142 -338.8 9:41
70% 0.31 0.30 0.31 0.29 0.25 4.957 0.548 4.921 3.972 3.557 0.161 - 0.100 0.142 -334.2 10:04
75% 0.31 0.30 0.31 0.29 0.25 4.957 0.549 4.921 3.973 3.559 0.161 - 0.101 0.142 -316.2 10:27
80% 0.31 0.31 0.31 0.29 0.25 4.958 0.549 4.923 3.975 3.560 0.161 - 0.101 0.142 -314.1 10:51
85% 0.31 0.31 0.31 0.29 0.25 4.960 0.550 4.924 3.975 3.561 0.161 - 0.101 0.141 -331.4 11:14
90% 0.31 0.30 0.31 0.29 0.25 4.961 0.550 4.925 3.977 3.562 0.161 - 0.101 0.142 -331.5 11:37
95% 0.31 0.31 0.31 0.29 0.25 4.957 0.549 4.921 3.973 3.559 0.161 - 0.101 0.142 -345.7 12:00
100% 0.31 0.31 0.31 0.29 0.25 4.956 0.548 4.920 3.971 3.557 0.161 - 0.101 0.142 -320.1 12:23
第一确保第二列到第六列在30%左右，如果小于15大于60则考虑修改burnin次数，增加迭代次数。
第二确保第7列到倒数第3列在在最后几轮计算要稳定下来（预估的几个节点分化时间），不会出现增大的现象，如果结束时还未稳定，增加迭代次数和nsample次数
检查输出文件outfile
第一次：
t_n50 0.0399 (0.0270, 0.0519) (0.0273, 0.0521) 0.0248 (Jnode 42)
t_n51 0.0027 (0.0013, 0.0044) (0.0012, 0.0042) 0.0030 (Jnode 41)
t_n52 4.5199 (3.1873, 5.2593) (3.3363, 5.3380) 2.0016 (Jnode 40)
t_n53 4.2802 (3.0241, 4.9904) (3.1261, 5.0532) 1.9271 (Jnode 39)
t_n54 1.3992 (0.9847, 1.6928) (1.0122, 1.7086) 0.6965 (Jnode 38)
t_n55 0.9347 (0.6529, 1.1481) (0.6740, 1.1593) 0.4853 (Jnode 37)
t_n56 0.3230 (0.2219, 0.4071) (0.2292, 0.4127) 0.1835 (Jnode 36)
t_n57 2.2284 (1.5726, 2.6683) (1.6091, 2.6898) 1.0807 (Jnode 35)
t_n58 0.8284 (0.5737, 1.0397) (0.5894, 1.0525) 0.4631 (Jnode 34)
t_n59 1.3289 (0.9244, 1.6125) (0.9576, 1.6269) 0.6694 (Jnode 33)
t_n60 1.0425 (0.7217, 1.2734) (0.7516, 1.2907) 0.5391 (Jnode 32)
t_n61 0.5531 (0.3783, 0.6975) (0.3918, 0.7071) 0.3153 (Jnode 31)
第二次
t_n50 0.0398 (0.0267, 0.0520) (0.0267, 0.0520) 0.0252 (Jnode 42)
t_n51 0.0027 (0.0013, 0.0044) (0.0012, 0.0043) 0.0031 (Jnode 41)
t_n52 4.5069 (3.1599, 5.2554) (3.3138, 5.3393) 2.0255 (Jnode 40)
t_n53 4.2680 (2.9935, 4.9956) (3.1746, 5.1140) 1.9394 (Jnode 39)
t_n54 1.3958 (0.9725, 1.6951) (1.0121, 1.7215) 0.7094 (Jnode 38)
t_n55 0.9338 (0.6449, 1.1422) (0.6702, 1.1577) 0.4874 (Jnode 37)
t_n56 0.3217 (0.2179, 0.4057) (0.2229, 0.4093) 0.1864 (Jnode 36)
t_n57 2.2226 (1.5485, 2.6566) (1.6267, 2.7178) 1.0911 (Jnode 35)
t_n58 0.8263 (0.5657, 1.0346) (0.5698, 1.0366) 0.4668 (Jnode 34)
t_n59 1.3248 (0.9226, 1.6068) (0.9629, 1.6345) 0.6716 (Jnode 33)
t_n60 1.0404 (0.7233, 1.2707) (0.7444, 1.2859) 0.5415 (Jnode 32)
t_n61 0.5514 (0.3779, 0.6970) (0.3831, 0.7008) 0.3177 (Jnode 31)
比较两次结果相关性。吧这个文件放到tracer中，确保ESS>200。
网友评论