MCMCtree 计算分化时间

小黑Black_007

发布于：2022年10月25日

目前我知道的可以做物种分化时间计算的软件主要有beast，beast2，r8s，paml的mcmctree，还有Nagy论文中的PhyloBayes (R package)。目前呢，我是尝试使用了beast2，非常的慢，我84个物种接近50万AA (amino acids) 每100万代运算时间差不多27h，但是一般设置是1亿代也就是需要100多天，年都过完了。暂且不提，所以开始尝试MCMCtree和r8s，本篇就是mcmctree的使用记录

1. 安装

本文主要介绍在Linux或者WSL上的安装，直接根据官网的教程走就可以了

# 下载
wget http://abacus.gene.ucl.ac.uk/software/paml4.9j.tgz
# 解压
tar xf paml4.9j.tgz
# 进入文件夹以后开始删除和编译工作
rm bin/*.exe  
cd src  
make -f Makefile  
ls -lF  
rm *.o  
mv baseml basemlg codeml pamp evolver yn00 chi2 mcmctree infinitesites ../bin  
cd ..  
ls -lF bin  
bin/baseml  
bin/codeml  
bin/evolver
echo 'PATH=$PATH:/home/dell/software/paml/bin/' >> ~/.bashrc
source ~/.bashrc

2. 文件的准备

网上常见的教程基本都在说MCMCtree用碱基序列来估算时间序列比较方便，陈连福也说了碱基序列更具有可靠性（这里我不是很理解，为什么碱基序列更具有可靠性），但是氨基酸建树的教程也有，在官方文档里边Toturial 4也是氨基酸的教程，主要参考的就是这俩。
需要准备三个文件
2.1 有根树的文件，不需要枝长，需要添加化石节点，化石节点的话是在相应的位置打个单引号，单位是MYA，也就是百万年，比如105个百万年到210个百万年，就是$’>1.05<2.10’$。

这个文件可以通过RaxML构建的树bipartition提取出来，我不知道这么提取出来没有枝长的树，所以是提取了文件以后输入：

1	sed "s/\:[0-9\.]//g" all.trees > all1.trees

2.2 序列文件，需要时phylip格式的
这一步很多都可以转，比如我有一个从钊哥那里继承的祖传代码：fas2phy.py

import re
with open('all.fa', 'r') as fin:
    sequences = [(m.group(1), ''.join(m.group(2).split()))
    for m in re.finditer(r'(?m)^>([^ \n]+)[^\n]*([^>]*)', fin.read())]
with open('all.phy', 'w') as fout:
    fout.write('%d %d\n' % (len(sequences), len(sequences[0][1])))
    for item in sequences:
        fout.write('%-20s %s\n' % item)

或者用R语言来弄：

library(devtools)
library(ape)
library(phylotools)

setwd("D:\\Agaricales5\\all\\divergencetime")
getwd()
data <- read.fasta("notall.fa")
dat2phylip(data, outfile = "notall1.phy")

然后因为有些序列之中包含了奇怪的氨基酸序列，比如B，U，这两个建议直接手动替换成X或者-，删除会缩减序列长度。因为肯定不会很多，批量删除的话容易把序列名称中的B和U改掉。
2.3 mcmctree.ctl文件
建议的是照着图中标红的区域对文件进行更改，其他的不建议动

mcmctree.ctl

        seed = -1
     seqfile = notall1.aa
    treefile = notall1.trees
    mcmcfile = mcmc.txt
     outfile = out.txt

       ndata = 1
     seqtype = 2  * 0: nucleotides; 1:codons; 2:AAs
     usedata = 3    * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
       clock = 2    * 1: global clock; 2: independent rates; 3: correlated rates
     RootAge = '<1.0'  * safe constraint on root age, used if no fossil for root.

       model = 0    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
       alpha = 0    * alpha for gamma rates at sites
       ncatG = 5    * No. categories in discrete gamma

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

     BDparas = 1 1 0.1    * birth, death, sampling
 kappa_gamma = 6 2      * gamma prior for kappa
 alpha_gamma = 1 1      * gamma prior for alpha

 rgene_gamma = 2 20 1   * gamma prior for overall rates for genes
sigma2_gamma = 1 10 1    * gamma prior for sigma^2     (for clock=2 or 3)

    finetune = 1: .1 .1 .1 .1 .1 .1  * auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr

       print = 1
      burnin = 2000
    sampfreq = 10
     nsample = 20000

3. 运行程序

基本上完成上面三个文件的准备就需要

1 2	# 最后使用 mcmctree mcmctree.ctl

在这里会发生一个报错：

这个是正常的，如下是官方教程所示：

我们需要手动的把tmp0001.ctl文件按照图中所示进行更改，删掉“out.BV”和“rst”两个文件，并且把“wag.dat“从dat文件夹中拷贝到这个文件夹里就可以执行下一步了。

1	codeml tmp0001.ctl

这样就使用WAG+Gamma生成了适当的Hessian矩阵，接下来将rst2重命名为in.BV,现在可以更改mcmctree.ctl 的 usedata = 2

回到上一目录，新建一个文件夹，将那三个文件和新生成的in.BV拷贝进去

接下来运行

1	mcmctree mcmctree.ctl

4. 补充

4.1 一条龙服务

整体来讲，按上述内容走就可以很好的运行了，但我还是觉得非常之麻烦，于是写了一个mcmctree.sh

mcmctree mcmctree.ctl;
rm out.BV rst tmp0001.ctl;
cp ../wag.dat .;
cp ../tmp0001.ctl .;
codeml tmp0001.ctl;
cp rst2 step3;
cd step3;
mv rst2 in.BV;
mcmctree mcmctree2.ctl;

就是我们在一个文件夹中放入蛋白序列文件(*.aa)、树文件(*.trees)和mcmctree.ctl文件，并创建一个step3文件夹，把这三个文件复制进去，并在这个文件夹的上一级文件夹中放入wag.dat和tmp0001.ctl就可以直接运行mcmctree.sh了。