mcmctree需要3個(gè)輸入文件
- 多序列比對的phy格式文件（注意：這個(gè)phy和一般軟件的phy不太一樣）使用腳本把多序列比對的fa格式轉(zhuǎn)為phy格式
- 準(zhǔn)備包含2個(gè)物種之間的化石時(shí)間的進(jìn)化樹
- mcmctree.ctl運(yùn)行文件
準(zhǔn)備phy格式文件supergene.phy
如果現(xiàn)有多序列比對文件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就是本次需要的文件
如果現(xiàn)有的是多序列比對是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
準(zhǔn)備進(jìn)化樹文件mcmc.input.tree
從單拷貝基因構(gòu)建的進(jìn)化樹开呐，可以手動(dòng)刪除除了物種名之外的信息。也可以使用下面的R腳本刪除
library(ape)
MLtree <- read.tree("all.fa.Jul12.contree") #讀取單拷貝進(jìn)化樹
MLtree$edge.length <- NULL #刪除分支長度信息
MLtree$node.label <- NULL
write.tree(MLtree,"MLtree.nobranch.tree") #輸出只有物種名的進(jìn)化樹
在手動(dòng)在第一行添加上物種的數(shù)量25和1中間是一個(gè)空格确买。
最終的進(jìn)化樹文件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));
中間的時(shí)間格式是百萬年喉脖。這個(gè)需要從http://www.timetree.org/ 查詢兩個(gè)物種間的分化時(shí)間
比如上圖搜索玉米和擬南芥的分開時(shí)間，可以看到結(jié)果是142.1和163.5MYA.
如果是進(jìn)化樹中祖驱，則寫為
>142.1<163.5书闸。此處的單位是MYA,百萬年紊选》却罚看到兩者是單子葉和雙子葉的分開時(shí)間姥敛。
準(zhǔn)備運(yùn)行的配置文件mcmctree.ctl
seed = -1 # 表示使用系統(tǒng)當(dāng)前時(shí)間為隨機(jī)數(shù)
seqfile = supergene.phy # 多序列比對文件
treefile = mcmc.input.tree # 帶有校準(zhǔn)信息有根樹
mcmcfile = mcmc.txt #輸出mcmc信息文本文件，可以用Tracer軟件查看
outfile = out.txt # 輸出文件瞎暑，記錄了超度量樹和進(jìn)化速率等參數(shù)信息
ndata = 1 # 輸入的多序列比對的數(shù)據(jù)個(gè)數(shù)彤敛，這是密碼子位置的數(shù)據(jù)；如果有一個(gè)了赌，則設(shè)置為1
seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs #數(shù)據(jù)類型
usedata = 3 * 0: no data; 1:seq like; 2:normal approximation; 3:out.BV (in.BV) # 設(shè)置是否利用多序列比對的數(shù)據(jù)：\
#0墨榄，表示不使用多序列比對數(shù)據(jù)，則不會(huì)進(jìn)行l(wèi)ikelihood計(jì)算勿她，雖然能得到mcmc樹且計(jì)算速度飛快袄秩，但是其分歧時(shí)間結(jié)果是有問題的；\
#1逢并，表示使用多序列比對數(shù)據(jù)進(jìn)行l(wèi)ikelihood計(jì)算之剧，正常進(jìn)行MCMC，是一般使用的參數(shù); \
#2砍聊，進(jìn)行正常的approximation likelihood分析猪狈，此時(shí)不需要讀取多序列比對數(shù)據(jù)，直接讀取當(dāng)前目錄中的in.BV文件辩恼。該文件是使用usedata = 3參數(shù)生成的out.BV文件重命名而來的。\
#此外谓形，由于程序BUG灶伊，當(dāng)設(shè)置usedata = 2時(shí)，一定要在改行參數(shù)后加 *寒跳，否則程序報(bào)錯(cuò) Error: file name empty.. \
#3聘萨，程序利用多序列比對數(shù)據(jù)調(diào)用baseml/codeml命令對數(shù)據(jù)進(jìn)行分析，生成out.BV文件童太。由于mcmctree調(diào)用baseml/codeml進(jìn)行計(jì)算的參數(shù)設(shè)置可能不太好（特別是對蛋白序列進(jìn)行計(jì)算時(shí)）米辐，\
#推薦自己修該軟件自動(dòng)生成的baseml/codeml配置文件胸完，然后再手動(dòng)運(yùn)行baseml/codeml命令，再整合其結(jié)果文件為out.BV文件翘贮。
clock = 1 * 1: global clock; 2: independent rates; 3: correlated rates # (1) 全局分子鐘,適用于近緣物種（兩物種分化時(shí)間小于20個(gè)million分化時(shí)間或者序列差異小于5%即近緣物種）\
# (2) 樹枝上進(jìn)化速率服從獨(dú)立同分布的對數(shù)正態(tài)分布（推薦使用）
# (3) 樹枝上的進(jìn)化速率自相關(guān)赊窥。
RootAge = '<54' * safe constraint on root age, used if no fossil for root. # 設(shè)置root節(jié)點(diǎn)的分歧時(shí)間，一般設(shè)置一個(gè)最大值狸页。
model = 0 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
# 堿基替換模型：
#（0）無參數(shù)（適用于近緣物種）
#（1）參數(shù)為轉(zhuǎn)換顛換比Kappa
#（2）參數(shù)為T锨能，C，A芍耘，G的頻率
#（3）最復(fù)雜的進(jìn)化模型
#（4）參數(shù)為T址遇，C，A斋竞，G的頻率+kappa（適用于遠(yuǎn)緣物種）
alpha = 0.5 * alpha for gamma rates at sites
# 核酸序列中不同位點(diǎn)倔约，其進(jìn)化速率不一致，其變異速率服從GAMMA分布坝初。一般設(shè)置GAMMA分布的alpha值為0.5浸剩。\
# 若該參數(shù)值設(shè)置為0，則表示所有位點(diǎn)的進(jìn)化速率一致脖卖。此外乒省，當(dāng)userdata = 2時(shí)，alpha畦木、ncatG袖扛、alpha_gamma、kappa_gamma這4個(gè)GAMMA參數(shù)無效十籍。因?yàn)閡serdata = 2時(shí)蛆封，不會(huì)利用到多序列比對的數(shù)據(jù)。
ncatG = 5 * No. categories in discrete gamma # 設(shè)置離散型GAMMA分布的categories值勾栗。
cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)?
BDparas = 1 1 0.1 * birth, death, sampling # 設(shè)置出生率惨篱、死亡率和取樣比例。若輸入有根樹文件中的時(shí)間單位發(fā)生改變围俘，則需要相應(yīng)修改出生率和死亡率的值砸讳。例如，時(shí)間單位由100Myr變換為1Myr界牡，則要設(shè)置成".01 .01 0.1"簿寂。
kappa_gamma = 6 2 * gamma prior for kappa #設(shè)置kappa（轉(zhuǎn)換/顛換比率）的GAMMA分布參數(shù)
alpha_gamma = 1 1 * gamma prior for alpha #設(shè)置GAMMA形狀參數(shù)alpha的GAMMA分布參數(shù)
rgene_gamma = 2 20 1 * gammaDir prior for rate for genes #設(shè)置序列中所所有位點(diǎn)平均[堿基/密碼子/氨基酸]替換率的Dirichlet-GAMMA分布參數(shù)：
sigma2_gamma = 1 10 1 * gammaDir prior for sigma^2 (for clock=2 or 3) #設(shè)置所有位點(diǎn)進(jìn)化速率取對數(shù)后方差（sigma的平方）的Dirichlet-GAMMA分布參數(shù)：alpha=1、beta=10宿亡、初始方差值為1常遂。
# 當(dāng)clock參數(shù)值為1時(shí)，表示使用全局的進(jìn)化速率挽荠，各分枝的進(jìn)化速率沒有差異克胳，即方差為0平绩，該參數(shù)無效；
# 當(dāng)clock參數(shù)值為2時(shí)漠另，若修改了時(shí)間單位捏雌，該參數(shù)值不需要改變；
# 當(dāng)clock參數(shù)值為3時(shí)酗钞，若修改了時(shí)間單位腹忽，該參數(shù)值需要改變。
finetune = 1: .1 .1 .1 .1 .1 .1 * auto (0 or 1): times, musigma2, rates, mixing, paras, FossilErr #冒號前的值設(shè)置是否自動(dòng)進(jìn)行finetune砚作，一般設(shè)置成1窘奏，然后程序自動(dòng)進(jìn)行優(yōu)化分析
print = 1 * 0: no mcmc sample; 1: everything except branch rates 2: everything # 設(shè)置打印mcmc的取樣信息
#0，不打印mcmc結(jié)果葫录；
#1着裹，打印除了分支進(jìn)化速率的其它信息（各內(nèi)部節(jié)點(diǎn)分歧時(shí)間、平均進(jìn)化速率米同、sigma2值）骇扇；
#2，打印所有信息面粮。
burnin = 2000 #將前2000次迭代burnin后少孝，再進(jìn)行取樣（即打印出該次迭代計(jì)算的結(jié)果信息，各內(nèi)部節(jié)點(diǎn)分歧時(shí)間熬苍、平均進(jìn)化速率稍走、sigma2值和各分支進(jìn)化速率等）
sampfreq = 10 #每10次迭代則取樣一次。
nsample = 20000 # 當(dāng)取樣次數(shù)達(dá)到該次數(shù)時(shí)柴底，則取樣結(jié)束（程序也將運(yùn)行結(jié)束）婿脸。
注意一定要修改的幾個(gè)參數(shù)：
seqfile = supergene.phy # 多序列比對文件
treefile = mcmc.input.tree # 帶有校準(zhǔn)信息有根樹
ndata = 1 # 輸入的多序列比對的數(shù)據(jù)個(gè)數(shù)，這是密碼子位置的數(shù)據(jù)柄驻；如果有一個(gè)狐树，則設(shè)置為1,指的是你的phy文件有多少組密碼子的位置，我的是只有1個(gè)鸿脓。一般用單拷貝的都是1個(gè)抑钟。
seqtype = 2 * 0: nucleotides; 1:codons; 2:AAs #數(shù)據(jù)類型 因?yàn)槲矣玫氖堑鞍仔蛄斜葘Φ模源颂幨?野哭，根據(jù)你建樹的實(shí)際情況選擇
usedata = 3 此處設(shè)置為3味赃，用于獲取輸出的out.BV,用于下一步分析。
RootAge = '<0.54' 這個(gè)一定要設(shè)置一個(gè)最大的值虐拓，就是你的進(jìn)化樹的根的物種和其他物種的分開時(shí)間，如果這個(gè)值設(shè)置的比里面物種的分開時(shí)間還小傲武，則會(huì)導(dǎo)致錯(cuò)誤蓉驹。注意單位
開始第一步分析
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
#再看一下是不是修改成功了城榛，如果失敗了，需要手動(dòng)修改
grep usedata mcmctree2.ctl
mcmctree mcmctree2.ctl
輸出結(jié)果解析
FigTree.tre 生成含有分歧時(shí)間的超度量樹文件
使用figtree打開生成的文件FigTree.tre态兴，即可直接看到進(jìn)化樹狠持，在左側(cè)工具欄勾選對應(yīng)的欄
繪制出來的圖的下面的標(biāo)尺的時(shí)間順序不對是反的。使用reverse后就都是負(fù)數(shù)瞻润。
還可以使用mcmctreeR這個(gè)R包來進(jìn)行可視化
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()
遇到的報(bào)錯(cuò)
Error: check and think about multificating trees..
解決方法把進(jìn)化樹從((a, b), (c, (d, e)),h);修改為(((a, b), (c, (d, e))),h);
原因是你的樹不是二叉樹喘垂，因?yàn)樾枰嵌鏄洌匀绻愕臉涞母殖鰜聿皇且粋€(gè)括號绍撞，就會(huì)報(bào)錯(cuò)正勒，此時(shí)需要在根的外邊，給其他的那一塊加上括號傻铣。
Error:mcmctree species speciesA not found in master tree
解決辦法是把樹里面的化石分化時(shí)間和前面物種之間加上:
參考資料
https://fish-evol.org/mcmctreeExampleVert6/text1Eng.html
存在的問題
據(jù)說mcmctree對于基于蛋白的序列的樹估計(jì)的分化時(shí)間不太準(zhǔn)章贞，所以可以考慮把蛋白轉(zhuǎn)為cds序列，然后再進(jìn)行分析非洲。
成對的序列比對結(jié)果鸭限，蛋白轉(zhuǎn)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败京，然后轉(zhuǎn)為對應(yīng)的cds的比對。
多個(gè)蛋白序列的多重比較的結(jié)果梦染，蛋白轉(zhuǎn)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比對的結(jié)果和cds的結(jié)果不匹配的位點(diǎn)
-outout fasta 指定輸出的格式是fasta格式
pep.aln是pep的多序列比對文件
nuc.fa 是cds的序列
輸出的是fa格式的cds的多重比對的結(jié)果nuc.aln
參考資料
http://www.reibang.com/p/b12e058c6597
http://www.reibang.com/p/46b28829b078
https://yanzhongsino.github.io/2021/10/29/bioinfo_align_pep2cds/
http://www.reibang.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
