1.3 使用treeio导入树及相关数据_R实战：系统发育树的数据集成操作及可视化-QQ阅读男生轻小说网

书名：R实战：系统发育树的数据集成操作及可视化
作者名：余光创
本章字数：5186字
更新时间：2023-12-12 18:56:23

1.3 使用treeio导入树及相关数据

我们经常使用进化树来展现物种间的进化关系，而与树中的物种、菌株等分类单元关联的信息则可以在进化树提供的进化历史背景下进行进一步的分析。例如，我们可以通过研究树中不同流感病毒株的宿主信息来探究某一病毒谱系的宿主范围。除此之外，诸如分离株宿主、感染时间、感染地点等与病毒株直接相关的元数据也常被应用于系统发育的进化分析及比较分析的建模中，以推断它们在病毒进化或传播过程中的动力学关系。我们通常将这类的元数据、表型数据或实验数据存储为与节点或分支关联的注释数据，使用不同的软件来存储此类树+注释的文件格式往往是不一致的。

目前来说，将进化树导入R还存在着许多限制。Newick和NEXUS格式尚可以被ape、phylobase、NeXML等多种包导入，NeXML格式可以使用RNeXML包来解析。系统发育分析领域中其他的一些常用软件的分析结果却没有得到很好的支持。例如，PHYLOCH包虽然支持导入BEAST及MrBayes的输出结果，但是只能解析内部节点属性，而忽略叶节点属性[4]。还有很多软件所输出的进化树及关联数据仍需要熟悉编程的人员才能将其导入R。除此之外，将一些外部数据（如实验数据及临床数据等）与进化树关联起来也是一件难事。

为了填补R［28］社区在这方面的空白，我们基于R语言开发了一款名为treeio［29］的R包，而它能解析大多数常用进化分析软件输出的不同格式的进化树文件。treeio能将进化树的拓扑结构及其关联的多种数据和进化推论一并解析，如NHX注释，BEAST或r8s［25］的分子钟速率推断，CODEML推断的同义替换与非同义替换，BASEML、HyPhy及CODEML重构的祖先序列等。我们在treeio［29］中开发了多种解析函数（见表1.1）来支持读取这些文件格式，目前已经支持读取的文件格式包括Newick、NEXUS、NHX、Jplace及Phylip，也支持读取ASTRAL、BEAST、EPA、HyPhy、MEGA、MrBayes、PAML、PHYLDOG、PPLACER、r8s、RAxML、RevBayes等多种软件的输出结果。

表1.1　解析函数名称及其说明

续表

treeio包为R中进化树的输入与输出定义了多种基本的类与函数。我们可以利用这些类与函数，将各类不同软件所推断的进化证据解析并导入R中。例如，CODEML［26］推断出的祖先序列或计算出的dn/ds值、BEAST［21］推断出的进化枝支持度（后验概率）、EPA［19］和PPLACER［18］推断出的短读序列放置信息等，这些进化证据都可以通过treeio导入R以供后续的各类分析并通过ggtree包［30］实现对进化树可视化上的注释。随着各类分析工具及模型的发展，将多种对相同进化树在不同软件分析中的结果整合起来共同分析也变得越来越困难。为了解决这个问题，treeio包［29］设计了merge_tree（）函数来整合不同来源的树数据。除此之外，treeio包还支持将外部数据与进化树关联到一起。

在解析完树后，treeio包使用在tidytree包中定义的S4类——treedata用来存储进化树及其相关数据。这些相关数据在treedata中会被映射到分支或节点上，以便后续能更有效地使用ggtree［30］和ggtreeExtra［31］来对树及其注释信息进行可视化。像这样一个专门为解析、整合及注释系统发育学数据而设计的可编程平台可以让我们更加容易地发现进化动力学关系及相关性模式。

图1.3所示为treeio包［29］概述及其与tidytree和ggtree的关系。treeio包用于解析多种文件格式、软件输出文件中的进化树及关联数据。treeio包使用treedata对象来存储系统发育树及其分支、节点关联的数据，也提供了多种函数来操作带有数据的树文件。用户可以将treedata对象转换为整洁数据框（tidy data frame，其中每一行表示树中的一个节点，每一列表示一个变量），从而使用tidytree中的整洁接口（tidy interface）来对树进行处理；同时，可以选择将树的拓扑结构从treedata中提取出来，存储到Newick文件或Nexus文件中，或者将树及关联数据一起导入同一个文件中（BEAST、Nexus或jtree），这些树关联的数据可以在可视化中使用ggtree对树进行注释。除了常见进化树的可视化，ggtree也支持其他树形对象的可视化，包括处理宏基因组学数据的phyloseq对象及处理病毒爆发数据的obkData对象。ggtree还支持对phylo、multiPhylo（ape包）、phylo4、phylo4d（phylobase包）、phylog（ade4包）、phyloseq（phyloseq包）obkData（OutbreakTools包）等多种R社区中所定义的用以存储树与特定领域数据的树对象的可视化，以及层次聚类结果（如hclust与dendrogram对象）的可视化。

图1.3　treeio包概述及其与tidytree和ggtree的关系

1.3.1 treeio简介

treeio包［29］针对进化树的输入和输出进行设计，定义了S4类对象treedata来存储进化树及各种软件包输出在内的不同种类的关联数据或协变量，同时设计了相应的解析器来将对应的进化树及其注释数据转换为R的对象，以便于后续的数据处理及分析。除此之外，treeio还设计了多种存取器（accessor），便于提取树文件中的注释信息。例如，get.fields（）函数用于识别并提取树中所能获取到的注释信息；get.placements（）函数用于提取PPLACER、EPA软件所得出的进化放置信息（phylogenetic placement）；get.subs（）函数用于提取从父节点到子节点的遗传替换信息；get.tips（）函数用于提取叶节点的序列信息。

目前，R社区内大部分进化树相关的R包中使用的仍是ape包［32］中定义的S3类对象phylo。为了利用treeio导入树后，能继续使用这些包进行分析，treeio中的as.phylo（）函数并将其生成的S4类对象treedata转换为S3类对象phylo。但需要注意的是，转换后的phylo对象只含有进化树的拓扑结构信息，而不包含任何注释信息。treeio中的as.treedata（）函数用于将含有进化分析结果（如利用ape计算出的自举值或利用phangorn［33］推断的祖先状态等）的phylo对象转换为S4类对象treedata。转换之后，用户可以更加便利地把注释信息映射到树结构上，以便于后续使用ggtree［30］实现对树的可视化。

为了便于将多种不同种类的数据整合到进化树，在treeio［29］中定义了merge_tree（）函数。而对于存储在用户自定义文件类型中的其他信息（如采样地点、分类学信息、实验结果、演化性状等进化证据），先通过treeio及tidytree中定义的full_join（）函数对这些数据进行整合，再将这些文件通过R中的I/O函数读取后，我们便可以通过full_join（）函数（或者使用ggtree中的％＜+％操作符）将这些信息与树关联起来，并将它们变为分支或节点的属性，以便于后续与其他数据进行比对，或者在可视化时将其展示于树上［34］。

同时，为了更好地存储树及其关联的复杂数据，treeio还提供了write.beast（）函数与write.jtree（）函数，用于将treedata对象输出至单一文件中。

1.3.2 treeio解析函数演示

1.3.2.1　解析BEAST输出文件

由于％在names中并不是一个合法的字符，treeio会自动将特征名中的百分数（x％）转换为小数形式（0.x），如将length_95％_HPD转换为length_0.95_HPD。

同时，除了树结构本身，treeio还能将所有BEAST软件的分析结果一起存储到S4对象中，以供在后续流程中提供对树的注释。

1.3.2.2　解析MEGA输出文件

MEGA（Molecular Evolutionary Genetics Analysis）［22］是一款可以进行序列比对及进化树构建的软件。在构建进化树后，MEGA支持将进化树输出为Newick、Tabular和Nexus这3种格式的文件。其中，treeio中的read.tree（）函数或read.newick（）函数用于解析Newick文件；MEGA与BEAST输出的都是非标准的Nexus文件。treeio［29］中的read.mega（）函数可用于解析MEGA输出的Nexus文件。

我们利用MEGA的Tabular输出文件将树及其关联信息（本示例中为分歧时间）存储在一个表格化的平面文本文件中。treeio中的read.mega_tubular（）函数用于解析MEGA输出的Tabular文件。

1.3.2.3　解析MrBayes输出文件

MrBayes软件输出的同样是非标准Nexus文件，与BEAST的输出文件稍有不同但大致相似。因此，在使用treeio中的read.mrbayes（）函数解析MrBayes输出文件时，会在内部调用read.beast（）函数。

1.3.2.4　解析PAML输出文件

PAML（Phylogenetic Analysis by Maximum Likelihood）是一个使用最大似然法来对DNA和蛋白质序列进行系统发育分析的软件包，在其子程序BASEML与CODEML中实现了多种树搜索算法。treeio中的read.paml_rst（）函数用于解析BASEML与CODEML输出的rst文件。两者生成的rst文件的主要区别在于序列间的空格位置不同。在BASEML生成的rst文件中，每10个碱基间会由一个空格隔开，而在CODEML生成的rst文件中则是每3个碱基（一个三联体）间会由一个空格隔开。

也可以通过read.paml_rst（）函数来解析CODEML中的rst文件。

BASEML与CODEML通过边缘最大似然法或联合最大似然法重建的祖先序列会被存储在S4类对象treedata中，并映射到对应的节点上。treeio［29］会自动确定每个分支两端节点上的序列间的替换信息。我们通过将核苷酸序列翻译为氨基酸序列，也可以确定氨基酸序列间的替换信息。这些经由计算的替换信息会被存储在treedata对象中，以便于后续对进化树进行注释与可视化。

CODEML可以推断选择压力，并估算dN/dS、dN及dS的值。这些信息会被存储于CODEML输出的mlc文件中，并通过read.codeml_mlc（）函数来解析。

treeio不仅能单独解析rst文件与mlc文件，而且能使用read.codeml（）函数同时读取rst文件和mlc文件。

所有rst文件与mlc文件中的特征数据都会被导入单一S4对象中，以便后续对进化树进行注释与可视化。例如，我们可以在同一棵进化树上［30］同时展示mlc文件中获取的dN/dS值及rst文件中获取的氨基酸替换信息。

1.3.2.5　解析HyPhy输出文件

HyPhy（Hypothesis testing using Phylogenies）是一个用于进行遗传序列分析的软件包。HyPhy会将推断的祖先序列与树的拓扑结构一起输出到Nexus文件中，并使用treeio中的read.hyphy.seq（）函数解析此文件。

如果想要将序列映射至树上，则需要提供树内部节点的标签；如果想要确定替换信息，则需要提供叶节点的序列信息。在本示例中，与解析CODEML输出文件时类似，替换信息会被treeio自动确定。

1.3.2.6　解析r8s输出文件

r8s包可以通过参数、半参数及非参数的方法来设定宽松分子钟，以便更加精准地估算分歧时间及进化速率［25］。在输出时，r8s会将3棵进化树存储到一个日志文件中。这3棵进化树分别为时间树（TREE）、速率树（RATO）与绝对替换树（PHYLO）。

时间树使用分歧时间作为树枝长的标尺，速率树使用替换速率作为树枝长的标尺，绝对替换树以替换的绝对值作为树枝长的标尺。在解析完文件后，treeio会将这3棵进化树存储到multiPhylo对象中。

1.3.2.7　解析RAxML自举分析输出文件

RAxML自举分析软件输出的是非标准的Newick文件，其在枝长后还添加了方括号来存储自举值。这导致此文件一般不能用常规的Newick解析器来解析，如ape包中的read.tree（）函数等。而treeio中的read.raxml（）函数则用于解析RAxML输出文件，利用该函数提取出自举值，存储为树的另一特征，方便后续进行可视化时将自举值显示于树上，或者根据自举值来为树的分支着色等。

1.3.2.8　解析NHX树

NHX是一种Newick的拓展格式，由于在Newick的基础上引进了NHX标签，因此NHX文件能存储更多信息。NHX文件在系统发育学的软件中比较常见，如PHYLDOG［14］与RevBayes［35］就是使用NHX格式来存储统计推论的。而treeio中的read.nhx（）函数则用于解析NHX格式的文件。下面的代码展示了如何使用treeio导入含有PHYLDOG推断出的关联数据的NHX树。

1.3.2.9　解析Phylip树

Phylip格式内包含了Phylip序列格式的分类单元多序列比对信息，以及根据这些分类单元序列所建的Newick树。treeio中的read.phylip（）函数用于解析Phylip格式的文件。ggtree中的msaplot（）函数用于将多序列比对结果根据树的结构重新排列，并将其展示于树的右侧。除此之外，我们也可以结合ggmsa包来对多序列比对结果进行分析（另请参见本章参考文献［36］中的基本流程5）。

1.3.2.10　解析EPA和pplacer的输出文件

EPA［19］与pplacer［18］输出的均为Jplace格式文件，而treeio中的read.jplace（）函数则用于解析Jplace格式的文件。

利用treeio计算每个分支上进化放置信息的数量，并将其存储于nplace特征中，以便在后续使用ggtree［30］进行可视化时用来设置分支线条的粗细或颜色。

1.3.2.11　解析jtree文件

jtree格式是在treeio包［29］中定义的，是基于JSON格式的文件，用于支持树数据在不同编程语言之间的流通。我们可以使用wirte.jtree（）函数来将进化树及其关联数据输出到同一个jtree文件中，所有支持JSON的解析器都可以解析jtree文件。其包含3个键（key）：tree、data及metadata。其中，tree的值包含由Newick格式拓展而来的树文本，在枝长后还添加了花括号，用于存储边编号；data的值包含了与分支或节点关联的数据；metadata的值包含了额外的一些元信息。

1.3.3 将其他树形对象转换为phylo对象或treedata对象

为了拓展treeio［29］、tidytree和ggtree的应用范围，treeio为as.phylo和as.treedata定义了多种函数，以将类似Phylo4d、pml等树形对象转换为phylo对象或treedata对象。这样，用户就可以通过treeio包来更加便利地实现一些操作。例如，将相关数据映射至树结构上、将含有或不含有相关数据的树输出到同一个文件中、对带有数据的树进行可视化。这些转换函数让我们可以使用tidytree对整洁接口进行数据整理，或者使用ggtree利用图形语法（grammar of graphic）对树进行可视化。表1.2所示为将树形对象转换为phylo对象或treedata对象的函数。

表1.2　将树形对象转换为phylo对象或treedata对象的函数

续表

这里使用pml对象来举例说明。pml对象是在phangorn包中定义的，其中，pml（）函数用于根据给定的序列比对信息及模型计算出进化树的似然性，计算完成后可以使用optm.pml（）函数对模型参数进行优化。最后输出一个pml对象，通过treeio［29］中的as.treedata（）函数来将其转换为treedata对象，而pml对象中的氨基酸替换信息（pml推断出的祖先序列）也会被存储到treedata对象中，以供后续分析中使用ggtree进行可视化。如图1.4所示，将pml对象转换为treedata对象后就可以使用tidytree来处理树及其相关数据，也可以使用ggtree和ggtreeExtre来对树及其关联数据进行可视化。