原文鏈接《 Data Integration, Manipulation and Visualization of Phylogenetic Trees》
1、樹文件輸入
1.1系統(tǒng)發(fā)育樹構(gòu)建綜述
系統(tǒng)發(fā)育樹是一種描述一群生物之間親緣關(guān)系的樹，它可以根據(jù)生物的遺傳序列來構(gòu)建。根狀的系統(tǒng)發(fā)育樹表示一種進(jìn)化歷史模型，該模型由樹節(jié)點(diǎn)之間的祖先-后代關(guān)系和不同親緣程度的姐妹或表親生物體的聚類關(guān)系所描述，如圖1.1所示。在傳染病研究中，系統(tǒng)發(fā)育樹通常由病原體基因或基因組序列構(gòu)建，以顯示哪個(gè)病原體樣本在遺傳上更接近另一個(gè)樣本，從而為未觀察到的潛在流行病學(xué)聯(lián)系和暴發(fā)的潛在來源提供見解。
系統(tǒng)發(fā)育樹可以用基于距離的方法或基于字符的方法由遺傳序列構(gòu)建?；诰嚯x的方法，包括UPGMA (unweighted pair group method with算術(shù)均值法)和NJ (Neighbor-joining法)，都是基于序列間的兩兩遺傳距離矩陣?；谧址姆椒ǎㄗ畲蠊?jié)約法(MP) (Fitch 1971)
、最大似然法(ML) (J. Felsenstein 1981)和貝葉斯馬爾科夫鏈蒙特卡羅法(BMCMC) (Rannala and Yang 1996)，基于描述遺傳特征演化的數(shù)學(xué)模型，并根據(jù)其最優(yōu)性準(zhǔn)則尋找最佳的系統(tǒng)發(fā)育樹。
最大簡約(MP)方法假設(shè)進(jìn)化變化是罕見的，并將特征狀態(tài)變化的數(shù)量最小化(例如，DNA替換的數(shù)量)。這一標(biāo)準(zhǔn)類似于奧卡姆的剃刀，即能夠解釋數(shù)據(jù)的最簡單的假設(shè)就是最好的假設(shè)。不加權(quán)簡約法假設(shè)不同字符(核苷酸或氨基酸)的突變是等可能的，而加權(quán)法假設(shè)突變的不均等可能(例如，第三個(gè)密碼子位置比其他密碼子位置更容易發(fā)生突變;而且過渡突變的頻率比反轉(zhuǎn)突變高)。MP方法的概念是直接和直觀的，這可能是它在生物學(xué)家中流行的原因，他們更關(guān)心研究問題，而不是分析的計(jì)算細(xì)節(jié)。然而，這種方法有幾個(gè)缺點(diǎn)，特別是樹推理可能會(huì)受到眾所周知的系統(tǒng)錯(cuò)誤——長分枝吸引(long-branch attraction, LBA)——的影響，它錯(cuò)誤地將遠(yuǎn)親譜系推斷為近親 (Joseph Felsenstein 1978)。這是因?yàn)镸P方法較差地考慮了從現(xiàn)有遺傳數(shù)據(jù)中很難觀察到的許多序列進(jìn)化因素(如逆轉(zhuǎn)和收斂)。
最大似然法(ML)和貝葉斯馬爾可夫鏈蒙特卡羅法(BMCMC)是構(gòu)建系統(tǒng)發(fā)育樹最常用的兩種方法，也是科學(xué)出版物中最常用的方法。ML和BMCMC方法需要序列演化的替代模型。不同的序列數(shù)據(jù)有不同的替代模型來表述DNA、密碼子和氨基酸的進(jìn)化過程。目前有幾種核苷酸取代模型，包括JC69、K2P、F81、HKY和GTR (Arenas 2015)。這些模型可結(jié)合各站點(diǎn)間的比率變化(記作+ Γ) (Z. Yang 1994)和不變站點(diǎn)的比例(記作+I) (Shoemaker and Fitch 1989)使用。之前的研究 (Lemmon and Moriarty 2004)表明，對(duì)替代模型的錯(cuò)誤說明可能會(huì)導(dǎo)致系統(tǒng)發(fā)育推斷的偏差。推薦對(duì)最適合的替換模型進(jìn)行過程測試。
ML方法的最優(yōu)準(zhǔn)則是找到給定序列數(shù)據(jù)的最大似然樹。ML方法的過程很簡單:計(jì)算一棵樹的可能性，并優(yōu)化其拓?fù)浜头种?以及替代模型參數(shù)，如果不是固定的)，直到找到最好的樹。啟發(fā)式搜索，例如那些在PhyML和RAxML中實(shí)現(xiàn)的搜索，通常用于基于似然準(zhǔn)則尋找最佳樹。貝葉斯方法基于給定的替換模型，通過MCMC對(duì)樹進(jìn)行抽樣，找到后驗(yàn)概率最大的樹。BMCMC的優(yōu)點(diǎn)之一是，在MCMC過程中，可以自然、方便地從采樣樹中獲得包含后驗(yàn)枝概率的參數(shù)方差和樹的拓?fù)洳淮_定性。此外，拓?fù)洳淮_定性對(duì)其他參數(shù)估計(jì)的影響也自然地集成到BMCMC系統(tǒng)發(fā)育框架中。
在一個(gè)簡單的系統(tǒng)發(fā)育樹中，與樹分支/節(jié)點(diǎn)相關(guān)的數(shù)據(jù)可以是分支長度(表明遺傳或時(shí)間發(fā)散)和譜系支持，如bootstrap程序估計(jì)的bootstrap值或BMCMC分析中從采樣樹中總結(jié)的后驗(yàn)進(jìn)化枝概率。
1.2系統(tǒng)發(fā)育樹格式
有幾種文件格式設(shè)計(jì)用于存儲(chǔ)系統(tǒng)發(fā)育樹以及與節(jié)點(diǎn)和分支相關(guān)的數(shù)據(jù)。三種常用的格式是Newick3, NEXUS (Maddison et al. 1997)和Phylip (Joseph Felsenstein 1989)。一些格式(例如，NHX)是從Newick格式擴(kuò)展而來的。大多數(shù)進(jìn)化生物學(xué)軟件都支持Newick和NEXUS格式作為輸入，而一些軟件工具則通過引入新的規(guī)則/數(shù)據(jù)塊來存儲(chǔ)進(jìn)化推理，從而輸出更新的標(biāo)準(zhǔn)文件(如 BEAST
和MrBayes)。在其他情況下(例如，PAML和 r8s)，輸出日志文件只能由它們自己的單一軟件識(shí)別。
1.2.1 Newich 格式
Newick樹格式是用計(jì)算機(jī)可讀形式表示樹的標(biāo)準(zhǔn)格式。
用于演示Newick文本編碼樹結(jié)構(gòu)的示例樹。尖端對(duì)齊到右邊，樹枝的長度被標(biāo)記在每個(gè)樹枝的中間。
圖1.2中所示的根樹可以用以下字符序列表示為Newick樹文本。
((t2:0.04,t1:0.34):0.89,(t5:0.37,(t4:0.03,t3:0.67):0.9):0.59);
Newick格式以分號(hào)結(jié)束。內(nèi)部節(jié)點(diǎn)由一對(duì)匹配的圓括號(hào)表示。括號(hào)之間是該節(jié)點(diǎn)的后代節(jié)點(diǎn)。例如(t2:0.04,t1:0.34)表示t2和t1的父節(jié)點(diǎn)，它們是直系后代。同級(jí)節(jié)點(diǎn)由逗號(hào)分隔，提示由名稱表示。分支長度(從父節(jié)點(diǎn)到子節(jié)點(diǎn))由子節(jié)點(diǎn)之后和冒號(hào)前面的實(shí)數(shù)表示。與內(nèi)部節(jié)點(diǎn)或分支相關(guān)聯(lián)的奇異數(shù)據(jù)(例如，引導(dǎo)值)可以編碼為節(jié)點(diǎn)標(biāo)簽，并用冒號(hào)前的簡單文本/數(shù)字表示。
Newick樹的格式是由Meacham在1984年開發(fā)的 PHYLIP (Retief 2000)包。Newick格式現(xiàn)在是最廣泛使用的樹格式,使用的是 PHYLIP, PAUP* (Wilgenbusch and Swofford 2003), TREE-PUZZLE (Schmidt et al. 2002), MrBayes和許多其他應(yīng)用。Phylip tree格式包含了一種基于MSA序列在同一文件中構(gòu)建的相應(yīng)Newick樹文本的Phylip多重序列對(duì)齊(MSA)。
1.2.2 NEXUS tree format
NEXUS格式 (Maddison et al. 1997)
將Newick樹形文本與相關(guān)信息組織成稱為塊的獨(dú)立單元。NEXUS塊具有以下結(jié)構(gòu)
#NEXUS
...
BEGIN characters;
...
END;
上面的示例樹可以保存為以下NEXUS格式:
#NEXUS
[R-package APE, Fri Dec 31 09:56:12 2021]
BEGIN TAXA;
DIMENSIONS NTAX = 5;
TAXLABELS
t2
t1
t5
t4
t3
;
END;
BEGIN TREES;
TRANSLATE
1 t2,
2 t1,
3 t5,
4 t4,
5 t3
;
TREE * UNTITLED = [&R] ((1:0.04,2:0.34):0.89,(3:0.37,(4:0.03,5:0.67):0.9):0.59);
END;
可以使用方括號(hào)放置注釋。大多數(shù)程序都能識(shí)別出一些塊，包括TAXA(包含分類單元信息)、DATA(包含數(shù)據(jù)矩陣，如序列比對(duì))和TREE(包含系統(tǒng)發(fā)育樹，即Newick樹文本)。值得注意的是，塊可以非常多樣化，其中一些塊只能被一個(gè)特定的程序識(shí)別。例如，由 PAUP*導(dǎo)出的NEXUS文件有一個(gè)包含PAUP*命令的PAUP塊，而 FigTree導(dǎo)出的NEXUS文件有一個(gè)包含可視化設(shè)置的FigTree塊。NEXUS將不同類型的數(shù)據(jù)組織成不同的塊，而支持讀取NEXUS的程序可以解析一些它們識(shí)別到的塊，并忽略那些它們無法識(shí)別的塊。這是一種很好的機(jī)制，允許不同的程序使用相同的格式，而不會(huì)在出現(xiàn)不支持的數(shù)據(jù)類型時(shí)崩潰。值得注意的是，大多數(shù)程序只支持TAXA、DATA和TREE塊的解析，因此，能夠從NEXUS中通用地讀取所有類型的數(shù)據(jù)塊的程序/平臺(tái)將對(duì)系統(tǒng)發(fā)育數(shù)據(jù)集成非常有用。
數(shù)據(jù)塊廣泛用于存儲(chǔ)序列對(duì)齊。為此，用戶可以以Phylip格式存儲(chǔ)樹和序列數(shù)據(jù)，這實(shí)際上分別是Phylip多序列對(duì)齊和Newick樹文本。它用于系統(tǒng)發(fā)育分析包(PHYLIP)。
1.2.3新Hampshire擴(kuò)展格式
Newick、NEXUS和phylip主要用于存儲(chǔ)系統(tǒng)發(fā)育樹和與內(nèi)部節(jié)點(diǎn)或分支相關(guān)的基本奇異數(shù)據(jù)（ basic singular data）。除了分支和節(jié)點(diǎn)上的奇異數(shù)據(jù)注釋(上面提到過)，基于Newick(也稱為New Hampshire)的新Hampshire擴(kuò)展格式(New Hampshire eXtended format, NHX)被開發(fā)出來，用來引入標(biāo)簽，將多個(gè)數(shù)據(jù)字段與樹節(jié)點(diǎn)(內(nèi)部節(jié)點(diǎn)和提示)關(guān)聯(lián)起來。標(biāo)簽放置在分支長度之后，必須在[&&NHX and]之間包裝，這使得它可能與NEXUS格式兼容，因?yàn)樗x了[和]之間的注釋字符。NHX也是 PHYLDOG (Boussau et al. 2013) 和 RevBayes (H?hna et al. 2016)。r (ATV) (Zmasek and Eddy 2001)是一個(gè)java工具，它支持顯示以NHX格式存儲(chǔ)的注釋數(shù)據(jù)，但是這個(gè)包不再更新了。
下面是一個(gè)NHX格式的文檔
(((ADH2:0.1[&&NHX:S=human], ADH1:0.11[&&NHX:S=human]):0.05
[&&NHX:S=primates:D=Y:B=100],ADHY:0.1[&&NHX:S=nematode],
ADHX:0.12[&&NHX:S=insect]):0.1[&&NHX:S=metazoa:D=N],(ADH4:0.09
[&&NHX:S=yeast],ADH3:0.13[&&NHX:S=yeast],ADH2:0.12[&&NHX:S=yeast],
ADH1:0.11[&&NHX:S=yeast]):0.1[&&NHX:S=Fungi])[&&NHX:D=N];
1.2.4 Jplace 格式
為了存儲(chǔ)映射到系統(tǒng)發(fā)育樹上的NGS短序列(用于宏基因組分類)，Matsen提出了這種系統(tǒng)發(fā)育位置的jplace格式 (Frederick A. Matsen et al. 2012)。Jplace格式基于JSON，包含四個(gè)鍵:樹、字段、位置、元數(shù)據(jù)和版本。樹值包含從Newick樹格式擴(kuò)展而來的樹文本，方法是將邊緣標(biāo)簽放在括號(hào)中(如果可用)，然后在分支長度之后，將邊緣編號(hào)放在花括號(hào)中，放在邊緣標(biāo)簽之后。字段值包含位置數(shù)據(jù)的頭信息。位置的值是一個(gè)pquery列表。每個(gè)pquery包含兩個(gè)鍵:p表示位置，n表示名稱或nm表示具有多重性的名稱。p的值是該pquery的位置列表。
下面是一個(gè)jplace格式文件
{
"tree": "(((((((A:4{1},B:4{2}):6{3},C:5{4}):8{5},D:6{6}):
3{7},E:21{8}):10{9},((F:4{10},G:12{11}):14{12},H:8{13}):
13{14}):13{15},((I:5{16},J:2{17}):30{18},(K:11{19},
L:11{20}):2{21}):17{22}):4{23},M:56{24});",
"placements": [
{"p":[24, -61371.300778, 0.333344, 0.000003, 0.003887],
"n":["AA"]
},
{"p":[[1, -61312.210786, 0.333335, 0.000001, 0.000003],
[2, -61322.210823, 0.333322, 0.000003, 0.000003],
[3, -61352.210823, 0.333322, 0.000961, 0.000003]],
"n":["BB"]
},
{"p":[[8, -61312.229128, 0.200011, 0.000001, 0.000003],
[9, -61322.229179, 0.200000, 0.000003, 0.000003],
[10, -61342.229223, 0.199992, 0.000003, 0.000003]],
"n":["CC"]
}
],
"metadata": {"info": "a jplace sample file"},
"version" : 2,
"fields": ["edge_num", "likelihood", "like_weight_ratio",
"distal_length", "pendant_length"
]
}
Jplace是 PPLACER (Frederick A. Matsen, Kodner, and Armbrust 2010b)和 (EPA) (Berger, Krompass, and Stamatakis 2011)的輸出格式。但是這兩個(gè)程序不包含可視化工具。PPLACER提供placeviz來將jplace文件轉(zhuǎn)換為系統(tǒng)演化xml或Newick格式，Archaeopteryx可以將其可視化。
1.2.5軟件輸出
RAxML (Stamatakis 2014)可以通過將bootstrap值存儲(chǔ)為內(nèi)部節(jié)點(diǎn)標(biāo)簽來輸出Newick格式。RAxML支持的另一種方式是將bootstrap值放在方括號(hào)內(nèi)和分支長度之后。這不能被大多數(shù)支持Newick格式的軟件支持，方括號(hào)將被忽略。
BEAST (Bouckaert et al. 2014)的輸出是基于NEXUS的，并且在樹塊中引入了方括號(hào)來存儲(chǔ)BEAST推斷出的進(jìn)化證據(jù)。在括號(hào)內(nèi)，如果特征值的長度超過1(例如HPD或替代率范圍)，也可以加入花括號(hào)。這些括號(hào)放在節(jié)點(diǎn)長度和分支長度之間(即，如果存在，則放在label之后，冒號(hào)之前)。括號(hào)不是Newick格式定義的，它是NEXUS注釋的保留字符。因此，對(duì)于標(biāo)準(zhǔn)的NEXUS解析器，這些信息將被忽略。
下面是一個(gè)BEAST模塊輸出格式
TREE * TREE1 = [&R] (((11[&length=9.47]:9.39,14[&length=6.47]:6.39)[&length=25.72]:25.44,4[&length=9.14]:8.82)[&length=3.01]:3.1,
(12[&length=0.62]:0.57,(10[&length=1.6]:1.56,(7[&length=5.21]:5.19,
((((2[&length=3.3]:3.26,(1[&length=1.34]:1.32,(6[&length=0.85]:0.83,
13[&length=0.85]:0.83)[&length=2.5]:2.49)[&length=0.97]:0.94)
[&length=0.5]:0.5,9[&length=1.76]:1.76)[&length=2.41]:2.36,
8[&length=2.19]:2.11)[&length=0.27]:0.24,(3[&length=3.33]:3.31,
(15[&length=5.29]:5.27,5[&length=3.29]:3.27)[&length=1.04]:1.04)
[&length=1.98]:2.04)[&length=2.83]:2.84)[&length=5.39]:5.37)
[&length=2.02]:2)[&length=4.35]:4.36)[&length=0];
BEAST輸出可以包含許多不同的進(jìn)化推理，這取決于在BEAUTi中定義的用于運(yùn)行的分析模型。例如在分子時(shí)鐘分析中，它包含速率、長度、高度、后驗(yàn)和相應(yīng)的HPD以及用于不確定度估計(jì)的范圍。速率是該分支的估計(jì)進(jìn)化速率。長度是指樹枝的長度，以年為單位。高度為節(jié)點(diǎn)到根的時(shí)間，后驗(yàn)為貝葉斯演化支可信度值。上面的例子是分子時(shí)鐘分析的輸出樹，應(yīng)該包含這些推論。為了節(jié)省空間，上面只顯示了長度估計(jì)。此外，分子進(jìn)化遺傳學(xué)分析 (MEGA) (Kumar, Stecher, and Tamura 2016)也支持以BEAST兼容的Nexus格式導(dǎo)出樹(見1.3.2節(jié))。
MrBayes (Huelsenbeck and Ronquist 2001)是一個(gè)程序，它使用馬爾可夫鏈蒙特卡羅方法從后驗(yàn)概率分布中采樣。它的輸出文件用兩組方括號(hào)分別注釋節(jié)點(diǎn)和分支。例如下面，節(jié)點(diǎn)的后驗(yàn)分支概率和分支長度估計(jì)
tree con_all_compat = [&U] (8[&prob=1.0]:2.94e-1[&length_mean=2.9e-1],10[&prob=1.0]:2.25e-1[&length_mean=2.2e-1],((((1[&prob=1.0]:1.43e-1
[&length_mean=1.4e-1],2[&prob=1.0]:1.92e-1[&length_mean=1.9e-1])[&prob=1.0]:
1.24e-1[&length_mean=1.2e-1],9[&prob=1.0]:2.27e-1[&length_mean=2.2e-1])
[&prob=1.0]:1.72e-1[&length_mean=1.7e-1],12[&prob=1.0]:5.11e-1
[&length_mean=5.1e-1])[&prob=1.0]:1.76e-1[&length_mean=1.7e-1],
(((3[&prob=1.0]:5.46e-2[&length_mean=5.4e-2],(6[&prob=1.0]:1.03e-2
[&length_mean=1.0e-2],7[&prob=1.0]:7.13e-3[&length_mean=7.2e-3])[&prob=1.0]:
6.93e-2[&length_mean=6.9e-2])[&prob=1.0]:6.03e-2[&length_mean=6.0e-2],
(4[&prob=1.0]:6.27e-2[&length_mean=6.2e-2],5[&prob=1.0]:6.31e-2
[&length_mean=6.3e-2])[&prob=1.0]:6.07e-2[&length_mean=6.0e-2])[&prob=1.0]:,
1.80e-1[&length_mean=1.8e-1]11[&prob=1.0]:2.37e-1[&length_mean=2.3e-1])
[&prob=1.0]:4.05e-1[&length_mean=4.0e-1])[&prob=1.0]:1.16e+000
[&length_mean=1.162699558201079e+000])[&prob=1.0][&length_mean=0];
為了節(jié)省空間，刪除了大部分推理，只包含分枝概率的prob和分枝長度的平均值的length_mean。此文件的完整版本還包含prob_stddev、prob_range、probb (percent)、probb +-sd(概率推斷)和length_median、length_95%_HPD(每個(gè)分支)。
BEAST和MrBayes的輸出預(yù)計(jì)將被支持NEXUS的軟件不進(jìn)行推理分析(作為注釋刪除)。FigTree支持對(duì)BEAST和MrBayes輸出進(jìn)行解析，并提供可用于在樹上顯示或注釋的推理。但是從那里提取這些數(shù)據(jù)進(jìn)行進(jìn)一步的分析仍然是一個(gè)挑戰(zhàn)。
HyPhy (Pond, Frost, and Muse 2005)可以進(jìn)行大量的系統(tǒng)發(fā)育分析，包括祖先序列重建。對(duì)于祖先序列重建，這些序列和Newick樹文本以NEXUS格式存儲(chǔ)，作為主要的分析輸出。它沒有完全遵循NEXUS的定義，只是將祖先節(jié)點(diǎn)標(biāo)簽而不是外部節(jié)點(diǎn)標(biāo)簽放在TAXA中。MATRIX塊包含祖先節(jié)點(diǎn)的序列比對(duì)，不能被轉(zhuǎn)回存儲(chǔ)在TREES塊中的樹，因?yàn)樗话?jié)點(diǎn)標(biāo)簽。下面是示例輸出(為了節(jié)省空間，只顯示前72bp的對(duì)齊位置)
#NEXUS
[
Generated by HYPHY 2.0020110620beta(MP) for MacOS(Universal Binary)
on Tue Dec 23 13:52:34 2014
]
BEGIN TAXA;
DIMENSIONS NTAX = 13;
TAXLABELS
'Node1' 'Node2' 'Node3' 'Node4' 'Node5' 'Node12' 'Node13' 'Node15'
'Node18' 'Node20' 'Node22' 'Node24' 'Node26' ;
END;
BEGIN CHARACTERS;
DIMENSIONS NCHAR = 2148;
FORMAT
DATATYPE = DNA
GAP=-
MISSING=?
NOLABELS
;
MATRIX
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATTGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAATGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
ATGGAAGACTTTGTGCGACAGTGCTTCAATCCAATGATCGTCGAGCTTGCGGAAAAGGCAATGAAAGAATAT
END;
BEGIN TREES;
TREE tree = (K,N,(D,(L,(J,(G,((C,(E,O)),(H,(I,(B,(A,(F,M)))))))))));
END;
還有其他一些應(yīng)用程序輸出豐富的信息文本，其中也包含與相關(guān)數(shù)據(jù)相關(guān)的系統(tǒng)發(fā)育樹。例如， r8s (Sanderson 2003)在其日志文件中輸出三棵樹，即TREE、RATE和PHYLO，分別用于按時(shí)間、替換率和絕對(duì)替換縮放的分支。
最大似然系統(tǒng)發(fā)育分析 (PAML) (Ziheng Yang 2007)是一個(gè)用于DNA或蛋白質(zhì)序列系統(tǒng)發(fā)育分析的程序包。兩個(gè)主要程序，BASEML和CODEML，實(shí)現(xiàn)了各種模型。BASEML使用許多可用的核苷酸替換模型(包括JC69、K80、F81、F84、HKY85、T92、TN93和GTR)來估計(jì)樹拓?fù)浣Y(jié)構(gòu)、分支長度和替換參數(shù)。CODEML估計(jì)了密碼子替代模型下的同義和非同義替換率、陽性選擇的似然比檢驗(yàn)(Goldman and Yang 1994)。
BASEML輸出mlb文件，其中包含輸入序列(taxa)比對(duì)和一個(gè)具有分支長度的系統(tǒng)發(fā)育樹，以及替換模型和其他估計(jì)參數(shù)。補(bǔ)充結(jié)果文件rst包含序列比對(duì)(若重建祖先序列則與祖先序列比對(duì))和比對(duì)中每個(gè)位點(diǎn)演化的樸素經(jīng)驗(yàn)貝葉斯(Naive Empirical Bayes, NBE)概率。 CODEML輸出mlc文件，該文件包含樹形結(jié)構(gòu)和同義詞和非同義詞替換率的估計(jì)。CODEML還輸出一個(gè)補(bǔ)充的結(jié)果文件rst，該文件與BASEML類似，不同之處在于該位點(diǎn)被定義為密碼子而不是核苷酸。解析這些文件可能是乏味的，通常需要許多后期處理步驟，并需要編程方面的專業(yè)知識(shí)(例如使用Python5或Perl6)。
引入方括號(hào)用于存儲(chǔ)額外信息是很常見的，包括用于存儲(chǔ)bootstrap值的RAxML、用于注釋的NHX格式、用于邊緣標(biāo)簽的jplace格式、用于進(jìn)化估計(jì)的BEAST格式等。但不同軟件中方括號(hào)的位置不一致，內(nèi)容使用不同的規(guī)則存儲(chǔ)關(guān)聯(lián)數(shù)據(jù)，這些屬性使關(guān)聯(lián)數(shù)據(jù)難以解析。對(duì)于大多數(shù)軟件來說，他們只會(huì)忽略方括號(hào)，只在文件是兼容的情況下解析樹結(jié)構(gòu)。其中一些包含無效字符(如jplace格式的tree字段中的花括號(hào))，甚至樹結(jié)構(gòu)無法被標(biāo)準(zhǔn)解析器解析。
從不同的進(jìn)化推理軟件產(chǎn)生的不同分析結(jié)果中提取有用的系統(tǒng)發(fā)育/分類單元相關(guān)信息顯示在同一棵系統(tǒng)發(fā)育樹上，并進(jìn)行進(jìn)一步的分析是很困難。FigTree支持BEAST輸出，但不支持其他大多數(shù)包含進(jìn)化推理或相關(guān)數(shù)據(jù)的軟件輸出。對(duì)于那些輸出豐富的文本文件(如r8s、pml等)，樹狀結(jié)構(gòu)無法被任何樹查看軟件解析，用戶需要專門知識(shí)從輸出文件中手工提取系統(tǒng)發(fā)育樹和其他有用的結(jié)果數(shù)據(jù)。但這種手工操作速度慢，容易出錯(cuò)。
從該領(lǐng)域常用軟件的不同分析輸出中提取進(jìn)化數(shù)據(jù)的系統(tǒng)發(fā)育樹并不容易。其中的一些(如PAML產(chǎn)出和jplace文件)不需要軟件或編程庫來解析文件的支持,而其他人(例如, BEAST
和MrBayes輸出)可以解析沒有進(jìn)化推理時(shí)存儲(chǔ)在方括號(hào)將省略的評(píng)論,大多數(shù)軟件。雖然FigTree支持可視化由BEAST和MrBayes推斷出的進(jìn)化統(tǒng)計(jì)數(shù)據(jù)，但是提取這些數(shù)據(jù)以進(jìn)行進(jìn)一步的分析是不支持的。不同的軟件包為不同的分析實(shí)現(xiàn)不同的算法(例如，PAML用于dN/dS, HyPhy用于祖先序列，BEAST用于skyline分析)。因此，在處理基因組序列數(shù)據(jù)時(shí)，需要高效靈活地整合不同的分析推理結(jié)果，以便全面理解、比較和進(jìn)一步分析。這促使我們開發(fā)程序庫來解析來自不同來源的系統(tǒng)發(fā)育樹和數(shù)據(jù)。
1.3使用treeio獲取樹數(shù)據(jù)
系統(tǒng)發(fā)育樹是描述物種進(jìn)化關(guān)系的常用方法。與分類單元種/株相關(guān)的信息可以在系統(tǒng)發(fā)育樹描述的進(jìn)化歷史的背景下進(jìn)一步分析。例如，可以研究樹中流感病毒株的宿主信息，以了解病毒系的宿主范圍。此外，與分類單元毒株直接相關(guān)的這些元數(shù)據(jù)(如隔離宿主、時(shí)間、地點(diǎn)等)也經(jīng)常受到進(jìn)一步進(jìn)化或比較系統(tǒng)發(fā)育模型和分析的影響，以推斷其與病毒進(jìn)化或傳播過程相關(guān)的動(dòng)態(tài)。所有這些元數(shù)據(jù)或其他表型數(shù)據(jù)或?qū)嶒?yàn)數(shù)據(jù)都存儲(chǔ)為與節(jié)點(diǎn)或分支相關(guān)的注釋數(shù)據(jù)，并且通常由不同的分析程序以不一致的格式產(chǎn)生。
R語言讀取樹文件仍然是有限制的。Newick和Nexus可以通過幾個(gè)包導(dǎo)入，包括 ape, phylobase。NeXML格式可以被 RNeXML解析。然而，該領(lǐng)域廣泛使用的軟件包的分析結(jié)果并沒有得到很好的支持。SIMMAP輸出可以被phyext2和 phytools解析。雖然 PHYLOCH可以導(dǎo)入BEAST和MrBayes的輸出文件，但只能解析內(nèi)部節(jié)點(diǎn)屬性，而tip屬性則會(huì)被忽略掉。許多其他軟件輸出主要是需要編程專家來導(dǎo)入帶有相關(guān)數(shù)據(jù)的樹。將外部數(shù)據(jù)(包括實(shí)驗(yàn)和臨床數(shù)據(jù))與系統(tǒng)發(fā)育聯(lián)系起來是進(jìn)化生物學(xué)家面臨的另一個(gè)障礙。
為了填補(bǔ)大多數(shù)樹格式或軟件輸出無法在同一軟件/平臺(tái)中解析的空白，開發(fā)了一個(gè)R包 treeio (Wang et al. 2020)，用于解析來自通用進(jìn)化分析軟件的各種樹文件格式和輸出。treeio包是用R編程語言(R Core Team 2016)開發(fā)的。不僅可以解析樹結(jié)構(gòu)，還可以解析相關(guān)的數(shù)據(jù)和進(jìn)化推斷，包括NHX注釋、時(shí)鐘速率推斷(來自BEAST)或r8s (Sanderson 2003)程序)、同義和非同義替換(來自CODEML)、祖先序列構(gòu)建(來自HyPhy、BASEML或CODEML)等。目前，treeio能夠讀取以下文件格式:Newick, Nexus, New Hampshire eXtended format (NHX)， jplace和Phylip，以及以下分析程序的數(shù)據(jù)輸出:ASTRAL, BEAST, EPA, HyPhy, MEGA, MrBayes, pml, PHYLDOG, PPLACER, r8s, RAxML和RevBayes等。這可以通過treeio中開發(fā)的幾個(gè)解析器函數(shù)實(shí)現(xiàn)(表1.1) (Wang et al. 2020)。
| Parser function | Description |
|---|---|
| read.astral | parsing output of ASTRAL |
| read.beast | parsing output of BEAST |
| read.codeml | parsing output of CodeML (rst and mlc files) |
| read.codeml_mlc | parsing mlc file (output of CodeML) |
| read.fasta | parsing FASTA format sequence file |
| read.hyphy | parsing output of HYPHY |
| read.hyphy.seq | parsing ancestral sequences from HYPHY output |
| read.iqtree | parsing IQ-Tree Newick string, with the ability to parse SH-aLRT and UFBoot support values |
| read.jplace | parsing jplace file including the output of EPA and pplacer |
| read.jtree | parsing jtree format |
| read.mega | parsing MEGA Nexus output |
| read.mega_tabular | parsing MEGA tabular output |
| read.mrbayes | parsing output of MrBayes |
| read.newick | parsing Newick string, with the ability to parse node label as support values |
| read.nexus | parsing standard NEXUS file (re-exported from ape) |
| read.nhx | parsing NHX file including the output of PHYLDOG and RevBayes |
| read.paml_rst | parsing rst file (output of BaseML or CodeML) |
| read.phylip | parsing phylip file (phylip alignment + Newick string) |
| read.phylip.seq | parsing multiple sequence alignment from phylip file |
| read.phylip.tree | parsing newick string from phylip file |
| read.phyloxml | parsing phyloXML file |
| read.r8s | parsing output of r8s |
| read.raxml | parsing output of RAxML |
| read.tree | parsing newick string (re-exported from ape) |
TABLE 1.1: treeio中定義的解析器函數(shù)
treeio包定義了用于系統(tǒng)發(fā)育樹輸入和輸出的基類和函數(shù)?？梢詫⒊Ｓ密浖@得進(jìn)化證據(jù)推斷數(shù)據(jù)用于R語言分析的基礎(chǔ)程序 。例如,dN / dS值或祖先序列推斷 CODEML (Ziheng Yang 2007),進(jìn)化枝支持值(后)推斷 BEAST (Bouckaert et al. 2014)和短期讀取位置由 EPA (Berger, Krompass, and Stamatakis 2011)和 PPLACER (Frederick A. Matsen, Kodner, and Armbrust 2010a)。這些進(jìn)化證據(jù)可以在R中進(jìn)一步分析，并使用 ggtree (Yu et al. 2017)來注釋系統(tǒng)發(fā)育樹。隨著分析工具和模型的不斷發(fā)展，整合來自不同來源的不同種類的數(shù)據(jù)和分析結(jié)果以在相同的系統(tǒng)發(fā)育樹背景下進(jìn)行整合分析是一個(gè)挑戰(zhàn)。treeio包 (Wang et al. 2020)提供了一個(gè)merge_tree函數(shù)，允許組合來自不同來源的樹數(shù)據(jù)。此外，treeio還允許將外部數(shù)據(jù)鏈接到系統(tǒng)發(fā)育樹結(jié)構(gòu)。
解析之后，將通過一個(gè)在tidytree包中定義的S4類treedata來存儲(chǔ)具有關(guān)聯(lián)數(shù)據(jù)的樹結(jié)構(gòu)。這些被解析的數(shù)據(jù)被映射到treedata對(duì)象內(nèi)部的樹分支和節(jié)點(diǎn)，因此可以有效地使用 ggtree (Yu et al. 2017)和 ggtreeExtra (Xu, Dai, et al. 2021)可視化地對(duì)樹進(jìn)行注釋。用于系統(tǒng)發(fā)育數(shù)據(jù)解析、集成和注釋的可編程平臺(tái)使我們更容易識(shí)別進(jìn)化動(dòng)力學(xué)和相關(guān)模式(圖1.3，圖1 (Wang et al. 2020))。
1.3.1 treeio概述
treeio包 (Wang et al. 2020)定義了S4類，用于存儲(chǔ)系統(tǒng)發(fā)育樹，這些樹包含來自不同來源的不同類型的相關(guān)數(shù)據(jù)或協(xié)變量，包括來自不同軟件包的分析輸出。它還定義了相應(yīng)的解析器函數(shù)，用于用注釋數(shù)據(jù)解析系統(tǒng)發(fā)育樹，并將它們作為數(shù)據(jù)對(duì)象存儲(chǔ)在R中，以便進(jìn)一步操作或分析(見表1.1)。定義了幾個(gè)訪問器函數(shù)來方便訪問樹注釋數(shù)據(jù)，包括get.field 用于獲取樹對(duì)象中可用的注釋特性的字段get.placements 用于獲得系統(tǒng)發(fā)育放置結(jié)果(即 PPLACER, EPA等的輸出)，get.subs用于獲取父節(jié)點(diǎn)到子節(jié)點(diǎn)的遺傳替換，get.tipseq用于得到 tip sequences。
S3類phlo是由 ape (Paradis, Claude, and Strimmer 2004)包中定義的，在R社區(qū)和許多包中被廣泛使用。由于treeio使用S4類，為了使那些可用的R包能夠分析treeio導(dǎo)入的樹，treeio提供了as.phylo函數(shù)將treeio生成的樹對(duì)象轉(zhuǎn)換為只包含樹結(jié)構(gòu)而不包含注釋數(shù)據(jù)的Phylo對(duì)象。另一方面，treeio也提供了as.treedata函數(shù)將phylo對(duì)象與進(jìn)化分析結(jié)果(例如,由ape計(jì)算的bootstrap values或 phangorn (Schliep 2011) 推斷的祖先狀態(tài)等)存儲(chǔ)treedata S4對(duì)象,因此很容易將數(shù)據(jù)映射到樹結(jié)構(gòu)和使用ggtree可視化 ggtree (Yu et al. 2017)。
為了在系統(tǒng)發(fā)育樹中集成不同類型的數(shù)據(jù)，treeio (Wang et al. 2020)提供了merge_tree函數(shù)(詳細(xì)信息見2.2.1節(jié))，用于結(jié)合從不同來源導(dǎo)入的進(jìn)化統(tǒng)計(jì)數(shù)據(jù)/證據(jù)，包括那些常見的樹文件和分析程序的輸出(表1.1)。還有其他信息，如采樣位置、分類信息、實(shí)驗(yàn)結(jié)果、進(jìn)化特征等，以用戶定義的格式存儲(chǔ)在單獨(dú)的文件中。在treeio中，我們可以使用標(biāo)準(zhǔn)的R IO函數(shù)從用戶文件中讀取這些數(shù)據(jù)，并通過在tidytree和treeio包中定義的full_join方法將它們附加到樹對(duì)象(請(qǐng)參閱ggtree中定義的 %<+% operator)。附加后的數(shù)據(jù)將成為節(jié)點(diǎn)或分支的關(guān)聯(lián)屬性，可以與合并的其他數(shù)據(jù)進(jìn)行比較，也可以在樹中直觀顯示 (Yu et al. 2018)。
為了便于存儲(chǔ)與系統(tǒng)發(fā)育樹相關(guān)的復(fù)雜數(shù)據(jù)，treeio應(yīng)用write.baset和write.jtree功能可以將一個(gè)treedata對(duì)象導(dǎo)出到單個(gè)文件中(參見第3章)。
1.3.2功能演示
#######1.3.2.1 BEAST輸出解析
file <- system.file("extdata/BEAST", "beast_mcc.tree", package="treeio")
beast <- read.beast(file)
beast
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/BEAST/beast_mcc.tree'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 14 internal nodes.
##
## Tip labels:
## A_1995, B_1996, C_1995, D_1987, E_1996, F_1997, ...
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'height', 'height_0.95_HPD', 'height_median',
## 'height_range', 'length', 'length_0.95_HPD',
## 'length_median', 'length_range', 'posterior', 'rate',
## 'rate_0.95_HPD', 'rate_median', 'rate_range'.
由于%不是names中的有效字符，所有包含x%的特性名都會(huì)轉(zhuǎn)換為0.x。例如，length_95%_HPD將被更改為length_0.95_HPD。
不僅樹結(jié)構(gòu)，而且BEAST推斷的所有特性都將存儲(chǔ)在S4對(duì)象中。這些特性可以用于樹的注釋(圖5.8)。
1.3.2.2 MEGA輸出解析
分子進(jìn)化遺傳學(xué)分析(MEGA)軟件(Kumar, Stecher，和Tamura 2016)支持導(dǎo)出三種不同格式的樹:Newick，tabular，和Nexus。Newick文件可以使用read.tree或read.newick來解析。newick功能。MEGA Nexus文件類似于BEAST Nexus, treeio (Wang et al. 2020)，可以使用read.mega功能來解析。
file <- system.file("extdata/MEGA7", "mtCDNA_timetree.nex",
package = "treeio")
read.mega(file)
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/MEGA7/mtCDNA_timetree.nex'.
##
## ...@ phylo:
##
## Phylogenetic tree with 7 tips and 6 internal nodes.
##
## Tip labels:
## homo_sapiens, chimpanzee, bonobo, gorilla,
## orangutan, sumatran, ...
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'branch_length', 'data_coverage', 'rate',
## 'reltime', 'reltime_0.95_CI', 'reltime_stderr'.
表格輸出包含表格平面文本文件中的樹和相關(guān)信息(本例中的分化時(shí)間)。ead.mega_tabular 功能可以同時(shí)用數(shù)據(jù)解析樹。
file <- system.file("extdata/MEGA7", "mtCDNA_timetree_tabular.txt",
package = "treeio")
read.mega_tabular(file)
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/MEGA7/mtCDNA_timetree_tabular.txt'.
##
## ...@ phylo:
##
## Phylogenetic tree with 7 tips and 6 internal nodes.
##
## Tip labels:
## chimpanzee, bonobo, homo sapiens, gorilla,
## orangutan, sumatran, ...
## Node labels:
## , , demoLabel2, , ,
##
## Rooted; no branch lengths.
##
## with the following features available:
## 'RelTime', 'CI_Lower', 'CI_Upper', 'Rate', 'Data
## Coverage'.
1.3.2.3解析MrBayes輸出
雖然MrBayes生成的Nexus文件與BEAST的輸出不同，但它們是相似的。treeio包提供了read.mrbayes()，它在內(nèi)部調(diào)用read.beast()來解析MrBayes的輸出。
file <- system.file("extdata/MrBayes", "Gq_nxs.tre", package="treeio")
read.mrbayes(file)
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/MrBayes/Gq_nxs.tre'.
##
## ...@ phylo:
##
## Phylogenetic tree with 12 tips and 10 internal nodes.
##
## Tip labels:
## B_h, B_s, G_d, G_k, G_q, G_s, ...
##
## Unrooted; includes branch lengths.
##
## with the following features available:
## 'length_0.95HPD', 'length_mean', 'length_median',
## 'prob', 'prob_range', 'prob_stddev', 'prob_percent',
## 'prob+-sd'.
1.3.2.4 PAML輸出解析
最大似然系統(tǒng)發(fā)育分析(PAML)是一套工具，用于DNA和蛋白質(zhì)序列的系統(tǒng)發(fā)育分析使用最大似然。樹搜索算法是在 BASEML和CODEML中實(shí)現(xiàn)的。treeio中提供的read.paml_rst()函數(shù)可以解析BASEML和CODEML中的rst文件。唯一的區(qū)別是序列中的空間。對(duì)于BASEML，每十個(gè)堿基用一個(gè)空格隔開，而對(duì)于CODEML，每三個(gè)堿基(三聯(lián)體)用一個(gè)空格隔開。
brstfile <- system.file("extdata/PAML_Baseml", "rst", package="treeio")
brst <- read.paml_rst(brstfile)
brst
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/PAML_Baseml/rst'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 13 internal nodes.
##
## Tip labels:
## A, B, C, D, E, F, ...
## Node labels:
## 16, 17, 18, 19, 20, 21, ...
##
## Unrooted; includes branch lengths.
##
## with the following features available:
## 'subs', 'AA_subs'.
類似地，我們可以從CODEML解析rst文件。
crstfile <- system.file("extdata/PAML_Codeml", "rst", package="treeio")
## type can be one of "Marginal" or "Joint"
crst <- read.paml_rst(crstfile, type = "Joint")
crst
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/PAML_Codeml/rst'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 13 internal nodes.
##
## Tip labels:
## A, B, C, D, E, F, ...
## Node labels:
## 16, 17, 18, 19, 20, 21, ...
##
## Unrooted; includes branch lengths.
##
## with the following features available:
## 'subs', 'AA_subs'.
樹節(jié)點(diǎn)。treeio (Wang et al. 2020)將自動(dòng)確定每個(gè)分支兩端序列之間的替換。氨基酸替換也將通過將核苷酸序列翻譯成氨基酸序列來確定。這些計(jì)算得到的替換也將存儲(chǔ)在S4對(duì)象中，以便以后進(jìn)行有效的樹注釋(圖5.10)。
CODEML推斷選擇壓力和估計(jì)dN/dS, dN和dS。這些信息存儲(chǔ)在輸出文件mlc中，該文件可由read.codeml_mlc()函數(shù)解析。
mlcfile <- system.file("extdata/PAML_Codeml", "mlc", package="treeio")
mlc <- read.codeml_mlc(mlcfile)
mlc
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/PAML_Codeml/mlc'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 13 internal nodes.
##
## Tip labels:
## A, B, C, D, E, F, ...
## Node labels:
## 16, 17, 18, 19, 20, 21, ...
##
## Unrooted; includes branch lengths.
##
## with the following features available:
## 't', 'N', 'S', 'dN_vs_dS', 'dN', 'dS', 'N_x_dN',
## 'S_x_dS'.
可以像前面演示的那樣分別解析rst和mlc文件，也可以使用read.codeml()函數(shù)將它們一起解析。
## tree can be one of "rst" or "mlc" to specify
## using tree from which file as base tree in the object
ml <- read.codeml(crstfile, mlcfile, tree = "mlc")
ml
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/PAML_Codeml/rst',
## '/home/ygc/R/library/treeio/extdata/PAML_Codeml/mlc'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 13 internal nodes.
##
## Tip labels:
## A, B, C, D, E, F, ...
## Node labels:
## 16, 17, 18, 19, 20, 21, ...
##
## Unrooted; includes branch lengths.
##
## with the following features available:
## 'subs', 'AA_subs', 't', 'N', 'S', 'dN_vs_dS', 'dN',
## 'dS', 'N_x_dN', 'S_x_dS'.
rst和mlc文件中的所有特性都被導(dǎo)入到一個(gè)單獨(dú)的S4對(duì)象中，因此可用于進(jìn)一步的注釋和可視化。例如，我們可以在同一棵系統(tǒng)發(fā)育樹上注釋和顯示dN/dS(來自mlc文件)和氨基酸替換(來自rst文件)(Yu et al. 2017)。
1.3.2.5 HyPhy輸出解析
HyPhy (Hypothesis testing using Phylogenies)是一個(gè)用于分析基因序列的軟件包。由HyPhy推斷出的祖先序列存儲(chǔ)在Nexus輸出文件中，其中包含樹形拓?fù)浜妥嫦刃蛄?。要解析這個(gè)數(shù)據(jù)文件，用戶可以使用read.hyphy.seq()函數(shù)。
ancseq <- system.file("extdata/HYPHY", "ancseq.nex", package="treeio")
read.hyphy.seq(ancseq)
## 13 DNA sequences in binary format stored in a list.
##
## All sequences of same length: 2148
##
## Labels:
## Node1
## Node2
## Node3
## Node4
## Node5
## Node12
## ...
##
## Base composition:
## a c g t
## 0.335 0.208 0.237 0.220
## (Total: 27.92 kb)
為了將序列映射到樹中，用戶還應(yīng)該提供一個(gè)內(nèi)部的節(jié)點(diǎn)標(biāo)記樹。如果用戶想要確定替換，他們還需要提供提示序列。在這種情況下，替換將自動(dòng)確定，就像我們解析 CODEML的輸出一樣。
nwk <- system.file("extdata/HYPHY", "labelledtree.tree", package="treeio")
tipfas <- system.file("extdata", "pa.fas", package="treeio")
hy <- read.hyphy(nwk, ancseq, tipfas)
hy
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/HYPHY/labelledtree.tree'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 13 internal nodes.
##
## Tip labels:
## K, N, D, L, J, G, ...
## Node labels:
## Node1, Node2, Node3, Node4, Node5, Node12, ...
##
## Unrooted; includes branch lengths.
##
## with the following features available:
## 'subs', 'AA_subs'.
1.3.2.6 r8輸出解析
r8s包使用參數(shù)、半?yún)?shù)和非參數(shù)方法松弛分子鐘，以便更好地估計(jì)散度時(shí)間和進(jìn)化速率 (Sanderson 2003)。日志文件中輸出三棵樹，時(shí)間樹、速率樹和絕對(duì)替換樹分別為TREE、RATO和PHYLO。
時(shí)間樹是由散度時(shí)間標(biāo)度的，速率樹是由替換率標(biāo)度的，絕對(duì)替換樹是由替換的絕對(duì)數(shù)量標(biāo)度的。解析文件后，這三棵樹都存儲(chǔ)在一個(gè)multiPhylo對(duì)象中(圖4.15)。
r8s <- read.r8s(system.file("extdata/r8s", "H3_r8s_output.log", package="treeio"))
r8s
## 3 phylogenetic trees
1.3.2.7 RAxML引導(dǎo)分析的解析輸出
RAxML引導(dǎo)分析輸出的Newick樹文本不是標(biāo)準(zhǔn)的，因?yàn)樗鼘⒁龑?dǎo)值存儲(chǔ)在分支長度后的方括號(hào)中。這個(gè)文件通常不能被傳統(tǒng)的Newick解析器解析，比如ape::read.tree()。函數(shù)read.raxml()可以讀取這樣的文件，并將引導(dǎo)程序存儲(chǔ)為一個(gè)額外的特性，它可以用于顯示在樹中，或用于為樹分支著色，等等。
raxml_file <- system.file("extdata/RAxML",
"RAxML_bipartitionsBranchLabels.H3",
package="treeio")
raxml <- read.raxml(raxml_file)
raxml
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/RAxML/RAxML_bipartitionsBranchLabels.H3'.
##
## ...@ phylo:
##
## Phylogenetic tree with 64 tips and 62 internal nodes.
##
## Tip labels:
## A_Hokkaido_M1_2014_H3N2_2014,
## A_Czech_Republic_1_2014_H3N2_2014,
## FJ532080_A_California_09_2008_H3N2_2008,
## EU199359_A_Pennsylvania_05_2007_H3N2_2007,
## EU857080_A_Hong_Kong_CUHK69904_2006_H3N2_2006,
## EU857082_A_Hong_Kong_CUHK7047_2005_H3N2_2005, ...
##
## Unrooted; includes branch lengths.
##
## with the following features available:
## 'bootstrap'.
1.3.2.8解析NHX樹
NHX (New Hampshire eXtended)格式是通過引入NHX標(biāo)簽對(duì)Newick的擴(kuò)展。NHX通常用于系統(tǒng)發(fā)育軟件，包括 PHYLDOG (Boussau et al. 2013), RevBayes (H?hna et al. 2014)，用于存儲(chǔ)統(tǒng)計(jì)推斷。下面的代碼導(dǎo)入了一個(gè)帶有由PHYLDOG推斷的相關(guān)數(shù)據(jù)的NHX樹(圖3.1A)。
nhxfile <- system.file("extdata/NHX", "phyldog.nhx", package="treeio")
nhx <- read.nhx(nhxfile)
nhx
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/NHX/phyldog.nhx'.
##
## ...@ phylo:
##
## Phylogenetic tree with 16 tips and 15 internal nodes.
##
## Tip labels:
## Prayidae_D27SS7@2825365, Kephyes_ovata@2606431,
## Chuniphyes_multidentata@1277217,
## Apolemia_sp_@1353964, Bargmannia_amoena@263997,
## Bargmannia_elongata@946788, ...
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'Ev', 'ND', 'S'.
1.3.2.9解析Phylip tree
Phylip格式包含多個(gè)序列對(duì)齊的Phylip序列格式與相應(yīng)的Newick樹文本，由分類單元序列構(gòu)建。多序列比對(duì)可通過ggtree的msaplot()函數(shù)或結(jié)合ggmsa包進(jìn)行樹狀結(jié)構(gòu)排序，并顯示在樹的右側(cè)(see also Basic Protocol 5 of (Yu 2020))。
phyfile <- system.file("extdata", "sample.phy", package="treeio")
phylip <- read.phylip(phyfile)
phylip
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/sample.phy'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 13 internal nodes.
##
## Tip labels:
## K, N, D, L, J, G, ...
##
## Unrooted; no branch lengths.
1.3.2.10解析EPA和pplacer輸出
EPA (Berger, Krompass, and Stamatakis 2011)和PPLACER (Frederick A. Matsen, Kodner, and Armbrust 2010a)具有共同的輸出文件格式j(luò)place，可以通過read.jplace()函數(shù)進(jìn)行解析。
jpf <- system.file("extdata/EPA.jplace", package="treeio")
jp <- read.jplace(jpf)
print(jp)
## 'treedata' S4 object that stored information
## of
## '/home/ygc/R/library/treeio/extdata/EPA.jplace'.
##
## ...@ phylo:
##
## Phylogenetic tree with 493 tips and 492 internal nodes.
##
## Tip labels:
## CIR000447A, CIR000479, CIR000078, CIR000083,
## CIR000070, CIR000060, ...
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'nplace'.
每個(gè)分支上的進(jìn)化放置數(shù)量將被計(jì)算并存儲(chǔ)為nplace特征，可以使用ggtree將其映射為行大小和/或顏色(Yu et al. 2017)。
1.3.2.11解析jtree格式
jtree是一種基于json的格式，在treeio包中定義(Wang et al. 2020)，以支持樹數(shù)據(jù)交換(見 session 3.3)。帶有相關(guān)數(shù)據(jù)的系統(tǒng)發(fā)育樹可以使用write.jtree()函數(shù)導(dǎo)出到單個(gè)jtree文件。jtree可以使用任何JSON解析器輕松解析。jtree格式包含三塊:tree、data和metadata。樹值包含從Newick樹格式擴(kuò)展而來的樹文本，方法是將邊號(hào)放在分支長度之后的花括號(hào)中。數(shù)據(jù)值包含特定于節(jié)點(diǎn)/分支的數(shù)據(jù)，而元數(shù)據(jù)值包含額外的元信息。
jtree_file <- tempfile(fileext = '.jtree')
write.jtree(beast, file = jtree_file)
read.jtree(file = jtree_file)
## 'treedata' S4 object that stored information
## of
## '/tmp/Rtmp5zDeZ2/file1eea82133a2fb.jtree'.
##
## ...@ phylo:
##
## Phylogenetic tree with 15 tips and 14 internal nodes.
##
## Tip labels:
## K_2013, N_2010, D_1987, L_1980, J_1983, G_1992, ...
##
## Rooted; includes branch lengths.
##
## with the following features available:
## 'height', 'height_0.95_HPD', 'height_range',
## 'length', 'length_0.95_HPD', 'length_median',
## 'length_range', 'rate', 'rate_0.95_HPD',
## 'rate_median', 'rate_range', 'height_median',
## 'posterior'.
1.3.3將其他樹狀對(duì)象轉(zhuǎn)換為phylo或treedata
為了擴(kuò)展treeio、tidytree和ggtree的應(yīng)用范圍，treeio (Wang et al. 2020)提供了幾個(gè)as.phylo和as.treedata方法，用于將其他樹狀對(duì)象(如phylo4d和pml)轉(zhuǎn)換為phylo或treedata對(duì)象。因此，用戶可以很容易地將相關(guān)數(shù)據(jù)映射到樹結(jié)構(gòu)，將有/沒有數(shù)據(jù)的樹導(dǎo)出到單個(gè)文件，操作和可視化有/沒有數(shù)據(jù)的樹。這些轉(zhuǎn)換函數(shù)(表1.2)創(chuàng)造了使用tidy接口來處理tree和使用圖形語法的ggtree來可視化tree的可能性。
|-|-|-|
| Convert function | Supported object | Description |
|---|---|---|
| as.phylo | ggtree | convert ggtree object to phylo object |
| as.phylo | igraph | convert igraph object (only tree graph supported) to phylo object |
| as.phylo | phylo4 | convert phylo4 object to phylo object |
| as.phylo | pvclust | convert pvclust object to phylo object |
| as.phylo | treedata | convert treedata object to phylo object |
| as.treedata | ggtree | convert ggtree object to treedata object |
| as.treedata | phylo4 | convert phylo4 object to treedata object |
| as.treedata | phylo4d | convert phylo4d object to treedata object |
| as.treedata | pml | convert pml object to treedata object |
| as.treedata | pvclust | convert pvclust object to treedata object |
表1.2:樹狀對(duì)象到phylo或 treedata對(duì)象的轉(zhuǎn)換
這里，我們以pml對(duì)象為例，它是在phangorn包中定義的。pml()函數(shù)計(jì)算給定序列比對(duì)和模型的系統(tǒng)發(fā)育樹的可能性，而 optim.pml()函數(shù)優(yōu)化不同的模型參數(shù)。輸出是一個(gè)pml對(duì)象，可以使用 treeio (Wang et al. 2020)的as.treedata將其轉(zhuǎn)換為一個(gè)treedata對(duì)象。存儲(chǔ)在樹狀數(shù)據(jù)對(duì)象中的氨基酸替換(pml估計(jì)的祖先序列)可以使用ggtree進(jìn)行可視化，如圖1.4所示。
library(phangorn)
treefile <- system.file("extdata", "pa.nwk", package="treeio")
tre <- read.tree(treefile)
tipseqfile <- system.file("extdata", "pa.fas", package="treeio")
tipseq <- read.phyDat(tipseqfile,format="fasta")
fit <- pml(tre, tipseq, k=4)
fit <- optim.pml(fit, optNni=FALSE, optBf=T, optQ=T,
optInv=T, optGamma=T, optEdge=TRUE,
optRooted=FALSE, model = "GTR",
control = pml.control(trace =0))
pmltree <- as.treedata(fit)
ggtree(pmltree) + geom_text(aes(x=branch, label=AA_subs, vjust=-.5))
圖1.4:將pml對(duì)象轉(zhuǎn)換為treedata對(duì)象這允許使用tidytree來處理樹數(shù)據(jù)，以及使用ggtree和ggtreeExtra來顯示與相關(guān)數(shù)據(jù)相關(guān)的樹。
1.3.4從樹數(shù)據(jù)對(duì)象中獲取信息
導(dǎo)入樹后，用戶可能希望提取存儲(chǔ)在樹數(shù)據(jù)對(duì)象中的信息。Treeio提供了幾種訪問器方法來提取樹結(jié)構(gòu)、存儲(chǔ)在對(duì)象中的特性/屬性以及它們對(duì)應(yīng)的值。
get.tree()或as.phylo()方法可以將treedata對(duì)象轉(zhuǎn)換為一個(gè)phylo對(duì)象，這是R社區(qū)中的基本樹對(duì)象，許多包都使用phylo對(duì)象。
beast_file <- system.file("examples/MCC_FluA_H3.tree", package="ggtree")
beast_tree <- read.beast(beast_file)
# or get.tree
as.phylo(beast_tree)
##
## Phylogenetic tree with 76 tips and 75 internal nodes.
##
## Tip labels:
## A/Hokkaido/30-1-a/2013, A/New_York/334/2004, A/New_York/463/2005, A/New_York/452/1999, A/New_York/238/2005, A/New_York/523/1998, ...
##
## Rooted; includes branch lengths.
get.fields方法返回存儲(chǔ)在對(duì)象中并與系統(tǒng)發(fā)育相關(guān)的特征/屬性向量。
get.fields(beast_tree)
## [1] "height" "height_0.95_HPD"
## [3] "height_median" "height_range"
## [5] "length" "length_0.95_HPD"
## [7] "length_median" "length_range"
## [9] "posterior" "rate"
## [11] "rate_0.95_HPD" "rate_median"
## [13] "rate_range"
get.data 方法返回所有關(guān)聯(lián)數(shù)據(jù)的浮點(diǎn)數(shù)。
get.data(beast_tree)
## # A tibble: 151 × 14
## height height_0.95_HPD height_median height_range
## <dbl> <list> <dbl> <list>
## 1 19 <dbl [2]> 19 <dbl [2]>
## 2 17 <dbl [2]> 17 <dbl [2]>
## 3 14 <dbl [2]> 14 <dbl [2]>
## 4 12 <dbl [2]> 12 <dbl [2]>
## 5 9 <dbl [2]> 9 <dbl [2]>
## 6 10 <dbl [2]> 10 <dbl [2]>
## 7 10 <dbl [2]> 10 <dbl [2]>
## 8 10.8 <dbl [2]> 10.8 <dbl [2]>
## 9 9 <dbl [2]> 9 <dbl [2]>
## 10 9 <dbl [2]> 9 <dbl [2]>
## # … with 141 more rows, and 10 more variables:
## # length <dbl>, length_0.95_HPD <list>,
## # length_median <dbl>, length_range <list>,
## # posterior <dbl>, rate <dbl>, rate_0.95_HPD <list>,
## # rate_median <dbl>, rate_range <list>, node <int>
如果用戶只想使用get.fields獲得features/attributes的部分子集數(shù)據(jù)，可以使用get.data或直接使用 [ 或 [[ 提取子集。
beast_tree[, c("node", "height")]
## # A tibble: 151 × 2
## node height
## <int> <dbl>
## 1 10 19
## 2 9 17
## 3 36 14
## 4 31 12
## 5 29 9
## 6 28 10
## 7 39 10
## 8 90 10.8
## 9 16 9
## 10 2 9
## # … with 141 more rows
head(beast_tree[["height_median"]])
## height_median1 height_median2 height_median3
## 19 17 14
## height_median4 height_median5 height_median6
## 12 9 10
1.4總結(jié)
用于推斷分子進(jìn)化的軟件工具(例如，祖先狀態(tài)、分子年代測定和選擇壓力等)正在激增，但沒有一種單一的數(shù)據(jù)格式可以被所有不同的程序使用，并能夠存儲(chǔ)不同類型的系統(tǒng)發(fā)育數(shù)據(jù)。大多數(shù)軟件包都有其獨(dú)特的輸出格式，并且這些格式彼此之間不兼容。解析軟件輸出具有挑戰(zhàn)性，這限制了使用不同工具進(jìn)行聯(lián)合分析。treeio包(Wang et al. 2020)提供了一組函數(shù)(表1.1)，用于解析各種類型的系統(tǒng)發(fā)育數(shù)據(jù)文件，以及一組轉(zhuǎn)換器(表1.2)，用于將樹狀對(duì)象轉(zhuǎn)換為門或樹狀數(shù)據(jù)對(duì)象。這些系統(tǒng)發(fā)育數(shù)據(jù)可以被整合，從而允許進(jìn)一步的探索和比較。到目前為止，分子進(jìn)化領(lǐng)域的大多數(shù)軟件工具都是孤立的，并且常常不能完全兼容彼此的輸入和輸出文件。這些軟件工具被設(shè)計(jì)用來進(jìn)行分析，而輸出通常在其他軟件中是不可讀的。目前還沒有設(shè)計(jì)任何工具來統(tǒng)一來自不同分析程序的推理數(shù)據(jù)。有效整合不同推理方法的數(shù)據(jù)，可以增強(qiáng)對(duì)研究對(duì)象的比較和理解，有助于發(fā)現(xiàn)新的系統(tǒng)模式，產(chǎn)生新的假設(shè)。
隨著系統(tǒng)發(fā)育樹在進(jìn)化背景下識(shí)別模式的應(yīng)用越來越廣泛，越來越多的不同學(xué)科在研究中使用系統(tǒng)發(fā)育樹。例如，空間生態(tài)學(xué)家可以將生物體的地理位置繪制到系統(tǒng)發(fā)育樹中，以了解該物種的生物地理學(xué) (Sch?n et al. 2015);疾病流行病學(xué)家可以將病原體采樣的時(shí)間和地點(diǎn)納入系統(tǒng)發(fā)育分析中，推斷疾病在時(shí)空空間的傳播動(dòng)態(tài) (Y.-Q. He et al. 2013);微生物學(xué)家可以確定不同病原菌株的致病性，并將其繪制到系統(tǒng)發(fā)育樹中，以確定致病性的遺傳決定因素(Bosi et al. 2016);基因組科學(xué)家可以使用系統(tǒng)發(fā)育樹來幫助對(duì)他們的宏基因組序列數(shù)據(jù)進(jìn)行分類((Gupta and Sharma 2015)。treeio等強(qiáng)大的工具可以將不同類型的數(shù)據(jù)導(dǎo)入并映射到系統(tǒng)發(fā)育樹中，這對(duì)于促進(jìn)這些與系統(tǒng)發(fā)育相關(guān)(也被稱為phylodynamics）的研究有著重要意義。這些工具還可以幫助在最高水平上整合不同的元數(shù)據(jù)(時(shí)間、地理、基因型、流行病學(xué)信息)和分析結(jié)果(選擇壓力、進(jìn)化速率)，并提供對(duì)研究生物體的全面了解。在流感研究領(lǐng)域，通過在同一種系統(tǒng)發(fā)育樹和進(jìn)化時(shí)間尺度上繪制不同元數(shù)據(jù)和分析結(jié)果來研究流感病毒的phylodynamics(Lam et al. 2015)。
