最近在學(xué)習(xí)構(gòu)建進(jìn)化樹，嘗試了用PHYML構(gòu)建ML樹和用MrBayes構(gòu)建貝葉斯樹，因?yàn)闃?gòu)建ML樹時(shí)用的是網(wǎng)站的PHYML進(jìn)行的设江，操作簡便暇屋，所以就未做整理似袁，只詳細(xì)整理了一下貝葉斯樹的構(gòu)建步驟及相關(guān)的知識。算是學(xué)習(xí)筆記，有錯(cuò)誤和不足之處或其它我還沒有g(shù)et到的功能還望道友指明并告知昙衅。
一屋彪、比對序列格式轉(zhuǎn)換
將fasta格式的序列先在MEGA軟件里進(jìn)行比對，然后導(dǎo)出比對結(jié)果為PAUP format（圖1）绒尊，用記事本打開導(dǎo)出的結(jié)果畜挥，將里面的“missing=?”刪掉，將里面的"datatype=nucleotide“改為"datatype=dna”婴谱，如果是蛋白序列就是“datatype=protein”蟹但。否則在用MrBayes程序打開時(shí)會報(bào)錯(cuò)。修改好后應(yīng)是如下所示的樣子：
至此就獲得了可用的“.nexus”格式文件谭羔。注意：要將文件與MrBayes程序放在同一目錄下华糖，否則會無法打開。（nexus格式：ntax=類群數(shù)瘟裸，nchar=序列長度客叉，datatype=數(shù)據(jù)類型）
二、MrBayes軟件下載安裝
直接在其官網(wǎng)下載就可以了（http://nbisweden.github.io/MrBayes/download.html）话告，注意根據(jù)自己的電腦系統(tǒng)及配置下載安裝相應(yīng)的版本兼搏。（似乎是不用安裝的，解壓后雙擊應(yīng)用程序就可以直接用了）
三沙郭、建樹
1佛呻、運(yùn)行MrBayes，輸入“EXE filename.nex”命令病线，打開第一步中準(zhǔn)備好的“filename.nexus”格式文件吓著。執(zhí)行成功，執(zhí)行窗口會輸出文件的簡單信息送挑。
2绑莺、參數(shù)設(shè)置
? ? ? lset 設(shè)置替換模型參數(shù)
輸入“Lset nst=6 rates=invgamma”命令，該命令設(shè)置進(jìn)化模型為with gamma-distributed rate variation across sites和a proportion of invariable sites的GTR模型惕耕。模型可根據(jù)需要更改纺裁，不過一般無須更改。
? ? ? Nst：核酸替代模型赡突，1=JC69对扶，2=F81，6=GTR惭缰，可以嘗試使用三個(gè)模型運(yùn)行浪南，選擇最優(yōu)的結(jié)果
? ? ? Rates：指定序列上每個(gè)位點(diǎn)的替換速率，equal表示替換速率都是一致的漱受，gamma表示 ? ? ? ?用?gamma來確定序列上的替換速率
? ? ? prset 設(shè)置模型的先驗(yàn)信息络凿，一般不確定就設(shè)置以下默認(rèn)參數(shù)骡送，依次輸入以下命令：
“prset statefreqpr=fixed(0.2612,?0.2180,0.2563,?0.2645)”? 設(shè)置GTR模型中核苷酸平衡頻率的先驗(yàn)概率
“prset revmatpr=fixed(1.0000,?3.6695,?0.4065,0.4065,?9.2409,?1.0000)”? 設(shè)置GTR模型中替換速率的先驗(yàn)分布
“prset shapepr=fixed(0.8448)”? 設(shè)置速率分布的尺度
“prset pinvarpr=fixed(0.5703)”? 設(shè)置不變位點(diǎn)的比例
? ? ? ? mcmc 設(shè)置抽樣信息，輸入如下命令
“mcmc Ngen=10000絮记，Samplefreq=10摔踱，Printfreq=100”
? ? ? ? Ngen：設(shè)置總抽樣次數(shù)
? ? ? ? Samplefreq：設(shè)置抽樣頻率
? ? ? ? Printfreq：設(shè)置打印頻率
注意：此處需要根據(jù)界面返回的分裂頻率分支頻率的標(biāo)準(zhǔn)偏差值判斷抽樣次數(shù)是否足夠，分裂頻率分支頻率的標(biāo)準(zhǔn)偏差需要小于0.01怨愤，否則需要增加抽樣次數(shù)派敷。即：如果分裂頻率分支頻率split frequencies的標(biāo)準(zhǔn)偏差standard deviation在100,000代generations以后低于0.01，則當(dāng)程序詢問：“Continue the analysis? (yes/no)”撰洗，回答no篮愉；如果高于0.01，回答yes并輸入抽樣次數(shù)差导，直到該值低于0.01试躏。（如下圖2所示）
? ? ? ? ?總結(jié)參數(shù)結(jié)果，輸入如下命令：
?“sump burnin=250”（在此為1000個(gè)樣品设褐，即任何相當(dāng)于你取樣的25％的值）颠蕴。
程序會輸出一個(gè)關(guān)于樣品（sample）的替代模型參數(shù)的總結(jié)表，包括mean助析，mode和95 % credibility interval of each parameter犀被，要保證所有參數(shù)PSRF（the potential scale reduction factor）的值接近1.0，如果不接近貌笨，分析時(shí)間要延長弱判。
? ? ? ? 總結(jié)系統(tǒng)樹襟沮，輸入如下命令：
“sumt burnin=250”锥惋，然后就等待結(jié)果出現(xiàn)。
程序會輸出一個(gè)具有每一個(gè)分支的posterior probabilities的樹以及一個(gè)具有平均枝長mean branch lengths的樹开伏。所有輸出文件均保存在與MrBayes程序同一目錄下膀跌。生成“ . t” 文件和“ . parts” 、 “ . con” 及“ . trprobs”文件各一個(gè)固灵。其中“ . con” 文件就是含有兩棵共有樹 ( consensus trees)的文件捅伤。
參數(shù)設(shè)置知識補(bǔ)充：
1、?lset 設(shè)置的參數(shù) --------使用help?lset查看?lset設(shè)置的參數(shù)
參數(shù)Nucmodel：指核酸的類型巫玻；其可選項(xiàng)中4by4指不區(qū)分序列上的位點(diǎn)丛忆；condon指使用密碼子模型（序列上每個(gè)位點(diǎn)的替換速率會依據(jù)密碼子模型來推斷）；doublet通常用于具有協(xié)同進(jìn)化效應(yīng)的序列仍秤；一般情況下都可以使用4by4熄诡，編碼序列最好使用condon。
參數(shù)code：指密碼子編碼的規(guī)律诗力；其可選項(xiàng)中universal指通用密碼子使用規(guī)律；metmt用于推測線粒體內(nèi)的基因；葉綠體用mycoplasma冲粤。
參數(shù)ploidy：指物種是單倍體還是二倍體。
2菜拓、?prset 設(shè)置模型的先驗(yàn)信息----------使用help?prset查看相關(guān)參數(shù)及其說明
參數(shù)tratiopr：指定轉(zhuǎn)換和顛換的比例，可以使用fixed來指定也可以通過beta分布來模擬產(chǎn)生笛厦。
參數(shù)revmatpr：指定GTR模型里替換速率的先驗(yàn)分布纳鼎。
參數(shù)aamodelpr：指定氨基酸替換模型中先驗(yàn)參數(shù)的分布。
參數(shù)statefreqpr：指定GTR模型中核苷酸平衡頻率的先驗(yàn)概率裳凸。
參數(shù)shapepr：設(shè)置速率分布的尺度參數(shù)喷橙。
3、?mcmc 設(shè)置抽樣信息----------使用help?mcmc查看相關(guān)參數(shù)
參數(shù)ngen：指的是總抽樣次數(shù)登舞。
參數(shù)nruns：制定獨(dú)立分析的次數(shù)贰逾。如果為2，表明程序從兩個(gè)獨(dú)立的樹形開始抽樣菠秒，分析完成后綜合兩個(gè)分析結(jié)果疙剑。
參數(shù)nchain：設(shè)置每次分析時(shí)運(yùn)行的chain的數(shù)量。
參數(shù)samplefreq：指定從總的樣本數(shù)中抽樣的頻率践叠。這個(gè)一般和參數(shù)ngen配合使用言缤，以保證最后用于分析的樣本量足夠。舉例：samplefreq設(shè)置為100,000禁灼，ngen設(shè)置為1000管挟，這樣100，000個(gè)隨機(jī)樣本中每隔1000個(gè)抽出一個(gè)樣本弄捕，最后一共可以得到1000個(gè)樣本僻孝。
參數(shù)burninfrac：該參數(shù)控制用于分析的樣本的數(shù)量。在mcmc抽樣初期的數(shù)據(jù)往往不可靠守谓，需要去掉穿铆。burninfrac控制去掉的比例，如為0.25則表示樣本的前25%的數(shù)據(jù)被去除斋荞，因此最終用來分析的總樣本數(shù)為1000*（1-0.25）=750荞雏。
運(yùn)行之后看方差，當(dāng)方差<0.01時(shí)即可認(rèn)為抽樣達(dá)到平穩(wěn)平酿，不需要進(jìn)行更多的抽樣分析凤优，否則應(yīng)繼續(xù)增加抽樣次數(shù)以達(dá)到平穩(wěn)為止。然后在屏幕輸出結(jié)果中找到chain swap information蜈彼，當(dāng)其顯示的四條鏈之間的交換頻率在0.1~0.8之間可以認(rèn)為結(jié)果是合理的筑辨，可以進(jìn)行下一步分析，否則需要重新設(shè)置參數(shù)柳刮，包括足夠長的ngen,適當(dāng)降低temp等挖垛。
MrBayes的高級功能：
1痒钝、可以直接在序列文件后輸入相關(guān)設(shè)置參數(shù)，然后讓程序自動(dòng)運(yùn)行就可以了痢毒，但不建議將“sump”和“sumt”這兩個(gè)命令寫入送矩，因?yàn)檫@兩個(gè)命令具有診斷的作用。
2哪替、partition功能栋荸，如果所分析的序列不均一，如既有編碼區(qū)又有非編碼區(qū)凭舶，或者想把編碼區(qū)分為密碼子第一晌块、第二和第三位堿基單獨(dú)分析的話就需要使用該功能。需要在序列文件中增加如下內(nèi)容：
begin mrbayes;
charset pos1=1./3;
charset pos1=1./3;
charset pos1=1./3;
partition region=3:pos1,pos2,pos3;
Iset applyto=(all) nst=6 rates=gamma;
prset applyto=(all) ratepr=variable;
end;
charset用于設(shè)置變量并賦值帅霜；1./3指從第一個(gè)位開始匆背，每隔三個(gè)位點(diǎn)取出一個(gè)值，并把這些值用變量pos1表示身冀，代表密碼子的第一位钝尸。其它類推。partition一行是告訴程序序列分為三部分搂根，prset一行是用來指定三個(gè)部分的參數(shù)是獨(dú)立估計(jì)的珍促。
如果序列分為編碼區(qū)和非編碼區(qū)可以按一下寫法：
charset noncoding=1-200;
charset coding=200-600;
partition region=2:?noncoding,?coding;
prset?ratepr=variable;
3、指定外群剩愧，在一組序列中可以指定外群猪叙，如果不指定，則以序列文件的第一個(gè)物種作為外群仁卷。
指定外群設(shè)置命令為：outgroup 5或者outgroup my_taxon 5（數(shù)字指要指定物種在序列中的位置）穴翩。
