寫在前面：我對建樹也是一知半解，這里只是想記錄一下自己跟別人學(xué)習(xí)的建樹方法谚咬，可能不具有普適性横殴。但畢竟寫在公眾平臺，大家選擇性參考茬缩。
主要的步驟包括用jModelTest來選擇核苷酸替代模型赤惊，用phylosuite進(jìn)行.nex文件的準(zhǔn)備，用在線建樹網(wǎng)站CIPRES（http://www.phylo.org/）進(jìn)行貝葉斯方法建樹凰锡。
1 jModelTest來選擇核苷酸替代模型
這一步具體怎么做參考我寫過的這篇文章：http://www.reibang.com/p/e5cfad89a1a4
2 文件的準(zhǔn)備
在這里我用自己的已經(jīng)進(jìn)行過多重比對（跑過mafft）的序列文件（.fasta）開始進(jìn)行建樹未舟。
用Phylosuite準(zhǔn)備用于貝葉斯法建樹的.nex文件圈暗。
這里簡單的說一下，貝葉斯法建樹要選擇特定的模型來建樹裕膀，首先一般會用jModeltest等軟件選出最合適的模型和各項參數(shù)员串。然后把要用的模型、參數(shù)信息和序列匯總在一起昼扛，生成一個.nex文件寸齐，提交到在線建樹網(wǎng)站CIPRES。
phylosuite是一個流程化的建樹軟件抄谐，從序列的下載到建樹渺鹦，整套流程都能在上面實現(xiàn)，具體的下載蛹含、安裝和使用毅厚，可以去官網(wǎng)看介紹：http://phylosuite.jushengwu.com/
打開phylosuite，點擊上方工具欄File-Import files or IDs浦箱，出現(xiàn)如下的對話框吸耿，Choose Files，選擇要導(dǎo)入的fasta文件憎茂。
導(dǎo)入后珍语，點擊工具欄Phylogeny-ModelFinder，出現(xiàn)如下的對話框：
基本的設(shè)定填一下竖幔，就可以按Start了板乙。
跑完后，會彈出一個對話框拳氢，關(guān)掉就可以了募逞。
點擊上方工具欄Phylogeny-MrBayes，會出現(xiàn)下面的對話框：
點擊Show MrBaves Data Block馋评，出現(xiàn)下面的對話框放接，
點擊Save to File,保存在你指定的位置，這就是我們需要的.nex文件留特。
然后纠脾，我們需要用記事本或其他文本編輯軟件打開該.nex文件，根絕jModelTest得到的結(jié)果進(jìn)行修改蜕青。
.nex文件的內(nèi)容主要分為三部分：
第一部分:
#NEXUS
BEGIN DATA;
dimensions ntax=86 nchar=12275;
format missing=?
datatype=DNA gap= - interleave;
這里的dimensions ntax是分類單元的數(shù)量苟蹈，就是一共對比了多少條序列；nchar是比對矩陣的總長度右核。
這一部分不需要進(jìn)行改動慧脱。
第二部分：
這一部分就是比對矩陣的信息，也不需要修改贺喝。
第三部分：
在用Phylosuite直接導(dǎo)出的.nex文件中菱鸥，第三部分是這樣的：
END;
begin mrbayes;
log start filename = log.txt;
lset nst=6 rates=invgamma Ngammacat=4;
prset statefreqpr = fixed(empirical);
mcmcp ngen=2000000 printfreq=1000 samplefreq=100 nchains=4 nruns=2 savebrlens=yes checkpoint=yes checkfreq=5000;
mcmc;
sumt conformat=Simple contype=Halfcompat relburnin=yes burninfrac=0.25;
sump relburnin=yes burninfrac=0.25;
end;
這里我們要根據(jù)jModelTest的結(jié)果進(jìn)行修改的部分是：
lset nst=6 rates=invgamma Ngammacat=4;
prset statefreqpr = fixed(empirical);
jModelTest根據(jù) BIC標(biāo)準(zhǔn)得到的進(jìn)化模型信息為：
這里計算出的模型是TPM2uf+I+G宗兼。
第一行中，lset nst 定義位點替換的模型（實際可以理解為 AC氮采、AT殷绍、AG、GC鹊漠、GT篡帕、TC 每對堿基的替換），用不同數(shù)字代替不同模型：1 代表所有替換率相同（如模型 JC69 或 F81）贸呢； 2 代表轉(zhuǎn)換和顛換可以有不同的替換率（如模型 K80 或 HKY85）； 6 代表所有替換率都不同（如模型 GTR）拢军。
常見的模型對應(yīng)的nst如下：
Jukes-Cantor (JC, nst=1): equal base frequencies, all substitutions equally likely (PAUP* rate classification: aaaaaa, PAML: aaaaaa) (Jukes and Cantor 1969)
Felsenstein 1981 (F81, nst=1): variable base frequencies, all substitutions equally likely (PAUP*: aaaaaa, PAML: aaaaaa) (Felsenstein 1981)
Kimura 2-parameter (K80, nst=2): equal base frequencies, one transition rate and one transversion rate (PAUP*: abaaba, PAML: abbbba) (Kimura 1980)
Hasegawa-Kishino-Yano (HKY, nst=2): variable base frequencies, one transition rate and one transversion rate (PAUP*: abaaba, PAML: abbbba) (Hasegawa et. al. 1985)
Tamura-Nei (TrN): variable base frequencies, equal transversion rates, variable transition rates (PAUP*: abaaea, PAML: abbbbf) (Tamura Nei 1993)
Kimura 3-parameter (K3P): variable base frequencies, equal transition rates, two transversion rates (PAUP*: abccba, PAML: abccba) (Kimura 1981)
transition model (TIM): variable base frequencies, variable transition rates, two transversion rates (PAUP*: abccea, PAML: abccbe)
transversion model (TVM): variable base frequencies, variable transversion rates, transition rates equal (PAUP*: abcdbe, PAML: abcdea)
symmetrical model (SYM): equal base frequencies, symmetrical substitution matrix (A to T = T to A) (PAUP*: abcdef, PAML: abcdef) (Zharkikh 1994)
general time reversible (GTR, nst=6): variable base frequencies, symmetrical substitution matrix (PAUP*: abcdef, PAML: abcdef) (e.g., Lanave et al. 1984, Tavare 1986, Rodriguez et. al. 1990)
注：K81uf+G = TPM1uf
TrNef is actually SYM
The TIM1, TIM2, TIM3, TPM1uf, TPM2uf, TPM3uf and TrN substitution models were replaced by the GTR model(nst=6)
rates 定義位點之間的替換率楞陷，有以下幾種選擇：
equal：位點的替換率無差異
?gamma：位點的替換率呈 gamma 分布，對應(yīng)+G
adgamma：位點的替換率自相關(guān)茉唉，邊緣位點替換率呈 gamma 分布固蛾，相鄰位點有相關(guān)的替換率
propinv：一定比例位點的替換率是恒定的，對應(yīng)+I
invgamma：一定比例位點的替換率是恒定的度陆，剩下位點的替換率呈 gamma 分布艾凯，對應(yīng)+I+G。
我們這里的TPM2uf可以用GTR來替代懂傀，也就是lset nst=6趾诗；I+G就是invagamma，也就是rates=invgamma蹬蚁。
第二行中的 prset statefreqpr = fixed(empirical)用的是經(jīng)驗值恃泪，我們要根據(jù)jModelTest的結(jié)果進(jìn)行修改：
首先：prset statefreqpr = fixed(0.3928,0.0901,0.0754,0.4416)
括號內(nèi)的四個數(shù)值就是freqA、freqC犀斋、freqG和freqT贝乎。
其次，shapepr 為shape parameter of gamma distribution of rate variation 叽粹，即gamma值览效，修改為shapepr = fixed(1.1630)，括號內(nèi)的為gamma值虫几；
pinvar為p-inv的值锤灿，即 pinvar = fixed(0.3690)
最后：revmat = fixed(0.7601,8.4971,0.7601,1.0000,8.4971,1.0000)，括號內(nèi)為R(a)持钉、R(b)衡招、R(c)、R(d)每强、R(e)和R(f)值始腾。
因此州刽，我們根據(jù)jModelTest結(jié)果修改后的是：
lset nst=6 rates=invgamma Ngammacat=4;
prset statefreqpr = fixed(0.3928,0.0901,0.0754,0.4416) pinvar = fixed(0.3690) shapepr = fixed(1.1630) revmat = fixed(0.7601,8.4971,0.7601,1.0000,8.4971,1.0000);
下面再來解釋一下進(jìn)行馬可夫鏈蒙特卡洛估算（Markov chain Monte Carlo simulation）的參數(shù)：
mcmc ngen=20000
samplefreq=100
diagnfreq=1000
以上舉例的參數(shù)設(shè)置，會保證你從后驗分布中得到至少200（mcmc ngen/samplefreq=20002/100）個取樣浪箭，diagnfreq指每隔多少代進(jìn)行一次診斷計算穗椅，在這里就是每隔1000代進(jìn)行一次。對于越大的數(shù)據(jù)奶栖，一般設(shè)定進(jìn)行運算分析的時間更長匹表，采樣頻率更小。默認(rèn)的采樣頻率（samplefreq）是500宣鄙，默認(rèn)的診斷頻率（diagnfreq）是5000袍镀，默認(rèn)的運行長度（mcmc ngen）是1,000,000。
3 在線網(wǎng)站建樹
我們要將.nex文件上傳到在線建樹網(wǎng)站CIPRES（http://www.phylo.org/）上用Mrbayes建樹冻晤。
關(guān)于具體的一些操作在之前說最大似然法建樹時提到過苇羡，大家可以參考一下：http://www.reibang.com/p/cdd3b3adc16f，也可以自己試一下鼻弧。
在Select Tool選擇如下：
在Set Parameters设江，要設(shè)定的參數(shù)如下：
注意：要勾選My Data Contains a MrBayes Data Block,就是說我們的.nex文件里已經(jīng)包括了很多參數(shù)信息。不需要另外設(shè)定了攘轩。Maximum Hours to Run那里的數(shù)值是根據(jù)自己的賬號還有的時間設(shè)定的叉存，因為這個網(wǎng)站的規(guī)定是每個賬號每年有限額的跑數(shù)據(jù)時間，這個數(shù)值要在你剩余的時間范圍內(nèi)度帮。
別的參數(shù)默認(rèn)即可（或者說我默認(rèn)了）歼捏。
跑完之后，在結(jié)果里笨篷，要下載的樹文件如下圖：
Download之后甫菠，用Figtree等軟件打開，就可以看結(jié)果了冕屯。
4 linux系統(tǒng)中MrBayes建樹
在linu新系統(tǒng)中用Mrbayes建樹也很簡單
在安裝了miniconda的前提下寂诱，用conda安裝MrBayes：
conda install -c bioconda mrbayes
然后一行命令就可以：
mpirun -np 8 mb -i /lustre/home/liurong/test11/81matk.nex
mpirun -np ：核心數(shù)目
mb ：.nex文件的絕對路徑
5 如何評價Mrbayes運算的結(jié)果
運算結(jié)束后的輸出文件有下面這些：
在輸出文件里面，查看名為log_mt_CDS.txt的log文件安聘，這個文件名是自己在輸入文件中通過log start filename = log_mt_CDS.txt設(shè)置的痰洒，自己想賦予log文件什么名字都可以。
如上圖所示：如果運算結(jié)束時的最后一個Average standard deviation of split frequences 值小于0.01浴韭，就代表結(jié)果是收斂的丘喻，是可以的，如果值太大念颈，需要調(diào)整參數(shù)泉粉，把輸入文件中的“mcmc ngen=20000”設(shè)置的更大，直到Average standard deviation of split frequences 值小于0.01。如果你關(guān)注的分支在系統(tǒng)發(fā)育樹中的支持率較高嗡靡，那么該值小于0.05也是可以接受的跺撼。
上圖展示的表格是對替換模型參數(shù)的采樣的匯總，
替換模型的所有參數(shù)的PSRF （potential scale reduction factor )值應(yīng)該接近1讨彼，代表運算收斂歉井。ESS（Estimated Sample Size）值應(yīng)該小于100的話代表對應(yīng)參數(shù)的采樣不足，如果參數(shù)的PSRF不接近1或ESS值小于100哈误，都應(yīng)該增加運算時間哩至。
分支和節(jié)點模型的PSRF （potential scale reduction factor )值應(yīng)該接近1，代表運算收斂蜜自，否則增加運算時間菩貌。
