mcmctree需要3个输入文件
- 多序列比对的phy格式文件（注意：这个phy和一般软件的phy不太一样）使用脚本把多序列比对的fa格式转为phy格式
- 准备包含2个物种之间的化石时间的进化树
- mcmctree.ctl运行文件
准备phy格式文件supergene.phy
如果现有多序列比对文件all.fa
cat all.fa |tr '\n' '\t'|sed 's/>/\n/g' |sed 's/\t/ /'|sed 's/\t//g'| awk 'NF > 0' > supergene.phy.tmp
awk '{print " "NR" "length($2)}' supergene.phy.tmp|tail -n 1 | cat - supergene.phy.tmp > supergene.phy
输出的supergene.phy就是本次需要的文件
如果现有的是多序列比对是phy格式all.fa.phy
sed -r 's/(.*) (.*)/\1 \2/g' all.fa.phy >supergene.phy 给序列和id之间添加空格
最终的supergene.phy格式如下
25 202587
speciesA MDTLLRTHSKLEFSVLLHGFSEKACKTHQ
specde4 MDTLLRTHSKLEFSVLLHGFSEKACKTHQ
spec3 MDTLLRTHSKLEFSVLLHGFSEKACKTHQ
spec2
MDTLLRTHSKLEFSVLLHGFSEKACKTHQ
...
spec1 MDTLLRTHSKLEFSVLLHGFSEKACKTHQ
准备进化树文件mcmc.input.tree
从单拷贝基因构建的进化树，可以手动删除除了物种名之外的信息。也可以使用下面的R脚本删除
library(ape)
MLtree <- read.tree("all.fa.Jul12.contree") #读取单拷贝进化树
MLtree$edge.length <- NULL #删除分支长度信息
MLtree$node.label <- NULL
write.tree(MLtree,"MLtree.nobranch.tree") #输出只有物种名的进化树
在手动在第一行添加上物种的数量25和1中间是一个空格。
最终的进化树文件mcmc.input.tree如下：
25 1
((speciesA :'>1.0<1.6', ((spec2, spec3), (specde4, ((((XXA, XXB), ((BBA, (BMW, WBA)), QWE)), ADA), ABA)))), (AKD:'>1.3<2.6', ((ADKD, ADLD), (DD3, ((((HPD, ((DMD2, DMD3), (TMD2, TMD1))), YID), NDD1), LD6D)))):'>3.5<35.1', (spec1));
中间的时间格式是百万年。这个需要从http://www.timetree.org/ 查询两个物种间的分化时间
比如上图搜索玉米和拟南芥的分开时间，可以看到结果是142.1和163.5MYA.
如果是进化树中，则写为
>142.1<163.5。此处的单位是MYA,百万年。看到两者是单子叶和双子叶的分开时间。
准备运行的配置文件mcmctree.ctl
seed = -1 # 表示使用系统当前时间为随机数
seqfile = supergene.phy # 多序列比对文件
treefile = mcmc.input.tree # 带有校准信息有根树
mcmcfile = mcmc.txt #输出mcmc信息文本文件，可以用Tracer软件查看
outfile = out.txt # 输出文件，记录了超度量树和进化速率等参数信息
ndata = 1 # 输入的多序列比对的数据个数，这是密码子位置的数据；如果有一个，则设置为1
seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs #数据类型
usedata = 3 * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV) # 设置是否利用多序列比对的数据：\
#0，表示不使用多序列比对数据，则不会进行likelihood计算，虽然能得到mcmc树且计算速度飞快，但是其分歧时间结果是有问题的；\
#1，表示使用多序列比对数据进行likelihood计算，正常进行MCMC，是一般使用的参数; \
#2，进行正常的approximation likelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的in.BV文件。该文件是使用usedata = 3参数生成的out.BV文件重命名而来的。\
#此外，由于程序BUG，当设置usedata = 2时，一定要在改行参数后加 *，否则程序报错 Error: file name empty.. \
#3，程序利用多序列比对数据调用baseml/codeml命令对数据进行分析，生成out.BV文件。由于mcmctree调用baseml/codeml进行计算的参数设置可能不太好（特别是对蛋白序列进行计算时），\
#推荐自己修该软件自动生成的baseml/codeml配置文件，然后再手动运行baseml/codeml命令，再整合其结果文件为out.BV文件。
clock = 1 * 1: global clock; 2: independent rates; 3: correlated rates # (1) 全局分子钟,适用于近缘物种（两物种分化时间小于20个million分化时间或者序列差异小于5%即近缘物种）\
# (2) 树枝上进化速率服从独立同分布的对数正态分布（推荐使用）
# (3) 树枝上的进化速率自相关。
RootAge = '<54' * safe constraint on root age, used if no fossil for root. # 设置root节点的分歧时间，一般设置一个最大值。
model = 0 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
# 碱基替换模型：
#（0）无参数（适用于近缘物种）
#（1）参数为转换颠换比Kappa
#（2）参数为T，C，A，G的频率
#（3）最复杂的进化模型
#（4）参数为T，C，A，G的频率+kappa（适用于远缘物种）
alpha = 0.5 * alpha for gamma rates at sites
# 核酸序列中不同位点，其进化速率不一致，其变异速率服从GAMMA分布。一般设置GAMMA分布的alpha值为0.5。\
# 若该参数值设置为0，则表示所有位点的进化速率一致。此外，当userdata = 2时，alpha、ncatG、alpha_gamma、kappa_gamma这4个GAMMA参数无效。因为userdata = 2时，不会利用到多序列比对的数据。
ncatG = 5 * No. categories in discrete gamma # 设置离散型GAMMA分布的categories值。
cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)?
BDparas = 1 1 0.1 * birth, death, sampling # 设置出生率、死亡率和取样比例。若输入有根树文件中的时间单位发生改变，则需要相应修改出生率和死亡率的值。例如，时间单位由100Myr变换为1Myr，则要设置成".01 .01 0.1"。
kappa_gamma = 6 2 * gamma prior for kappa #设置kappa（转换/颠换比率）的GAMMA分布参数
alpha_gamma = 1 1 * gamma prior for alpha #设置GAMMA形状参数alpha的GAMMA分布参数
rgene_gamma = 2 20 1 * gammaDir prior for rate for genes #设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数：
sigma2_gamma = 1 10 1 * gammaDir 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，然后程序自动进行优化分析
print = 1 * 0: no mcmc sample; 1: everything except branch rates 2: everything # 设置打印mcmc的取样信息
#0，不打印mcmc结果；
#1，打印除了分支进化速率的其它信息（各内部节点分歧时间、平均进化速率、sigma2值）；
#2，打印所有信息。
burnin = 2000 #将前2000次迭代burnin后，再进行取样（即打印出该次迭代计算的结果信息，各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等）
sampfreq = 10 #每10次迭代则取样一次。
nsample = 20000 # 当取样次数达到该次数时，则取样结束（程序也将运行结束）。
注意一定要修改的几个参数：
seqfile = supergene.phy # 多序列比对文件
treefile = mcmc.input.tree # 带有校准信息有根树
ndata = 1 # 输入的多序列比对的数据个数，这是密码子位置的数据；如果有一个，则设置为1,指的是你的phy文件有多少组密码子的位置，我的是只有1个。一般用单拷贝的都是1个。
seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs #数据类型 因为我用的是蛋白序列比对的，所以此处是2，根据你建树的实际情况选择
usedata = 3 此处设置为3，用于获取输出的out.BV,用于下一步分析。
RootAge = '<0.54' 这个一定要设置一个最大的值，就是你的进化树的根的物种和其他物种的分开时间，如果这个值设置的比里面物种的分开时间还小，则会导致错误。注意单位
开始第一步分析
mcmctree mcmctree.ctl
开始第2步分析
cp mcmctree.ctl mcmctree2.ctl
mv out.BV in.BV
##把原来的3修改为2
sed -i 's/usedata = 3/usedata = 2/' mcmctree2.ctl
#再看一下是不是修改成功了，如果失败了，需要手动修改
grep usedata mcmctree2.ctl
mcmctree mcmctree2.ctl
输出的中间文件的意义的解析
5% 0.40 0.38 0.49 0.52 0.48 192.342 182.765 159.115 145.363 125.924 112.838 - 34.259 0.001 -13238.5 0:03
10% 0.40 0.37 0.49 0.52 0.48 174.265 165.586 144.140 131.675 114.058 102.202 - 31.031 0.001 -13238.4 0:05
15% 0.40 0.37 0.49 0.52 0.48 168.111 159.732 139.043 127.020 110.026 98.588 - 29.934 0.001 -13238.4 0:08
20% 0.40 0.38 0.49 0.52 0.48 163.435 155.295 135.174 123.484 106.959 95.838 - 29.099 0.001 -13238.5 0:10
25% 0.40 0.38 0.49 0.52 0.48 159.276 151.344 131.729 120.334 104.230 93.393 - 28.357 0.001 -13238.5 0:12
30% 0.40 0.37 0.49 0.52 0.48 156.503 148.711 129.431 118.233 102.409 91.761 - 27.860 0.001 -13238.5 0:14
35% 0.40 0.37 0.49 0.52 0.49 154.218 146.540 127.539 116.505 100.914 90.420 - 27.452 0.001 -13238.4 0:16
40% 0.40 0.37 0.49 0.52 0.49 153.412 145.775 126.871 115.894 100.384 89.945 - 27.307 0.001 -13238.4 0:19
45% 0.40 0.37 0.49 0.52 0.49 151.942 144.380 125.654 114.782 99.422 89.082 - 27.045 0.001 -13238.4 0:21
50% 0.40 0.37 0.49 0.52 0.49 150.540 143.050 124.496 113.724 98.506 88.261 - 26.795 0.001 -13238.4 0:23
55% 0.40 0.37 0.49 0.52 0.48 149.618 142.174 123.733 113.027 97.902 87.718 - 26.631 0.001 -13238.5 0:25
60% 0.40 0.37 0.49 0.52 0.49 148.767 141.366 123.028 112.382 97.342 87.216 - 26.479 0.001 -13238.5 0:27
65% 0.40 0.37 0.49 0.52 0.49 148.103 140.736 122.479 111.881 96.907 86.827 - 26.360 0.001 -13238.5 0:30
70% 0.40 0.37 0.49 0.52 0.49 147.605 140.262 122.066 111.503 96.581 86.534 - 26.272 0.001 -13238.5 0:32
75% 0.40 0.37 0.49 0.52 0.49 147.376 140.045 121.876 111.330 96.431 86.399 - 26.230 0.001 -13238.5 0:34
80% 0.40 0.37 0.49 0.52 0.49 146.964 139.652 121.534 111.017 96.160 86.156 - 26.156 0.001 -13238.5 0:36
85% 0.40 0.37 0.49 0.52 0.49 146.665 139.369 121.287 110.791 95.965 85.981 - 26.103 0.001 -13238.5 0:38
90% 0.40 0.37 0.49 0.52 0.49 146.751 139.451 121.357 110.855 96.020 86.031 - 26.119 0.001 -13238.5 0:41
95% 0.40 0.37 0.49 0.52 0.49 146.869 139.564 121.455 110.945 96.097 86.100 - 26.140 0.001 -13238.5 0:43
100% 0.40 0.37 0.49 0.52 0.48 147.188 139.867 121.719 111.186 96.307 86.288 - 26.197 0.001 -13238.5 0:45
- 第1列的百分比表示进度
- 第2到6列表示参数的接受比例，一般是30%左右，20-40%是很好的结果，15-70%是可以接受的范围。如果这些值在刚开始时变化波动大，则说明是burnin数设置太小。
- 第7列到第x列，是各内部节点的平均分歧时间，第7列是root节点的平均分歧时间。如果有y个物种，则总共有y-1个内部节点，从第7列开始后的y-1列对应这些内部节点。
- 倒数第1-2列：likehood值和消耗的时间
Posterior means (95% Equal-tail CI) (95% HPD CI) HPD-CI-width
t_n31 147.1834 (136.5031, 190.0313) (134.5728, 173.8694) 39.2966 (Jnode 58)
t_n32 139.8726 (129.7714, 180.7106) (128.2654, 165.4282) 37.1628 (Jnode 57)
t_n33 121.7225 (112.9277, 157.2861) (111.6496, 144.0022) 32.3526 (Jnode 56)
t_n34 111.1893 (103.1524, 143.6049) (101.9706, 131.6090) 29.6384 (Jnode 55)
t_n35 96.3094 (89.3654, 124.4121) (88.4921, 114.1726) 25.6805 (Jnode 54)
t_n36 86.2896 (80.0571, 111.4612) (79.0442, 101.9939) 22.9497 (Jnode 53)
.......
mu 0.0014 ( 0.0011, 0.0015) ( 0.0012, 0.0016) 0.0004
lnL -13238.4546 (-13246.8960, -13232.0000) (-13246.0200, -13231.3900) 14.6300
mean 147.1834 139.8726 121.7225 111.1893 96.3094 86.2896 71.1097 49.5276 42.3463 37.7958 33.8800 9.6853 35.6002 18.4681 10.0232 38.9533 65.2759 33.7331 29.9372 21.5186 13.4677 9.4206 6.5970 43.0280 31.0316 14.7923 5.9181 4.2413 26.1976 0.0014 -13238.4546
Eff 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0006 0.0005 1.0067
time prior: Birth-Death-Sampling
rate prior: Log-Normal
arcsin transform is used in approx like calculation.
上面给出了各个内部节点的分歧时间，变异速率(mu)，变异速率方差（sigma2）
-t_n36这种以t开头的行，第2列的数字86.2896表示平均时间（86.2896MYA），(80.0571, 111.4612) 是95%双侧置信区间(95% Equal-tail CI)， (79.0442, 101.9939) 是95% HPD置信区间(95% HPD CI)。 22.9497是HPD-CI-width。
-
mu是变异速率。当上面的mcmctree.ctl里设置clock设置为1时，即设置为固定速率，所以此处就只输出1个mu. 上面的示例的结果中，物种间的分歧时间超过20MYA,即不适合再设置为固定速率，所以需要修改clock为2或3.
输出结果解析
FigTree.tre 生成含有分歧时间的超度量树文件
使用figtree打开生成的文件FigTree.tre，即可直接看到进化树，在左侧工具栏勾选对应的栏
绘制出来的图的下面的标尺的时间顺序不对是反的。使用reverse后就都是负数。
还可以使用mcmctreeR这个R包来进行可视化
https://cran.r-project.org/web/packages/MCMCtreeR/MCMCtreeR.pdf
https://github.com/PuttickMacroevolution/MCMCtreeR
library(MCMCtreeR)
MCMCtree <- readMCMCtree("FigTree.tre", forceUltrametric = TRUE, from.file = TRUE)
pdf("mcmctimetee.pdf")
MCMC.tree.plot(phy=MCMCtree, analysis.type="MCMCtree",
plot.type="phylogram", cex.tips=0.5)
dev.off()
输出的结果文件
mcmctimetee.pdf,这个就是MCMCtree.R输出的最终的绘图的结果文件
FigTree.tre 这个就是最终的进化树文件
mcmc.txt 这个文件可以使用Tracker软件（http://beast.bio.ed.ac.uk/tracer） 在windows端打开，查看ESS值。
SeedUsed 这个是随机数文件,设置相同的随机数，会输出相同的结果
out.txt 包含屏幕输出的结果文件。
判断MCMC结果是否达到收敛状态 (这一步是必须的，否则无法证明你的结果是否可靠)
1.使用Tracer打开输出的mcmc.txt文件，如果ESS值大于200，则说明结果良好。如果ESS<200,则还需要调整 burnin,sampfreq, nsample的值，需要更多的迭代次数
2.使用不同的seed值运行两次，比较两次的结果树的内部节点的Height值。
library(coda)
# 读取两次 MCMC 输出
run1 <- read.table("run1/mcmc.txt", header=TRUE)
run2 <- read.table("run2/mcmc.txt", header=TRUE)
# 转换为 mcmc 对象
mcmc1 <- as.mcmc(run1)
mcmc2 <- as.mcmc(run2)
# 合并链
combined <- mcmc.list(mcmc1, mcmc2)
# Gelman-Rubin 诊断（R-hat）
#R-hat (Potential Scale Reduction Factor)：R-hat < 1.1 表示链收敛。 如果 R-hat ≈ 1，说明两次运行结果一致.
gelman.diag(combined)
#查看2条曲线重叠度，如果重叠度高则说明链收敛。
plot(density(run1$t_n1), col="red", main="Node age comparison")
lines(density(run2$t_n1), col="blue")
##检测ACF图的自相关性，理想情况：ACF 在滞后 k=5~20 步内降至接近 0（即 |ACF| < 0.1）。
acf(run1$t_n31, main="Autocorrelation (Run 1)")
acf(run2$t_n31, main="Autocorrelation (Run 2)")
理想情况：
- ACF 在滞后 k=5~20 步内降至接近 0（即 |ACF| < 0.1）。说明样本间自相关性低，samplefreq 设置合理。
- 如果 ACF 下降极快（如 k=1~2 步后就接近 0），可能 samplefreq 设置过大，导致 采样过于稀疏，
-
如果ACF一直不下降，如下图所示，则说明samplefreq设置过小。
ACF图
调整参数的一般原则
- burnin：根据 Trace plot 调整，通常 ≥10,000。 丢弃前面多少代的数据
- samplefreq：通过 ACF 图优化，通常 10-100。每多少代取一个样
-
nsample：确保 ESS > 200，通常 10,000-50,000。取样的样本数
burnin的值一般是总代数的20%-40%之间。
总代数=nsample * samplefreq
例如：burnin = 100,000 , samplefreq =10, nsample = 500,000 .
即总代数=500k*10=5000K, 20%即为 100K.
可以先按照我上面的示例，运行一次。然后使用Tracer打开mcmc.txt,查看Trace的曲线，可以看到前面的大概多少代不稳定，burnin值即需要大于这个代数。
image.png
上图中可以看到ESS值显然小于200，所以需要增大运行的代数和取样数。同时右边可以看到前面的100,000代都不稳定，所以下一步运行，至少需要提升burnin到100,000。
| 优化收敛中问题 | 可能原因 | 解决方案 |
|---|---|---|
| ESS 过低 | samplefreq 太小 | 增大 samplefreq 或 nsample |
| 链不收敛（R-hat > 1.1） | burnin 不足 | 增加 burnin 或调整先验 |
| 计算时间过长 | nsample 过大 | 降低 nsample，提高 samplefreq |
遇到的报错
Error: check and think about multificating trees..
解决方法把进化树从((a, b), (c, (d, e)),h);修改为(((a, b), (c, (d, e))),h);
原因是你的树不是二叉树，因为需要是二叉树，所以如果你的树的根分出来不是一个括号，就会报错，此时需要在根的外边，给其他的那一块加上括号。
Error:mcmctree species speciesA not found in master tree
解决办法是把树里面的化石分化时间和前面物种之间加上:
参考资料
https://fish-evol.org/mcmctreeExampleVert6/text1Eng.html
存在的问题
据说mcmctree对于基于蛋白的序列的树估计的分化时间不太准，所以可以考虑把蛋白转为cds序列，然后再进行分析。
成对的序列比对结果，蛋白转cds，使用ParaAT2.0
cat sample.collinearity |grep "species_prefix"|cut -f2,3 >sample.id
echo "16" >proc
ParaAT.pl -g -t -h sample.id -n cds.fa -a pep.fa -m mafft -p proc -f axt -o paraat 2>paraat.log
sample.id 是两列，同源基因的序列id的文件
cds.fa 是上面的基因的cds文件
pep.fa 是上面的基因的pep文件
比对使用的是mafft,先比对pep，然后转为对应的cds的比对。
多个蛋白序列的多重比较的结果，蛋白转cds，使用pal2nal.
wget http://www.bork.embl.de/pal2nal/distribution/pal2nal.v14.tar.gz
tar -zxvf pal2nal.v14.tar.gz
pal2nal.pl -h
pal2nal.pl -nogap -nomismatch pep.aln nuc.fa -output fasta >nuc.aln
-nogap 是移除序列中的gap
-nomismatch移除pep比对的结果和cds的结果不匹配的位点
-outout fasta 指定输出的格式是fasta格式
pep.aln是pep的多序列比对文件
nuc.fa 是cds的序列
输出的是fa格式的cds的多重比对的结果nuc.aln
参考资料
https://www.jianshu.com/p/b12e058c6597
https://www.jianshu.com/p/46b28829b078
https://yanzhongsino.github.io/2021/10/29/bioinfo_align_pep2cds/
https://www.jianshu.com/p/d61de6d47782
mcmctree的官方教程
http://abacus.gene.ucl.ac.uk/software/MCMCtree.Tutorials.pdf
https://link.springer.com/protocol/10.1007/978-1-4939-9074-0_10
https://github.com/dosreislab/mcmc3r
