目前我知道的可以做物種分化時間計算的軟件主要有beast，beast2，r8s，paml的mcmctree，還有Nagy論文中的PhyloBayes (R package)。目前呢，我是嘗試使用了beast2，非常的慢，我84個物種接近50萬AA (amino acids) 每100萬代運算時間差不多27h，但是一般設置是1億代也就是需要100多天，年都過完了。暫且不提，所以開始嘗試MCMCtree和r8s，本篇就是mcmctree的使用記錄。
1. 安裝
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 有根樹的文件，不需要枝長，需要添加化石節(jié)點，化石節(jié)點的話是在相應的位置打個單引號，單位是MYA，也就是百萬年，比如105個百萬年到210個百萬年，就是'>1.05<2.10'。
[站外圖片上傳中...(image-b5b0e1-1666680122635)]
這個文件可以通過RaxML構建的樹bipartition提取出來，我不知道這么提取出來沒有枝長的樹，所以是提取了文件以后輸入：
sed "s/\:[0-9\.]//g" all.trees > all1.trees
2.2 序列文件，需要時phylip格式的
這一步很多都可以轉，比如我有一個從釗哥那里繼承的祖?zhèn)鞔a：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文件
建議的是照著圖中標紅的區(qū)域對文件進行更改，其他的不建議動
[圖片上傳失敗...(image-1a8a-1666680122635)]
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. 運行程序
基本上完成上面三個文件的準備就需要
# 最后使用
mcmctree mcmctree.ctl
在這里會發(fā)生一個報錯：
[圖片上傳失敗...(image-18a00d-1666680122635)]
這個是正常的，如下是官方教程所示：
[站外圖片上傳中...(image-50eb32-1666680122635)]
我們需要手動的把tmp0001.ctl文件按照圖中所示進行更改，刪掉“out.BV”和“rst”兩個文件，并且把“wag.dat“從dat文件夾中拷貝到這個文件夾里就可以執(zhí)行下一步了。
codeml tmp0001.ctl
這樣就使用WAG+Gamma生成了適當的Hessian矩陣，接下來將rst2重命名為in.BV,現在可以更改mcmctree.ctl 的 usedata = 2
回到上一目錄，新建一個文件夾，將 那三個文件和新生成的in.BV拷貝進去
接下來運行
mcmctree mcmctree.ctl
參考資料：
PAML的官方網站：
Phylogenetic Analysis by Maximum Likelihood (PAML) (ucl.ac.uk)
這一篇介紹了原理：
估算系統(tǒng)樹分歧時間 —— paml.mcmctree,r8s | 生信技工
(yanzhongsino.github.io)
這幾篇是介紹了MCMC估算的pepline：
mcmctree 定年 —— 使用氨基酸序列 (gitee.io)
mcmctree估算物種分歧時間 - 簡書 (jianshu.com)
使用PAML進行分歧時間計算 | 陳連福的生信博客 (chenlianfu.com)
這一篇介紹了很多參數
PAML軟件中的mcmctree命令估算物種分歧時間 - BPSO_mynotes - 博客園 (cnblogs.com)
這一篇介紹的比較全面，但是都沒講清楚：
【比較基因組】從進化樹到分化時間 - 簡書 (jianshu.com)
這一篇是介紹r8s的：
使用 r8s 估算物種分歧時間 | 陳連福的生信博客 (chenlianfu.com)
