多功能系統(tǒng)發(fā)育包設計和實現(xiàn)

摘要1.在這里，我提出了一種新的多功能系統(tǒng)發(fā)育包，phytools，對于R統(tǒng)計計算環(huán)境。2.該軟件包的重點是系統(tǒng)發(fā)育比較生物學的方法；但是，它還包括用于樹木推斷，系統(tǒng)發(fā)育輸入/輸出，繪圖，操縱和其他一些任務的工具。3.我描述并列出了在phytools中實現(xiàn)的主要方法，并以兩個說明性示例的形式提供了其用法的一些演示。4.最后，我通過簡要描述一個有效的Web日志作為結束語，該日志用于記錄phytools的當前和未來開發(fā)。我還注意到R計算環(huán)境中用于系統(tǒng)發(fā)育的其他Web資源。介紹在最近的幾十年中，系統(tǒng)發(fā)育在進化生物學中發(fā)揮了中心作用（Felsenstein1985，2004

;

Harvey＆Pagel1991

;

Losos2011）。在系統(tǒng)進化論者中，科學計算環(huán)境R（RDevelopmentCoreTeam2011）的普及度突飛猛進，特別是自開發(fā)多功能“猿”（系統(tǒng)發(fā)育和進化分析）R軟件包以來（Paradis，Claude和Strimmer，2004年））以及自Paradis的“UseR！”發(fā)布以來

系統(tǒng)發(fā)育書（Paradis2006）。近年來，以許多貢獻的軟件包的形式見證了R的系統(tǒng)發(fā)育能力的迅速擴展。大多數(shù)，例如流行的軟件包“geiger”（Harmon?等。2008年）和“phangorn”（?Schliep2011年），是通過構建在猿類中開發(fā)的功能和數(shù)據(jù)結構來實現(xiàn)的。在本說明中，我將描述一個用R語言編寫的新的系統(tǒng)發(fā)育軟件包。該庫稱為“phytools”（用于比較生物學的系統(tǒng)發(fā)育工具等），可以從CRAN的綜合R存檔網(wǎng)絡安裝。在此一攬子計劃中，我主要致力于通過自己的工作來實現(xiàn)幾種系統(tǒng)發(fā)育比較方法（Revell?2008、2009、2010；Revell，Harmon＆Collar2008；Revell＆Harrison2008；Revell＆Collar2009；Lindenfors，Revell＆Nunn2010；Revell＆Collar2009；Revell＆Collar2009。Mahler?等人，?2010年；Revell?等人，印刷中）;

但是，我還介紹了一些以前在R中未實現(xiàn)的對系統(tǒng)發(fā)育社區(qū)感興趣的新功能（例如，Cavalli-Sforza和Edwards1967；Baum1992；Ragan1992；Nielsen2002；Huelsenbeck，Nielsen＆Bollback2003；Bollback）2006年；O'Meara?等人?2006年；Ives，Midford和Garland2007年；Sidlauskas2008年），以及許多用于讀取，書寫，繪圖和操縱特殊類型的系統(tǒng)發(fā)育樹的簡單實用函數(shù)。我一直都在努力使phytools與猿類包裝的交互性最大化。例如，phytools的新功能之一是生成，繪制，讀取和寫入隨機字符映射樹的能力（Nielsen2002）。我沒有創(chuàng)建一種新的R對象來存儲隨機映射的系統(tǒng)發(fā)育樹，而是直接在為猿而開發(fā)并用于許多其他R系統(tǒng)發(fā)育軟件包的現(xiàn)有“系統(tǒng)”結構上構建。目前，phytools無法與“phylobase”軟件包互操作（RHackathon?等人，?2011），盡管將來會增加此功能。在以下各節(jié)中，我將描述phytools庫的主要功能。我將提供兩個說明性示例，以演示phytools的某些功能。最后，由于phytools尚在開發(fā)中，因此我將描述一個Web日志，以使phytools用戶及時了解該軟件包的錯誤，更新和未來的軟件開發(fā)。描述phytools庫完全使用科學計算語言R（RDevelopmentCoreTeam2011）編寫。它利用了其他軟件包中開發(fā)的功能，特別是核心系統(tǒng)發(fā)育軟件包猿（Paradis，Claude＆Strimmer2004）來進行多種類型的系統(tǒng)發(fā)育樹輸入和操作。它還使用系統(tǒng)發(fā)育推斷包phangorn（Schliep2011）進行推斷和某些其他計算。另外，phytools包直接或通過其依賴項依賴，導入或建議其他幾個R庫。其中包括以下軟件包：動畫：謝（2011）

;

校準：Graffelman（2010）；圖：Csardi＆Nepusz（2006）;

主持人：Genz＆Azzalini（2011）；女士：杰克遜（2011）；numDeriv：Gilbert（2011）；quadprog：Turlach＆Weingessel（2010）。到目前為止，我已經(jīng)為phytools包實現(xiàn)了許多功能。但是，我還應該指出，phytools仍在開發(fā)中，我希望phytools的功能在未來幾年中會大大擴展。在表1中，我提供了迄今為止在phytools庫中實現(xiàn)的主要功能的帶注釋的列表。這些功能涵蓋了系統(tǒng)發(fā)育生物學的幾個不同領域中的方法，稍后將進行介紹。比較生物學系統(tǒng)發(fā)育比較生物學中的幾種方法已在phytools中實現(xiàn)。這些涵蓋了廣泛的領域，包括祖先字符估計（例如anc.trend），用于研究字符特征隨時間演變的基于似然的方法（例如brownie.lite，evol.vcv，fitDiversityModel和phylosig），貝葉斯方法。檢測樹中速率變化的位置（evol.rate.mcmc），估算系統(tǒng)發(fā)育信號，包括采樣誤差（phylosig）以及在系統(tǒng)發(fā)育背景下進行統(tǒng)計假設檢驗的各種方法（例如phyl.cca，phyl。pairedttest，phyl.pca和phyl.resid）。模擬在phytools中實現(xiàn)了幾種模擬方法。這些功能包括在各種條件下的布朗運動仿真（fastBM），離散字符演化的仿真（sim.history），隨機字符映射的仿真（make.simmap）以及多種演化速率的仿真（sim.rates）以及其他功能（表1）。系統(tǒng)發(fā)生推斷在phytools軟件包中實現(xiàn)了幾種不同的系統(tǒng)發(fā)育推斷程序。這些功能通常高度依賴Schliep（2011）的phangorn庫中的計算和算法。其中一些功能包括矩陣表示簡約樹估計（mrp.supertree）和最小二乘系統(tǒng)發(fā)育推斷（optim.phylo.ls；表1）。圖形方法在phytools中實現(xiàn)了幾種圖形方法。其中有一個樹進二元morphospace的投影（phylomorphospace

;

圖1），繪制隨機字符地圖和歷史（plotSimmap通過時間），譜系與滅絕譜系（繪制LTT），布朗運動和形態(tài)（的動畫branching.diffusion）和其他功能（表1）。從Mahler?等人的數(shù)據(jù)獲得的兩個主要成分軸的形態(tài)空間圖。（2010）針對加勒比Anolis蜥蜴。每個點代表一個物種的表型值（彩色）或假設的祖先表型（黑色）。線通過假設的祖先連接相關物種（Sidlauskas2008）。不同的顏色代表了不同生態(tài)形態(tài)學類別的物種，這些物種已根據(jù)微生境利用，生態(tài)學和總體表型相似性進行了描述（Losos2009；圖例插圖中指示的顏色代碼）。系統(tǒng)發(fā)育已被修剪為僅包括這些所謂的“生態(tài)形態(tài)”物種。實用功能除上述科學功能外，phytools還包括許多用于系統(tǒng)發(fā)育輸入，輸出和操作的實用功能。這些是為了補充和補充ape和phangorn軟件包中現(xiàn)有的各種實用功能。表1中列出了其中一些功能。例子為了演示phytools的用法，我創(chuàng)建了兩個簡短的說明性示例，讀者可以輕松復制它們。首先，我使用模擬數(shù)據(jù)和phytools函數(shù)evol.rate.mcmc來確定進化速率隨時間變化的位置（Revell?等人，印刷中）。在第二步中，我模擬了隨機離散字符歷史和以離散字符狀態(tài)為條件的具有不同速率的連續(xù)字符，然后使用phytools函數(shù)brownie.lite擬合了多速率布朗字符演化模型（O'Meara?等（?2006年）。示例1：檢測速率偏移的位置在此示例中，我首先模擬了隨機的純出生系統(tǒng)發(fā)育；接下來，我在樹的不同部分以兩種不同的進化速率模擬系統(tǒng)發(fā)育上單個連續(xù)值字符的進化變化。我使用貝葉斯MCMC方法分析樹和數(shù)據(jù)，以確定進化速率隨時間變化的位置（Revell?等，印刷中）。最后，我分析了MCMC結果，以估計移動的位置以及該點前后的進化速率。首先，我加載了phytools軟件包。這還將在首次實例化時加載ape和其他必需的軟件包：>＃加載phytools軟件包（和猿）

>require（phytools）加載所需的軟件包：phytools正在加載所需的包裹：猿...接下來，我設置隨機數(shù)種子以提高可重復性（在這里，它只是設置為1）：>set.seed（1）＃設置種子我使用猿函數(shù)rbdtree來模擬隨機的純樹。在這種情況下，樹有91個分類單元。>＃模擬一棵樹（使用猿）

>tree<-rbdtree（b=log（50），d=0，Tmax=1）現(xiàn)在，出于仿真目的，我在預定位置分割了樹-在此指定了后代節(jié)點的數(shù)量以及沿邊緣的距離與根的距離。為此，我使用phytools函數(shù)splitTree。應該注意的是，下面使用的節(jié)點和邊緣位置只能保證在將隨機數(shù)種子設置為1（見上文）的條件下工作，否則應選擇不同的分割點。>樹<-splitTree（樹，列表（節(jié)點=153，bp=0·09））現(xiàn)在，我在樹的一部分中拉伸分支，然后重新附加子樹以進行仿真。為了重新連接樹的兩個部分，我使用phytools函數(shù)paste.tree，如下所示：>＃將樹枝伸到樹的一部分>樹[[2]]$edge.length<-樹[[2]]$edge.length*10>＃重新附加兩個子樹>sim.tree<-paste.tree（trees[[1]]，trees[[2]]）然后，我可以使用phytools函數(shù)plotTree繪制生成樹以進行仿真（其分支拉伸成與進化速率乘以時間成比例；圖2），并使用phytools函數(shù)fastBM在該拉伸樹上進行仿真：>＃繪制生成樹以進行仿真>plotTree（sim.tree，fsize=0·5）>x<-fastBM（sim.tree）＃在樹上模擬現(xiàn)在，我使用phytools函數(shù)evol.rate.mcmc進行貝葉斯MCMC分析（Revell?等人，印刷中）。在運行3·20GHz的Delli5650CPU上，此分析大約需要20分鐘。>＃執(zhí)行MCMC>res<-evol.rate.mcmc（tree，x，ngen=100000）控制參數(shù)（由用戶或默認設置）：清單11$sig1：數(shù)字4·34$sig2：數(shù)字4·34$a：num0·0222$sd1：編號0·867$sd2：編號0·867$sda：編號0·00444$kloc：num0·2$sdlnr：num1$rand.shift：數(shù)字0·05$打?。簄um100$樣本：數(shù)量100開始運行MCMC：狀態(tài)sig1sig2節(jié)點bp可能性04·33524·3352·03502220·345

-

130·7022狀態(tài)sig1sig2節(jié)點bp可能性1004·85170·9447

?

0·0248930·044

?

116·3554狀態(tài)sig1sig2節(jié)點bp可能性2006·79151·2872

?

0·02051530·057

?

101·0501完成MCMC運行。MCMC功能首先打印控制參數(shù)（用戶可以設置，盡管上面已經(jīng)給它們提供了默認值，請參見下文），然后以控制參數(shù)print（給定的頻率）打印MCMC鏈的狀態(tài)。在這里，每100代；200以后的幾代未在上面顯示）。接下來，我可以通過尋找與所有其他樣品的最小距離相加后樣品中的分裂估計移點的位置（這是多個可能的標準中的一個;參見雷維爾等人在壓）。為了進行此分析，我使用了phytools函數(shù)minSplit并將前20000代排除為老化：>est.split<-minSplit（tree，res$mcmc[201：nrow（res$mcmc），]）該分析花了大約6秒鐘才能在與之前所述相同的硬件上運行。最后，我需要預先處理后的樣品，以獲得平均偏移的采樣率tipward和rootward，對于每個樣品（參見雷維爾等人在壓）。我使用phytools函數(shù)posterior.evolrate進行此操作。然后，我可以將結果（估計的轉變點和進化速率）打印到屏幕上：>pp.result<-posterior.evolrate（tree，est.split，res$mcmc[201：nrow（res$mcmc），]，res$tips[201：nrow（res$mcmc）]））>est.sig1<-mean（pp.result[，“sig1”]）>est.sig2<-mean（pp.result[，“sig2”]）>est.split$node[1]153$bp[1]0·1003489>est.sig1[1]1·010709>est.sig2[1]10·24914該分析花費了大約9s的時間完成。這里，參數(shù)估計非常接近的產(chǎn)生換檔點[153，0·09]和的產(chǎn)生進化速率和。應當指出的是，在實際操作中，作者應該比這里給出的更多地關注MCMC的控制參數(shù)，尤其是模型參數(shù)的建議分布。可以通過調用evol.rate.mcmc的幫助文件來獲得有關函數(shù)控制的更多信息：>？evol.rate.mcmc或參考Revell?等。（印刷中）。此外，用戶應評估收斂性，并根據(jù)后驗分布計算樣本的有效樣本量。這可以使用MCMC診斷程序包“coda”（Plummer?等人?2006）來完成。請參考Revell?等。（印刷中）以獲取有關此方法的更多信息。示例2：模擬和分析多速率布朗演化在此示例中，我首先通過系統(tǒng)進化來模擬具有三個狀態(tài)的離散值字符特征的字符歷史。然后，我模擬連續(xù)性狀的演化，其速率取決于離散性狀的值。最后，我使用O'Meara?等人的似然方法將單速率和多速率演化模型擬合到數(shù)據(jù)和樹。（2006）。加載phytools之后，我首先設置種子（任意設置為10；在此僅出于可再現(xiàn)性進行此操作）：>set.seed（10）現(xiàn)在，我使用猿模擬隨機的純出生樹：>樹<-rbdtree（b=log（50），d=0，Tmax=1）這棵樹包含129個分類單元。接下來，我使用phytools函數(shù)sim.history模擬樹上具有三個狀態(tài)A，B和C的角色的隨機角色歷史，如下所示：>＃這是我們的轉移矩陣

>Q<

-matrix

（c（

?

2,1,1,1，

?

2,1,1,1，

?

2），3,3）>行名（Q）<-colnames（Q）<-c（“A”，“B”，“C”）>mtree<-sim.history（tree，Q）我可以使用phytools函數(shù)plotSimmap繪制模擬歷史記錄，以查看其外觀：>＃設置顏色>cols<-c（“紅色”，“藍色”，“綠色”）；

>名稱（列）<-行名（Q）>＃關閉標簽的繪圖樹

>plotSimmap（mtree，cols，ftype=“off”）此可視化效果如圖3所示。圖3在圖形查看器中打開PowerPoint示例2的模擬隨機歷史。在這里，不同的彩色分支表示在樹上演化的離散值字符特征的不同模擬狀態(tài)。顏色代碼為A：紅色，B：藍色和C：綠色。接下來，我使用phytools函數(shù)sim.rates使用三種不同的速率來模擬連續(xù)字符的演變：>＃設置速率

>sig2<-c（1,10,100）

>名稱（sig2）<-行名（Q）>X<‐sim.rates（mtree，sig2）＃模擬最后，我使用O'Meara?等人的似然方法擬合了多比率布朗模型。（2006年）具有phytools函數(shù)brownie.lite。這種可能性優(yōu)化大約需要3秒鐘才能在前面所述的相同硬件上運行。>fit.bm<-brownie.lite（mtree，X，maxit=4000）>fit.bm$sig2.single[1]47·15693$a。單[1]

-

5·242448$var.single[1]34·47709$logL1[1]

-

329·7582$k1[1]2$sig2.multipleBCA10·464639999·95633140·8247106...$logL.