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/ 查询两个物种间的分化时间
image.png
比如上图搜索玉米和拟南芥的分开时间，可以看到结果是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
输出结果解析
FigTree.tre 生成含有分歧时间的超度量树文件
使用figtree打开生成的文件FigTree.tre，即可直接看到进化树，在左侧工具栏勾选对应的栏
Figtree的界面
绘制出来的图的下面的标尺的时间顺序不对是反的。使用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()
mcmctreeR绘制出的图
遇到的报错
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
网友评论