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

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

;

Harvey＆Pagel1991

;

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

系統(tǒng)發(fā)育書（Paradis2006）。近年來，以許多貢獻(xiàn)的軟件包的形式見證了R的系統(tǒng)發(fā)育能力的迅速擴(kuò)展。大多數(shù)，例如流行的軟件包“geiger”（Harmon?等。2008年）和“phangorn”（?Schliep2011年），是通過構(gòu)建在猿類中開發(fā)的功能和數(shù)據(jù)結(jié)構(gòu)來實(shí)現(xiàn)的。在本說明中，我將描述一個(gè)用R語言編寫的新的系統(tǒng)發(fā)育軟件包。該庫稱為“phytools”（用于比較生物學(xué)的系統(tǒng)發(fā)育工具等），可以從CRAN的綜合R存檔網(wǎng)絡(luò)安裝。在此一攬子計(jì)劃中，我主要致力于通過自己的工作來實(shí)現(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中未實(shí)現(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í)用函數(shù)。我一直都在努力使phytools與猿類包裝的交互性最大化。例如，phytools的新功能之一是生成，繪制，讀取和寫入隨機(jī)字符映射樹的能力（Nielsen2002）。我沒有創(chuàng)建一種新的R對象來存儲隨機(jī)映射的系統(tǒng)發(fā)育樹，而是直接在為猿而開發(fā)并用于許多其他R系統(tǒng)發(fā)育軟件包的現(xiàn)有“系統(tǒng)”結(jié)構(gòu)上構(gòu)建。目前，phytools無法與“phylobase”軟件包互操作（RHackathon?等人，?2011），盡管將來會(huì)增加此功能。在以下各節(jié)中，我將描述phytools庫的主要功能。我將提供兩個(gè)說明性示例，以演示phytools的某些功能。最后，由于phytools尚在開發(fā)中，因此我將描述一個(gè)Web日志，以使phytools用戶及時(shí)了解該軟件包的錯(cuò)誤，更新和未來的軟件開發(fā)。描述phytools庫完全使用科學(xué)計(jì)算語言R（RDevelopmentCoreTeam2011）編寫。它利用了其他軟件包中開發(fā)的功能，特別是核心系統(tǒng)發(fā)育軟件包猿（Paradis，Claude＆Strimmer2004）來進(jìn)行多種類型的系統(tǒng)發(fā)育樹輸入和操作。它還使用系統(tǒng)發(fā)育推斷包phangorn（Schliep2011）進(jìn)行推斷和某些其他計(jì)算。另外，phytools包直接或通過其依賴項(xiàng)依賴，導(dǎo)入或建議其他幾個(gè)R庫。其中包括以下軟件包：動(dòng)畫：謝（2011）

;

校準(zhǔn)：Graffelman（2010）；圖：Csardi＆Nepusz（2006）;

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

;

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

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

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

-

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

?

0·0248930·044

?

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

?

0·02051530·057

?

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

>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繪制模擬歷史記錄，以查看其外觀：>＃設(shè)置顏色>cols<-c（“紅色”，“藍(lán)色”，“綠色”）；

>名稱（列）<-行名（Q）>＃關(guān)閉標(biāo)簽的繪圖樹

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

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

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