劉小澤寫于18.11.21
聚類樹(shù)是個(gè)有意思的東西，可以十分直觀地看到每個(gè)物種的親緣關(guān)系
之前沒(méi)有涉及到相關(guān)的分析，相關(guān)的工具只是聽(tīng)說(shuō)過(guò)大家都在用的MEGA，另外還有Y叔的ggtree。很好奇但是一直沒(méi)時(shí)間學(xué)，但是昨天被問(wèn)到一個(gè)關(guān)于Bayes建樹(shù)的問(wèn)題，于是我想試試。這次先簡(jiǎn)單介紹下貝葉斯方法，從難入易循序漸進(jìn)
聚類樹(shù)什么用？
我們經(jīng)常聽(tīng)到：系統(tǒng)發(fā)育樹(shù)、系統(tǒng)進(jìn)化樹(shù)、系統(tǒng)發(fā)生樹(shù)，其實(shí)都是為了推測(cè)物種進(jìn)化機(jī)制、預(yù)測(cè)機(jī)制背后的關(guān)鍵作用位置
概率名稱
- 先驗(yàn)概率：基于背景知識(shí)、歷史數(shù)據(jù)等經(jīng)驗(yàn)，在事件發(fā)生前預(yù)先判斷的概率
- 后驗(yàn)概率（posterior probability）：貝葉斯統(tǒng)計(jì)中，一個(gè)隨機(jī)事件或者一個(gè)不確定事件的后驗(yàn)概率是在考慮和給出相關(guān)證據(jù)或數(shù)據(jù)后所得到的條件概率（來(lái)自：維基百科）它是基于先驗(yàn)概率求得的反向條件概率
建樹(shù)方法
基于距離的方法包括：UPGMA（現(xiàn)在很少用了）、ME（Minimum Evolution 最小進(jìn)化法）、NJ（Neighbor-Joining 鄰接法）
基于特征：MP（Maximum parsimony 最大簡(jiǎn)約法）、ML(Maximum likelihood 最大似然法)、BI（Bayesian Inference貝葉斯法）
建樹(shù)完成，一般需要用bootstrap（自展率）或者Posterior probability（后驗(yàn)概率）來(lái)評(píng)估
貝葉斯建樹(shù)
這種方法是最準(zhǔn)確，但同時(shí)也最慢【貌似這是軟件的通病，比如比對(duì)軟件STAR也是如此】
軟件安裝
軟件的官網(wǎng)在：https://nbisweden.github.io/MrBayes/download.html
提供了pdf的官方教程【后臺(tái)回復(fù)“mb”獲取】
軟件安裝可以選擇三大系統(tǒng)平臺(tái)，但是我還是比較推薦linux，能用服務(wù)器就用服務(wù)器運(yùn)行，一是節(jié)省時(shí)間，二是減少自己電腦消耗。
安裝也很簡(jiǎn)單，直接用conda安裝就好，為了避免軟件安裝過(guò)程中出現(xiàn)一些潛在的沖突問(wèn)題，可以新建一個(gè)conda 環(huán)境
conda create -n mrbayes
source activate mrbayes
conda install -y mrbayes
快速開(kāi)始
快速教程分為四部分：讀取Nexus文件、設(shè)置進(jìn)化模型、運(yùn)行程序、歸納樣本
-
第一部分：先下載好一個(gè)測(cè)試文件
mkdir ~/mrbayes && cd ~/mrbayes wget http://ib.berkeley.edu/courses/ib200b/labs/primates.nex然后輸入
mb就會(huì)運(yùn)行mrbayes命令【需要注意的是，每一個(gè)mb的命令都是以MrBayes >開(kāi)頭】使用命令
execute primates.nex就把數(shù)據(jù)讀進(jìn)來(lái)（如果你的文件不在當(dāng)前操作目錄，最好使用文件全路徑表示）
-
第二部分：在
MrBayes >后面，輸入lset nst=6 rates=invgamma這是將進(jìn)化模型設(shè)置成GTR（General Time Reversible）模型，包含了變異位點(diǎn)的gamma分布率和非變異位點(diǎn)的比率【只是針對(duì)此處的演示DNA數(shù)據(jù)】
如果數(shù)據(jù)不是DNA或者RNA，或者想調(diào)用（invoke）其他模型，又或者不想用默認(rèn)參數(shù)，需要用到后面進(jìn)階的知識(shí)
-
第三部分：在
MrBayes >后面，輸入mcmc ngen=20000 samplefreq=100 printfreq=100 diagnfreq=1000這段代碼的結(jié)果會(huì)從后驗(yàn)概率分布中得到至少200個(gè)sample，然后每1000 generation計(jì)算一次結(jié)果；如果數(shù)據(jù)集比較大時(shí)，可以將樣本頻率調(diào)低一些，運(yùn)行的長(zhǎng)度更長(zhǎng)?！灸J(rèn)的設(shè)置是sample和print的frequency是500，diagnfreq是5000，ngen為1,000,000】每一行的末尾是估算運(yùn)行剩余時(shí)間
運(yùn)行結(jié)束后，會(huì)給出一個(gè)結(jié)果：average standard deviation of split frequencies。如果得到設(shè)置的20000 generations后，split frequency 是小于0.01的，就可以停止運(yùn)行了（輸入
no）；否則繼續(xù)運(yùn)行，直到這個(gè)值變成小于0.01。當(dāng)然，0.01是一個(gè)比較嚴(yán)格的閾值，如果只對(duì)樹(shù)中構(gòu)建的比較好的部分感興趣，那么split frequency小于0.05也是可以的 -
第四部分：輸入
sump就是summary of parameter，包括了mean, mode, 95% HPD (Highest Posterior Density) interval需要注意的是，PSRF (Potential Scale Reduction Factor) 應(yīng)該在所有參數(shù)中接近1；如果差別比較大，需要設(shè)置分析長(zhǎng)度再長(zhǎng)一些
結(jié)果接下來(lái)輸入
sumt，會(huì)根據(jù)后驗(yàn)概率對(duì)每個(gè)分支生成一個(gè)進(jìn)化分支圖，也會(huì)根據(jù)平均分枝長(zhǎng)度生成系統(tǒng)發(fā)育圖，這兩個(gè)樹(shù)圖可以輸出成被FigTree、TreeView、Mesquite等畫樹(shù)軟件識(shí)別的文件
詳細(xì)點(diǎn)的解讀
第一步 讀取數(shù)據(jù)
數(shù)據(jù)格式需要是Nexus，包含比對(duì)好的核苷酸或氨基酸序列、形態(tài)學(xué)數(shù)據(jù)、限制性位點(diǎn)數(shù)據(jù)或者這四種數(shù)據(jù)的混合?！娟P(guān)于Nexus格式的介紹：http://wiki.christophchamp.com/index.php?title=NEXUS_file_format】
Nexus數(shù)據(jù)中，不同的數(shù)據(jù)支持的字符是規(guī)定好的，例如
- DNA數(shù)據(jù)
A, C, G, T, R, Y, M, K, S, W, H, B, V, D, N
-
RNA數(shù)據(jù)
A, C, G, U, R, Y, M, K, S, W, H, B, V, D, N -
蛋白數(shù)據(jù)
A, R, N, D, C, Q, E, G, H, I, L, K, M, F, P, S, T, W, Y, V, X -
限制性位點(diǎn)數(shù)據(jù)
0, 1 -
形態(tài)學(xué)數(shù)據(jù)
0, 1, 2, 3, 4, 5, 6, 5, 7, 8, 9
如果出現(xiàn)某些模糊的位點(diǎn)，可以用圓括號(hào)或者花括號(hào)注釋，例如
# 不例如清楚氨基酸是A還是F
(AF), (A,F), {AF}, {A,F}
如果使用自己的真實(shí)數(shù)據(jù)，在導(dǎo)入數(shù)據(jù)后，需要將自己的數(shù)據(jù)調(diào)成和示例一樣的格式，就是類似這種。自己需要修改
ntax（總共要聚類的sample數(shù)）,nchar（每個(gè)sample的序列字符數(shù)），datatype（DNA/RNA/protein等）。如果有missing或者gap，需要說(shuō)明他們用什么字符表示（比如：missing=*gap=-）
使用execute或者exe + 文件名（有路徑的需要寫路徑） 來(lái)加載.nex文件
第二步 指定模型【重要的部分】
在加載數(shù)據(jù)以后，我們可以使用showmodel來(lái)看看針對(duì)自己的數(shù)據(jù)類型，有哪些模型可選，然后需要用到的命令：lset定義模型結(jié)構(gòu)，prset定義模型參數(shù)的先驗(yàn)概率分布
lset：使用help lset 看一下怎么設(shè)置lsetlset <parameter>=<option> ... <parameter>=<option>
還有一個(gè)表格，表示了參數(shù)及當(dāng)前設(shè)置
行首：Model setting for partition 1，表示的是你要比較的序列都是同一類型，如果不同類型就會(huì)分成不同的partition；
第一個(gè)參數(shù) Nucmodel： 設(shè)置核苷酸替換類型。默認(rèn)是4by4（意思是核苷酸的四種形式ACGT/U）；Doublet是核糖體DNA（即編碼rRNA的序列）的成對(duì)莖區(qū) [paired stem regions of ribosomal DNA]； Codon是利用密碼子分析DNA序列；Protein也就是DNA轉(zhuǎn)變的氨基酸序列；
第二個(gè)參數(shù)Nst：用數(shù)字設(shè)置替換的種類。1表示所有替換率都一樣（比如JC69 or F81 model）；2表示所有轉(zhuǎn)換和顛換的比率由一些差別（比如K80 or HKY85 model）；6允許所有的替換率存在差別（如GTR model）；mixed表示在所有可能的可逆替換模型中進(jìn)行“馬爾科夫鏈”抽樣，并組合出不同的模型；
第三個(gè)參數(shù)Code：只有當(dāng)Nocmodel設(shè)置為Codon時(shí)才有效，默認(rèn)就是全部密碼子；
第四個(gè)參數(shù)Ploidy：設(shè)置染色體倍數(shù)，"Haploid", "Diploid" or "Zlinked"；
第五個(gè)參數(shù)Rates: 設(shè)置模型的位點(diǎn)變異率。默認(rèn)是equal（即所有位點(diǎn)沒(méi)有變異）；gamma表示位點(diǎn)變異符合gamma分布；lnorm表示正態(tài)分布的對(duì)數(shù)；propinv表示變異位點(diǎn)的分布是常數(shù)；invgamma表示變異位點(diǎn)的分布是常數(shù)，但是其他位點(diǎn)分布是gamma分布；
第六、七個(gè)參數(shù)Ngammacat 、Nbetacat一般用默認(rèn)值4、5就好（大多數(shù)情況下，Ngammacat取4個(gè)rate categories就夠了，增加rate categories數(shù)量可以增加準(zhǔn)確度，但同時(shí)速度會(huì)變慢，4可以理解為一個(gè)折中的值）；
第三步 設(shè)置模型先驗(yàn)值及檢查模型
有六種參數(shù)類型：拓?fù)洌╰he topology）, 分枝長(zhǎng)度（the branch lengths）, 4個(gè)stationary frequencies of the nucleotides, 6個(gè)different nucleotide substitution rates, the proportion of invariable sites 以及 the shape parameter of the gamma distribution of rate variation
默認(rèn)值可以應(yīng)對(duì)絕大多數(shù)分析，使用help prset可以看到設(shè)置
使用shomodel會(huì)看到目前模型的參數(shù)設(shè)置
第四步 準(zhǔn)備分析
在運(yùn)行mcmc 之前，可以用help mcmc看下
- Ngen：generation的數(shù)量，可以先設(shè)置個(gè)小一點(diǎn)的值確保整個(gè)流程是正確的，然后還可以根據(jù)這個(gè)結(jié)果估算更大的Ngen使用時(shí)間。這里的測(cè)試數(shù)據(jù)我們使用20000即可，但是大數(shù)據(jù)集可以適當(dāng)調(diào)小一些。使用
mcmcp ngen=20000進(jìn)行調(diào)整，然后再用help mcmc確認(rèn) - Diagnfreq：每隔多少generation進(jìn)行一次評(píng)估，對(duì)大數(shù)據(jù)集來(lái)講設(shè)為5000比較合適，這里測(cè)試的12個(gè)樣本設(shè)置1000比較合適
mcmcp diagnfreq=1000
這個(gè)結(jié)果會(huì)保存在.mcmc文件中。每次評(píng)估結(jié)果出來(lái)，都會(huì)從馬爾可夫鏈開(kāi)頭根據(jù)burnin（固定的要被丟棄的樣本數(shù)量）或者burninfrac（樣本百分比）丟棄樣本。默認(rèn)情況下，cold chain的前25%會(huì)被丟棄（通過(guò)relburnin=yes and burninfrac=0.25設(shè)置） - Nchains：默認(rèn)情況下為4（3個(gè)heated chains，1個(gè)cold chain）一般來(lái)講，增加3個(gè)heated chains對(duì)復(fù)雜的大型數(shù)據(jù)集有幫助，當(dāng)然鏈數(shù)越多對(duì)內(nèi)存消耗越大。在服務(wù)器上，可以使用MPI version的mrbayes軟件來(lái)分配鏈給不同的CPU，提高運(yùn)行速度
- Samplefreq: 鏈被取樣的頻率，默認(rèn)是鏈上每隔500單位被取樣一次。小數(shù)據(jù)增加頻率，也就是減小數(shù)值比如測(cè)試的數(shù)據(jù)中12樣本只需要設(shè)置成
mcmcp samplefreq=100；大的數(shù)據(jù)集可以適當(dāng)增加數(shù)值。當(dāng)鏈被取樣結(jié)束后，參數(shù)信息被儲(chǔ)存在.p的文件中（tab分割）; 拓?fù)浜头种﹂L(zhǎng)度被輸出到.t文件(這些是Nexus tree file，可以導(dǎo)入TreeView) - Printfreq: 就是輸出到屏幕的結(jié)果間隔，默認(rèn)是500
第五步 開(kāi)始分析
輸入mcmc ，首先看到這樣一個(gè)表格，包括了需要用到的方法及占比
然后就開(kāi)始運(yùn)行，就像這樣
第一列是 generation number，后面8列負(fù)數(shù)代表內(nèi)存的物理位置
2-5列是run1的數(shù)據(jù)（之前看到Nruns=2，說(shuō)明有兩個(gè)完全獨(dú)立并同時(shí)進(jìn)行的run），一列的每行對(duì)應(yīng)一個(gè)chain。
最后一列是預(yù)計(jì)剩余時(shí)間
什么時(shí)候停止呢？
當(dāng)程序運(yùn)行到設(shè)置的Ngen時(shí)，會(huì)問(wèn)你要不要繼續(xù)，并給出目前已經(jīng)計(jì)算的結(jié)果Average standard deviation of split frequencies。一般來(lái)講這個(gè)數(shù)小于0.01就可以退出了。但是在數(shù)據(jù)集比較大時(shí)，并不是很理想，下降的很慢，而且開(kāi)始還有時(shí)會(huì)上升，如果到了Ngen還離0.01很遠(yuǎn)，比如最后還是1點(diǎn)幾，那么還需要增加generation。如果達(dá)到了0.01到0.05之間，那么基本可以結(jié)束。
第六步 總結(jié)
首先看下samplefreq生成的.p文件，不同的模型結(jié)果信息不同【這里以測(cè)試數(shù)據(jù)為例】
- 第一行：每個(gè)運(yùn)行結(jié)果獨(dú)立的ID號(hào)；
- 第二行：從左往右依次為：
Gen (generation number);
LnL (log likelihood of the cold chain);
TL (total tree length/ the sum of all branch lengths);
r(A<->C)等 (six GTR rate parameters);
pi (four stationary nucleotide frequencies);
alpha (shape
parameter of the gamma distribution of rate variation);
pinvar (proportion of invariable sites)
然后使用sump 得到一個(gè)圖一個(gè)表。生成的表格中，主要看PSRF這一列數(shù)值是否接近1.0【1.00-1.02是最理想的情況，但很難做到】
tree和branch數(shù)據(jù)存儲(chǔ)在.t 的文件中，也是Nexus格式
使用sumt relburnin =yes burninfrac = 0.25 ，直接返回的結(jié)果包括：兩個(gè)統(tǒng)計(jì)表兩個(gè)圖
-
左上：
obs表示number of times the partition was sampledProbab表示probability of the partition
Sd表示standard deviation of the partition frequency
Min/Max表示min and max of the standard deviation
Nruns表示number of runs 左下：記錄了length數(shù)據(jù)，其中l(wèi)ength的編號(hào)對(duì)應(yīng)左上圖中ID號(hào)，依然關(guān)注PSRF列（Potential Scale Reduction Factor）
右上：clade credibility tree: probability of each partition or clade in the tree
右下：phylogram tree：branch lengths measured in expected substitutions per site
另外，sumt還會(huì)返回5個(gè)額外的文件
.parts文件包含了二分法分類的key值；
.tstat和.vstat 包含了partition statistics和branch length statistics
.con 包含了consensus trees，這個(gè)文件可以在FigTree中打開(kāi)，展示一下后驗(yàn)概率以及每個(gè)枝的標(biāo)準(zhǔn)差【這個(gè)應(yīng)該是比較有用的】
.trprobs包含了在mcmc搜索中找到的樹(shù)，并用后驗(yàn)概率排序
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩，打造生信星球，想讓它成為一個(gè)不拽術(shù)語(yǔ)、通俗易懂的生信知識(shí)平臺(tái)。需要幫助或提出意見(jiàn)請(qǐng)后臺(tái)留言或發(fā)送郵件到Bioplanet520@outlook.com
