MrBayes軟體使用教學(完整版)

閱讀完整版之前，建議先使用一次懶人版。原文請看MrBayes的wiki。

Getting Data into MrBayes

MrBayes使用的資料格式稱為Nexus，內容則可以是DNA序列、胺基酸序列、型態特徵的資料或是restriction site二態資料，以及上述四種資料的組合都可以。Nexus檔案有一個檔頭宣稱#NEXUS，以及一些基礎描述，像是有幾個taxa、特徵數量(序列長度)、檔案類型、gap用什麼符號表是等等的，然後會分成一些blocks，以begin開始，以end結束，例如這個教學所使用的primates.nex，可以用Windows的記事本軟體打開來看，就呈現下面這個樣子。

像是這個檔案就只有一個block，稱為DATA，有的則會把taxa和character分成兩個不同的block。每一個block會有三個敘述：dimensions statement, the format statement, and the matrix statement。其中dimensions statement要包含ntax，也就是taxa的數量，以及nchar，在這裡指的就是已經align好的序列長度。而format statement則有五種可能：datatype=DNA (or RNA or Protein or Standard or Restriction)。而format statement則大概會包括gap=-、missing=?、interleave=yes、match=. 這些，當然所使用的符號不一定要這些，不過要記得更改成自己使用的符號就是了，而interleave的yes或no請看以下例子（也就是所謂的"interleaved" format和"sequential" format）：

interleave=yes

Turkey    AAGCTNGGGC ATTTCAGGGT

Salmo gairAAGCCTTGGC AGTGCAGGGT

H. SapiensACCGGTTGGC CGTTCAGGGT

Chimp     AAACCCTTGC CGTTACGCTT

Gorilla   AAACCCTTGC CGGTACGCTT

GAGCCCGGGC AATACAGGGT AT

GAGCCGTGGC CGGGCACGGT AT

ACAGGTTGGC CGTTCAGGGT AA

AAACCGAGGC CGGGACACTC AT

AAACCATTGC CGGTACGCTT AA

interleave=no（interleave=no是預設值，primates.nex這個例子就是interleave=no）

Turkey    AAGCTNGGGC ATTTCAGGGT

GAGCCCGGGC AATACAGGGT AT

Salmo gairAAGCCTTGGC AGTGCAGGGT

GAGCCGTGGC CGGGCACGGT AT

H. SapiensACCGGTTGGC CGTTCAGGGT

ACAGGTTGGC CGTTCAGGGT AA

Chimp     AAACCCTTGC CGTTACGCTT

AAACCGAGGC CGGGACACTC AT

Gorilla   AAACCCTTGC CGGTACGCTT

AAACCATTGC CGGTACGCTT AA

而接在 format statement後面的則是 matrix statement，也就是align好的序列。每一條序列最前面是序列的名稱，名稱和序列之間至少有一個空格(one space)區隔。

過去Nexus檔案通常使用MacClade或是Mesquite軟體所產生，不過MrBayes version 3並不支援全部的Nexus檔案標準，所以你可能需要按照上面所述作一些調整，讓MrBayes可以接受你的Nexus檔案。此外，MrBayes對每一種資料格式只接受固定的符號，例如{A, C, G, T, R, Y, M, K, S, W, H, B, V, D, N} for DNA data, {A, C, G, U, R, Y, M, K, S, W, H, B, V, D, N} for RNA data, {A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V, X} for protein data, {0, 1} for restriction (binary) data, and {0, 1, 2, 3, 4, 5, 6, 5, 7, 8, 9} for standard (morphology) data。而如果特徵不確定，MrBayes也可以接受，例如型態特徵在2和3之間，可以用(23), (2,3), {23}, or {2,3}等方式描述，DNA也是。不過如果你的Nexus有"Equate" and "Symbols" statements ，請記得先移除，因為這些敘述MrBayes不接受。

檔案準備好以後，請記得先把檔案跟MrBayes放在同一個資料夾裡面。在MrBayes>命令提示字元後面，輸入execute primates.nex，或是更簡單一點輸入exe primates.nex。這個動作會把資料帶進程式。注意如果你要執行的檔案不是這個範例，就要改成正確的名稱，例如是Fish.txt或是Bird.nex等，檔案的檔名不能有空白。如果一切正常，應該會出現這樣的畫面：

Specifying a Model

以下所有的指令都輸入在MrBayes>的命令提示字元後面。在做分析的時候，至少有兩個指令：lset and prset，會被用來設定演化的模型（evolutionary model ）。分析之前，可以用"showmodel"的指令來檢視一下model的設定。一般來說，lset用來定義model的結構，而prset則用來定義model參數的事前機率分佈(prior probability distributions)。接下來，我們會使用GTR + I + Γ model（a General Time Reversible model with a proportion of invariable sites and a gamma-shaped distribution of rates across sites）作為粒線體序列的演化模型。如果你對分子演化的隨機模型(stochastic models)不太熟悉，我們建議你可以先看一些文章，例如Felsenstein (2004)。

在設定之前，建議先輸入help lset，這時後會跑出一堆解釋各種設定的文字，可以先忽略，或是如果你要知道每個設定的細節，也可以仔細閱讀那些文字。不過在這裡我們把注意力放在最後面的那個表格，表格的內容就是目前的設定，看起來應該會是這樣：

首先，請注意Model settings的最前面一行表示這些設定是指定給partition 1，因為在預設情況下，MrBayes會把你的資料全部分成一個partition。如果你的資料不只一個partition，請另外設定，關於partition的設定請看partitioned bayesian analyses。

表格中的第一行是Nucmodel設定，可以讓你設定鹼基取代模式(nucleotide substitution model)。如果你的資料是paired stem regions of ribosomal DNA，請設定為Doublet，例如12S或16S等基因。如果你的資料是編碼序列，可設定為Codon。而4by4設定則是一般標準的DNA取代模型，也就是只有A,T,C,G四個狀態，這也是預設值。

而接下來則是設定取代模型的結構，也就是Nst。預設值是所有鹼基的取代都有相同的速率(Nst=1)，也就是F81或JC69的model，而Nst=2則表示transitions和transversions有不同的速率，也就是K80或是HKY85的model。不過在這個教學示範裡面，我們設定為Nst=6，也就是採用GTR的model，請輸入lset nst=6。

Code的設定就看DNA是來自於哪裡，例如脊椎動物粒線體就設定為Code=Vertmt，不過如果Nucmodel不是設定為Codon，那Code的地方就不必設定。而Ploidy也不必設定(This option is used when a coalescence prior is used on trees)。不過在Rates的部份，我們就需要從預設的Equal(no rate variation across sites)改成Invgamma(gamma-shaped rate variation with a proportion of invariable sites)，請輸入lset rates=invgamma。或是我們也可以同時輸入lset nst=6 rates=invgamma一次更改兩個設定。

而Ngammacat的部份是用來設定gamma distribution的categories有幾個，預設值為4，一般而言，4是一個合理而有效的設定，因此不必更改。當然，如果增加rate categories的數量，那麼likelihood的計算換更精確，但是運算時間會增加非常的多，而對大多數的例子來說並不影響最後的結果。

而Usegibbs則是用來設定是否使用Gibbs sampler來對機率採樣，通常設定為Yes，因為這種採樣方法非常快又很節省記憶體，但是有少數問題需要注意，如果有疑慮的人可以仔細看軟體的manual。

Omegavar是用來設定nonsynonymous/synonymous rate ratio (omega)，一般設定為Equal即可。而Coding一般也是設定為all，因為各種資料都適用。而剩下的一些設定，因為我們不會使用parsimony model以及covariotide model，所以這些都維持預設值。設定完成後，再輸入help lset，就可以看到更改後的設定：

Setting the Priors

接下來我們要設定model的priors（事前參數、先驗參數），總共有六個參數：

1. the topology拓樸

2. the branch lengths分支長度

3. the four stationary frequencies of the nucleotides四種鹼基的固定頻率

4. the six different nucleotide substitution rates六種不同的鹼基取代速率

5. the proportion of invariable sites不變位點的比例

6. the shape parameter of the gamma distribution of rate variation速率變化的gamma分佈型態參數

預設的事前參數其實已經適用於大多數的分析，而我們在這個教學裡面也不會改變他們。請輸入help prset，你就可以看到在你的model裡面目前那些參數的設定是如何：

我們應該專注在以下這些參數：Revmatpr (for the six substitution rates of the GTR rate matrix-六種不同的鹼基取代速率), Statefreqpr (for the stationary nucleotide frequencies of the GTR rate matrix-四種鹼基的固定頻率), Shapepr (for the shape parameter of the gamma distribution of rate variation-速率變化的gamma分佈型態參數), Pinvarpr (for the proportion of invariable sites-不變位點的比例), Topologypr (for the topology-拓樸), and Brlenspr (for the branch lengths-分支長度). 

Revmatpr和Statefreqpr預設的事前機率密度(prior probability density)都是1.0 (flat Dirichlet)。如果我們要從資料去估計這些參數，那麼這樣的預設是合理的，因為這種設定假設對於這些參數並沒有事前知識。鹼基的頻率和速率是可以被固定的，不過通常不建議這麼作。然而偶而我們必須固定鹼基的頻率，例如使用JC或是SYM等model的時候，這時候可以輸入prset statefreqpr=fixed(equal)。

如果我們要強調鹼基頻率的相同，我們可以設定prset statefreqpr = Dirichlet(10,10,10,10)，這樣一來還是equal frequencies，但是被加強了，或是更加強設定為prset statefreqpr=Dirichlet(100,100,100,100)。Dirichlet分佈的數值總和決定了這個分佈是否被強調，而四個數字之間的平衡則決定的每個鹼基（A, C, G, T）之間的期望比例。通常在Dirichlet分佈和觀察值之間會有關連，例如Dirichlet (150,100,90,140) distribution是來自於某些reference sequences被觀察到有150 A's, 100 C's, 90 G's and 140 T's。如果你的序列和reference sequences是各自獨立的，但是reference sequences又明顯和你序列的分析有關，那麼使用這些reference sequences的Dirichlet分佈數值就可能是合理的。

在這個教學範例中，我們會維持預設參數的設定，如果你已經按照上面說的改過設定，可以輸入prset statefreqpr=Dirichlet(1,1,1,1) 或 prs st = Dir(1,1,1,1)把設定再改回來。同樣的，我們也不改變substitution rates以及其他的事前參數設定。

Checking the Model

輸入showmodel，就可以看到分析前的設定，檢查看看是否正確。

Setting up the Analysis

要進行分析的指令為mcmc，不過在分析之前，我們建議輸入help mcmc看一下設定：

所謂的seed其實指的是亂數產生器的起始值，而swapseed則是用來產生鏈(chain)轉換序列的亂數產生器起始seed。除非使用者有自己要設定的值，不然這些數值都是由系統時鐘所產生的，所以你的數值會和範例不一樣，而每次分析也都不一樣。而Ngen則是分析會跑幾個世代(generations)，也就是MCMC運算幾個cycles。建議先跑一個比較少的generations，可以先確定分析的設定都是正確的，而且也可以用來估算完整的分析大概需要花多久的時間。我們先從10000個generations開始。為了不啟動分析，我們使用mcmcp的指令來更改數值，請輸入mcmcp ngen=10000，這樣就可以設定generations為10000。你可以輸入help mcmc來確認這個設定是否已經改變。

而在預設中，MrBayes會從不同的亂數樹(random tree)跑兩個完全獨立而同時的分析，設定為Nruns=2。同時多跑一次分析對於判斷你是否從事後機率分佈取得很好的樣本是很有幫助的，所以我們建議你保留這個預設值為2。主要的概念就是從不同的亂數樹裡面開始跑每一個分析(run)，在分析的早期，兩個分析會去採樣不同的樹，但是當他們達到一致的時候（當他們從事後機率分佈產生一個好的樣本），兩顆樹的樣本(samples)就會非常相似。

為了確認MrBayes會比較來自兩個不同分析的樣本樹，確認Mcmcdiagn設定為yes，而Diagnfreq合理的值，例如每第一千個generation。MrBayes將會在每一個Diagnfreq的世代時計算許多次的診斷，並且把結果放在.mcmc的檔案裡面。而最重要的診斷，也就是測量不同分析(runs)的樣本樹相似性(similarity)也會在每一個Diagnfreq的世代顯示在螢幕上。每次當診斷被計算時，從鏈(chain)一開始，不論是一個固定數量的樣本(burnin)或是一定百分比的樣本(burninfrac)會被丟棄。而relburnin的設定則是用來決定是要採用固定的burnin值(relburnin=no)還是採用burnin百分比(relburnin=yes)。MrBayes預設為丟棄從冷鏈(cold chain)來的最開始25%樣本(relburnin=yes and burninfrac=0.25)。

在預設中，MrBayes使用Metropolis coupling來改進目標分佈的MCMC採樣，Swapfreq, Nswaps, Nchains, 和Temp同時控制Metropolis coupling的運算行為。當Nchains設定為1時，就不會加熱(no heating is used)。當Nchains設定為n且大於1時，會有n-1個熱鏈(heated chains)被使用。預設的Nchains=4，也就是MrBayes會使用3個熱鏈(heated chains)和一個冷鏈(cold chain)。在我們的經驗裡面，大於50個taxa時分析會有問題出現，此時加熱是必要的，反之對於比較小的問題通常不需要加熱就可以成功分析（In our experience, heating is essential for problems with more than about 50 taxa, whereas smaller problems often can be analyzed successfully without heating.）。因此對於大型以及難以分析的資料來說，增加超過3個熱鏈也許會有利於分析。不過增加熱鏈後，分析的時間也會隨之而等比例增加。（除非MrBayes用光了所有記憶體的空間，這時候分析會突然變得非常慢）如果你的電腦是cluster的CPU，可以使用MPI版本的MrBayes，把冷鏈和熱鏈用不同的CPU跑分析，可以大大地加快運算速度。

MrBayes使用遞增加熱(incremental heating)的方式，當chain i被他自己的事後機率提升到power 1 / (1 + iλ)，這裡的λ是Temp設定的參數所控制的溫度。如此一來加熱的效果在事後機率就會扁平化，而熱鏈就更容易在事後分佈找到隔離的峰值，也可以幫助冷鏈更快速在這些峰值之間移動。每一個Swapfreq所設定的generation，會隨機選到兩個鏈，並且企圖轉換(swap)他們的狀態。對於許多分析來說，預設的設定應該可以執行得很好。然而如果你跑遠超過三個熱鏈，你也許應該要增加轉換的數量(Nswaps)。

而Samplefreq則是用來設定鏈多常被採樣。預設值是每第100個generation就採樣一次鏈，而這樣的設定也是用於大多數的分析。然而可能是我們的分析很小，所以很快就得到收斂，因此如果更頻繁對鏈採樣也是合理的，例如每第10個generation（這樣就確保了我們在ngen設定為10000時至少可以取得1000個樣本）。要改變採樣頻率，請輸入mcmcp samplefreq=10。

當chain被採樣時，當時model參數的值會被顯示，而取代模型(substitution model)的參數可以在一個副檔名為.p的檔案中找到，這個檔案是以tab區隔的文字檔，可以被大多數的統計和繪圖軟體打開。而樹型以及分支長度的資料則在副檔名為.t的檔案裡面，這是Nexus格式的樹檔案，可以被匯入像是PAUP*或是TreeView等軟體。預設狀態下，MrBayes會把branch lengths也存在.t的檔案裡面，不過如果你跑得是很大的分析（taxa很多），而且不需要branch lengths的資訊，可以不要存branch lengths而把空間謄出來給運算使用。

設定Printfreq的參數可以控制被顯示在螢幕上的chain頻率，可以維持為預設值，也就是每第100個generation就會被顯示在螢幕上。而Startingtree則預設為每一個chain會用一顆不同的亂數樹開始，建議一般使用的話維持預設值就好。

Running the Analysis

終於，我們準備好要開始分析了。輸入mcmc，MrBayes就會開始在螢幕顯示關於model的資訊，然後列出從事後分佈採樣將會被用到的建議機制。在教學範例，看起來應該像是這樣：

請記得MrBayes將會把大多數的努力擺在改變樹型(topology, Tau)以及分支長度(branch lengths, V)的參數。在我們的經驗裡面，樹型和分支長度是最難整合的參數，因此我們讓MrBayes把大多數的運算時間花在這些參數上。如果要改變建議機率(proposal probabilities)，可以用propset指令，不過請小心不恰當的改變可能會造成的後果。在起始的log likelihood之後，MrBayes會開始顯示每第100個generation的狀態，像是這樣：

左邊第一欄是generation number，而接下來的四欄負值則是第一次分析(run)的數值，數值則是那些chains的log likelihood values，如果是中括號[]，裡面的數值就是冷鏈，而小括號()裡面的數值則是熱鏈的值。當兩種鏈成功改變狀態，他們會互換欄的位置，如果Metropolis coupling運作得很好，冷鏈應該在不同的欄位移動，這也表示冷鏈和熱鏈成功轉換(swap)狀態。如果冷鏈在欄中卡住了，熱鏈就無法成功貢獻他的狀態給冷鏈，而Metropolis coupling就無法有效運作。這時候應該作更長的分析，或是冷熱鏈之間的溫度差異應該要更低。

而欄位中間的*號區隔了不同的兩個分析(runs)，最後面的欄則顯示分析的時間，例如大約每100個generations要花一秒鐘，因為在每一個generation會有不同的移動，因此分析的時間不會一樣，而尤其是一開使分析的時後會非常不穩定。不過到後面，時間的預估會越來越準確。

When to Stop the Analysis

分析結束後，MrBayes會問你要不要繼續跑分析，在回答之前，請先看看average standard deviation of split frequencies的數值。如果兩個分析(runs)趨近於平穩分佈(stationary distribution)，我們預期分頻的平均標準差(average standard deviation of split frequencies)會接近於零，顯示兩個樣本數非常相似。在範例裡面，average standard deviation大約在1000generations以後就接近0.07，而到分析的最後，減少到只剩0.000001。你自己的資料進行分析，也許值會不一樣，不過如果如果average standard deviation的值在分析的最後階段相對來說非常的低，那麼就沒有必要繼續分析了，因此當MrBayes問"Continue with analysis? (yes/no):"時，就輸入no來結束分析吧！

雖然我們建議使用收斂診斷(convergence diagnostic)，像是standard deviation of split frequencies等來決定分析的長度，我們還是要指出其實還有更簡單但是比較不強大的方法來決定什麼時候停止分析。其實最簡單的方法就是去看冷鏈的log likelihood values，或是更精確地說，是給定參數值資料的指數機率(log probability of the data given the parameter values)。也就是在中括號裡面的數值。在分析的初期，數值通常會快速增加（絕對值減少，不過因為數值是複數，所以其實數值是增加的），這個階段稱為"burn-in"(預燒)，而通常在分析過程中，這個階段的樣本都會被捨棄。而當冷鏈的然似值(likelihood)停止增加，而開始在一個穩定的區間任意漂移的時候，分析才算達到平穩階段，而這個階段所產生的樣本才是從事後機率分佈而來的好樣本。在平穩階段，我們也期待看到不同而立的分析有近似的然似值。不過然似值的趨勢也會騙人，所以最好還是使用split frequencies的收斂，而不要只看然似值的趨勢。

當你停止分析，MrBayes會顯示出許多有助於把分析最佳化的資訊，尤其是如果你的分析很難達到收斂的狀態時。不過因為範例是很好的樣本，所以這邊就不討論那些資訊了，這個部份可以看menual的[[[FAQ 3.2|section 5]] 。

Summarizing Samples of Substitution Model Parameters

在分析過程中，取代模型參數的樣本(samples of the substitution model parameters)被寫在.p檔案裡面，看起來就像這樣：

第一個中括號裡面的數值是一個亂數產生的ID號碼，用來讓你知道這是那一個分析，而下一行則是每一欄的標題，再下去則是樣本的數值。從左到右，欄位包括了(1) the generation number (Gen)世代數; (2) the log likelihood of the cold chain (LnL)冷鏈的指數然似值; (3) the total tree length (the sum of all branch lengths, TL樹長); (4) the six GTR rate parameters (r(A<->C), r(A<->G) etc六個GTR速率的參數); (5) the four stationary nucleotide frequencies (pi(A), pi(C) etc)四種固定鹼基頻率; (6) the shape parameter of the gamma distribution of rate variation (alpha)速率變化的gamma分佈型態參數; and (7) the proportion of invariable sites (pinvar)不變位點的比例. 如果你使用的model不一樣，那麼.p檔案也會有所不同。

MrBayes提供了sump指令來對採樣的參數值作摘要。在使用這個指令以前，我們必須訂burn-in。既然我們之前使用收斂診斷來決定捨棄前面25%的樣本，並且知道在10000個generations後就達到收斂了，那麼在最前面的10000個generations捨棄25%也很合理。而我們每第十個generations採樣一次，也就是有1000個樣本(其實是1001個)，其中有250個樣本，也就是25%是要捨棄的。所以我們可以輸入指令為sump burnin=250來對.p檔案作摘要，通常這個指令會對最新產生的.p檔案執行，或是也可以指定檔案。

sump指令會先產生一個圖，橫軸是generations，縱軸是log likelihood，如果已經達到平穩狀態，那麼這個圖上的點會呈現亂數分佈的樣子(white noise)，而不應該可以看出上升或下降的趨勢。圖看起來應該是這樣：

如果你可以從你的圖中看出明顯的上升或下降，你也許應該增加你的分析長度來從事後機率分佈獲得足夠的樣本。而在sump輸出資料的最後，這樣包括有樣本參數的表格：

這個表顯示出每一個參數的平均值、變異數、95%信賴區間的最高值和最低值，以及中位數。而參數則和.p檔案的一樣：TL、r、pi、alpha、pinvar等。而最後一欄則是收斂診斷的參數：Potential Scale Reduction Factor (PSRF)，如果我們的樣本夠好，那麼PSRF值應該接近於1.0，如果你的樣本太少，那麼PSRF的值可以會顯得很分散，表示你應該增加你分析的長度。在範例裡面，看起來應該可以再跑更長一點，以獲得更好的結果。

Summarizing Samples of Trees and Branch Lengths

樹和分支長度被儲存在.t檔案裡面，這種檔案是Nexus格式，看起來像是這樣：

為了要對樹和分支長度的資訊作摘要，請輸入sumt burnin=250。因為sumt和sump是不同的指令，因此還要重新指定一次burn-in。sumt會輸出兩個樹：有可信度數值的支序樹，在樹上的每一個節點你可以看到那個部份的可能性(credibility values)，以及有著分支長度的譜系樹。而統計的摘要如下圖：

上面的樹是有可信度數值的支序樹，在樹上的每一個節點你可以看到那個部份的可能性。而下面的樹是有著分支長度的譜系樹。

sumt指令會創造出三個額外的檔案：首先是.parts檔案，包刮了taxon bipartitions的list、他們的事後機率以及分支長度。第二個檔案是.con檔案，包含了兩棵consensus樹。第一棵同時有每個clades的事後機率以及分支長度，而這顆樹可以用TreeView軟體看。第二棵樹則是只有分支長度，不過這棵樹可以被很多編輯樹的軟體編輯，像是MacClade和Mesquite。第三個產生的檔案則是.trprobs檔案，包含了在MCMC搜尋期間的所有樹，以事後機率作為排序。