ixxmu
commented
5 months ago 使用 R+ggtree 可视化和注释系统发育树 by 生信七点半
本课程展示了如何使用ggtree，这是ggplot2包的扩展，用于可视化和注释系统发育树。
ggtree包
ggtree是一个R包，扩展了ggplot2以用于可视化和注释系统发育树及其共变量和其他相关数据。它可从Bioconductor获得。Bioconductor是一个项目，提供用于分析和注释各种基因组数据的工具。您可以在这里搜索和浏览Bioconductor包。
ggtreeBioconductor页面: bioconductor.org/packages/ggtree。ggtree主页: guangchuangyu.github.io/ggtree（包含有关包的更多信息、更多文档、美丽图像的图库，以及相关资源的链接）。ggtree出版物: Yu, Guangchuang, et al. “ggtree: an r package for visualization and annotation of phylogenetic trees with their covariates and other associated data.” Methods in Ecology and Evolution (2016) DOI:10.1111/2041-210X.12628。
Bioconductor包通常有很好的文档，以使用说明的形式呈现。查看ggtree的登录页面——在页面中部“文档”标题下有多个针对ggtree不同应用和功能的教程，充满了可运行的示例和解释。
就像从CRAN安装的R包一样，您只需要安装一次Bioconductor包（安装说明在这里），然后在每次开始新的R会话时加载它们。首先加载tidyverse包。
library(tidyverse)
library(ggtree)
关于被屏蔽的函数的说明：加载ggtree后，浏览一下您看到的一些输出。首次安装ggtree可能需要一段时间，因为ggtree依赖许多其他R包。这些包又可能依赖其他包。当您加载ggtree时，这些都会被加载到您的工作环境中。还请注意以“The following objects are masked from 'package:....”开头的行。这表示ggtree加载了它自己的某些函数，如dplyr的collapse()函数。现在，如果您想使用dplyr的collapse函数，您必须使用这种语法显式调用它：dplyr::collapse()。
导入树
ggtree是一个在R中用于可视化和注释系统发育树及其相关数据的包，作为ggplot2的扩展。它支持多种文件格式，这允许用户从不同的系统发育分析软件中导入树数据。以下是ggtree支持的一些主要文件格式及其对应的解析函数：
支持的文件格式
Newick: 一种广泛使用的表示树结构的文本格式。 Nexus: 一种包含丰富元数据的复杂格式，常用于进化分析。 Phylip: 另一种用于进化分析数据的格式。 Jplace: 用于存储植物放线菌类基因定位分析结果的文件格式。 NHX (New Hampshire eXtended format): Newick格式的扩展，允许在树中包含额外的注释信息。 BEAST, EPA, HYPHY, PAML, PHYLDOG, pplacer, r8s, RAxML, RevBayes: ggtree还支持这些系统发育分析软件的输出格式。
解析函数
read.tree: 用于读取Newick文件。 read.phylip: 用于读取Phylip文件。 read.jplace: 用于读取Jplace文件，包括来自EPA和pplacer的输出。 read.nhx: 用于读取NHX文件，包括来自PHYLODOG和RevBayes的输出。 read.beast: 用于解析BEAST软件的输出。 read.codeml 和 read.codeml_mlc: 用于解析CODEML软件的rst和mlc文件的输出。 read.hyphy: 用于解析HYPHY软件的输出。 read.paml_rst: 用于解析BASEML和CODEML软件的rst文件的输出。 read.r8s: 用于解析r8s软件的输出。 read.raxml: 用于解析RAxML软件的输出。
这些解析函数使ggtree能够处理来自不同来源和格式的系统发育树数据，进而允许用户在R中对这些树进行灵活的可视化和注释。ggtree的这些功能特别适合那些需要在系统发育分析中整合多种数据类型和来源信息的研究人员。利用ggtree，用户可以高效地探索系统发育树的结构，标注感兴趣的节点，以及将外部数据映射到树上，从而在生物学研究中获得深入的见解。本节演示了如何使用ggtree包来导入和可视化系统发育树数据。ggtree是基于ggplot2的一个扩展，专门用于系统发育树的可视化和注释。这里首先从导入系统发育树数据开始，展示了如何使用ggtree中的函数来构建和自定义树的基础图表。
导入系统发育树数据
library(tidyverse)
library(ggtree)
tree <- read.tree("data/tree_newick.nwk")
tree
这段代码首先加载了tidyverse和ggtree包，然后使用read.tree()函数从Newick格式的文件中导入系统发育树数据。tree对象包含了树的结构信息，但直接展示该对象通常不会提供有用的视觉信息。
基础树可视化
# 使用ggplot构建基础图形并添加geom_tree图层
ggplot(tree) + geom_tree() + theme_tree()
# 使用ggtree函数的便捷方式
ggtree(tree)
就像在ggplot2中使用ggplot()函数创建基础画布并通过+geom_???()添加图层一样，在ggtree中，我们可以使用geom_tree()函数添加树图层。ggtree()函数是为了方便使用而设计的，它相当于ggplot(...) + geom_tree() + theme_tree()的快捷方式。
添加比例尺
# 添加比例尺
ggtree(tree) + geom_treescale()
# 或使用theme_tree2()添加x轴的整体比例尺
ggtree(tree) + theme_tree2()
这部分代码展示了如何为树图添加比例尺，以表示遗传变化的量度。默认情况下，绘制的是进化距离或遗传变化显示在x轴上的系统发育图（phylogram）。如果您希望关闭分支长度的缩放并产生一个分类树（cladogram），可以在ggtree()调用中设置branch.length="none"选项。
自定义树的外观
# 绘制一个分类树，使用粗蓝色点线
ggtree(tree, branch.length="none", color="blue", size=2, linetype=3)
此代码段演示了如何使用ggtree自定义树图的外观，例如更改分支的颜色、大小和线型。在这个例子中，通过设置branch.length="none"生成了一个没有分支长度缩放的分类树，并使用蓝色、粗细为2、线型为点线的样式绘制分支。
为了完成这些练习，我们将使用ggtree包来创建不同样式的系统发育树。首先，确保你已经加载了ggtree和tidyverse包，并且已经导入了一个系统发育树对象tree。
练习 1：创建不同布局的系统发育树
创建倾斜的系统发育树
ggtree(tree, layout="slanted")
创建圆形的系统发育树
ggtree(tree, layout="circular")
创建圆形、未缩放的分类树，使用粗红线
ggtree(tree, layout="circular", branch.length="none") +
geom_treescale(color="red", size=2)
其他树的几何对象
添加额外的层。就像在ggplot2课程中一样，我们可以创建一个图表对象（例如p），来存储ggplot的基础布局，并根据需要添加更多层。下面我们将添加节点和尖端点，并最终标记尖端。
# 创建基础图表
p <- ggtree(tree)
# 添加节点点
p + geom_nodepoint()
# 添加尖端点
p + geom_tippoint()
# 标记尖端
p + geom_tiplab()
练习 2：自定义系统发育树的美学特征
按照要求创建具有以下美学特征的系统发育树：
尖端标签为紫色 紫色菱形尖端点 大型半透明黄色节点点 添加标题
# 创建基础图表并自定义
p <- ggtree(tree) +
geom_tiplab(color="purple") + # 尖端标签为紫色
geom_tippoint(shape=18, color="purple", size=3) + # 紫色菱形尖端点
geom_nodepoint(color="yellow", size=5, alpha=0.5) + # 大型半透明黄色节点点
ggtitle("Customized Phylogenetic Tree") # 添加标题
# 展示图表
p
在这里，shape=18代表菱形点，size参数控制点的大小，alpha参数控制透明度。这些练习展示了ggtree在系统发育树可视化和自定义中的灵活性和强大功能。树的注释是系统发育分析中一个重要的部分，它可以帮助我们理解树的结构和演化关系。ggtree提供了强大的工具来进行树的注释。geom_tiplab()函数添加了一些基本的注释，但是我们可以进一步进行注释。
内部节点编号
在深入了解ggtree如何内部处理树结构之前，了解内部节点编号是必要的。ggtree的某些函数，比如用于注释分支的函数，需要指定内部节点编号作为参数。要获取内部节点编号，可以使用geom_text函数，并通过aes(label=node)映射来显示，这里的"label"是树对象内部存储的"节点变量"（可以将其视为gapminder对象内的"continent"变量）。我们还提供hjust选项，以确保标签不会正好位于节点上。
ggtree(tree) + geom_text(aes(label=node), hjust=-.3)
另一种获取内部节点编号的方式是使用MRCA()函数，通过提供一个包含类群名的向量（使用c("taxon1", "taxon2")创建）。该函数将返回输入类群最近共同祖先（MRCA）的节点编号。首先，重新创建绘图，以便您可以选择您想要从中获取MRCA的类群。
ggtree(tree) + geom_tiplab()
让我们获取类群C+E和类群G+H的最近共同祖先。我们可以使用MRCA()来获取内部节点编号。之后，回到之前带节点标签的绘图以确认这一点。
MRCA(tree, tip=c("C", "E"))
# 输出可能是：[1] 17
MRCA(tree, tip=c("G", "H"))
# 输出可能是：[1] 21
通过这个过程，我们可以清楚地了解到系统发育树的结构，并且可以根据需要对特定的分支或节点进行注释。这对于突出显示树中的特定演化关系或分支群体特别有用。通过ggtree提供的工具，用户可以以直观和有效的方式探索和呈现他们的系统发育数据。
标记分支群（Clades）是在系统发育树注释中的一个重要步骤，它可以帮助我们高亮显示和标注树中特定的演化分支。ggtree提供了geom_cladelabel()函数，用于在选定的分支群上添加一个标签和相应的条形指示，以此来注释所选的分支群。你可以使用内部节点编号来选择分支群，这个编号是连接该分支群中所有分类单元的节点的编号。
标记特定分支群
下面的代码展示了如何为最近共同祖先在分类单元C和E之间的分支群（内部节点17）添加注释，并将标签设置为红色：
ggtree(tree) +
geom_cladelabel(node=17, label="Some random clade", color="red")
接下来，我们再次添加尖端标签。但你可能会注意到分支群标签和尖端标签之间距离太近了。我们可以通过添加offset参数来调整位置：
ggtree(tree) +
geom_tiplab() +
geom_cladelabel(node=17, label="Some random clade",
color="red2", offset=.8)
添加另一个分支群标签
现在，我们为连接分类单元G和H的分支群（内部节点21）添加另一个标签：
ggtree(tree) +
geom_tiplab() +
geom_cladelabel(node=17, label="Some random clade",
color="red2", offset=.8) +
geom_cladelabel(node=21, label="A different clade",
color="blue", offset=.8)
如果标签没有很好地对齐或者标签超出了图表的边缘，我们可以通过设置align=TRUE参数让标签对齐，并通过调整图表的x轴范围来确保所有标签都能完整显示：
ggtree(tree) +
geom_tiplab() +
geom_cladelabel(node=17, label="Some random clade",
color="red2", offset=.8, align=TRUE) +
geom_cladelabel(node=21, label="A different clade",
color="blue", offset=.8, align=TRUE) +
theme_tree2() +
xlim(0, 70) +
theme_tree()
使用geom_hilight()高亮整个分支群
另一个选择是使用geom_hilight()函数来高亮整个分支群，这为我们提供了一个更加直观的方式来强调树中的特定部分：
ggtree(tree) +
geom_tiplab() +
geom_hilight(node=17, fill="gold") +
geom_hilight(node=21, fill="purple")
通过使用ggtree包中的这些功能，我们可以灵活地对系统发育树进行详细的注释和定制，从而使得演化分析的结果更加清晰易懂。为了完成这个练习，我们首先需要确定两组分类单元（B+C和L+J）的最近共同祖先（MRCA）。接着，我们将使用ggtree包来绘制系统发育树，并通过各种geom图层来注释树，以展示演化事件和特定的分支群。
确定最近共同祖先（MRCA）
首先，使用MRCA()函数找出分类单元B+C和L+J的MRCA：
MRCA_bc <- MRCA(tree, tip=c("B", "C"))
MRCA_lj <- MRCA(tree, tip=c("L", "J"))
绘制和注释树
接下来，绘制树并添加注释：
# 绘制树并添加尖端标签
p <- ggtree(tree) +
geom_tiplab()
# 高亮特定分支群
p <- p + geom_hilight(node=MRCA_bc, fill="skyblue") +
geom_hilight(node=MRCA_lj, fill="lightgreen")
# 添加分支群标签
p <- p + geom_cladelabel(node=17, label="Superclade A-E",
color="red", offset=3, align=TRUE)
# 连接特定分类单元
p <- p + geom_taxalink("C", "E", color="gray", linetype=2) +
geom_taxalink("G", "J", color="gray", linetype=2, curvature=-.9)
# 添加比例尺条
p <- p + geom_treescale()
# 添加标题
p <- p + ggtitle("Evolutionary Events on Phylogenetic Tree")
# 改变布局（可选）
# p <- ggtree(tree, layout="circular") + ... # 在此处添加前面的所有注释代码
# 展示图表
p
这段代码展示了如何使用ggtree包的各种功能来创建一个带有详细注释的系统发育树。通过geom_hilight()我们可以高亮显示特定的分支群；geom_cladelabel()添加分支群标签；geom_taxalink()在特定的分类单元之间绘制连接线，展示可能的演化事件；最后，geom_treescale()添加了比例尺条。你可以根据需要调整offset和align参数以优化标签的位置，并使用ggtitle()添加标题。如果愿意，还可以通过设置layout="circular"在ggtree()函数中改变树的布局为圆形。
读取数据
在这个高级树注释的部分，我们将使用一组之前发表的数据集，来自Liang et al.的论文："Expansion of genotypic diversity and establishment of 2009 H1N1 pandemic-origin internal genes in pigs in China." 这个数据集包含了76个H3血凝素基因序列，这个谱系包含了猪和人类的流感A病毒。该数据集通过使用BEAST软件进行了重新分析，BEAST能够提供使用分子钟模型推断的有根、时间测量的系统发育树。
接下来，我们将导入BEAST输出的系统发育树，并为其添加注释。
# 读取BEAST生成的树文件
tree <- read.beast("data/flu_tree_beast.tree")
数据下载连接：https://4va.github.io/biodatasci/data.html这里使用read.beast()而不是之前用的read.tree()来读取树数据，因为BEAST软件输出的树文件包含了更多与时间相关的信息。
添加比例尺和调整布局
我们通过设置最新的采样日期mrsd="YYYY-MM-DD"来向ggtree()函数提供信息，这样可以在x轴上显示实际的日期。同时，我们添加比例尺条来表示遗传距离。
ggtree(tree, mrsd="2013-01-01") +
theme_tree2()
添加尖端标签和调整轴范围
为了添加尖端标签并使它们右对齐，同时减少标签连接线的粗细，我们可以使用geom_tiplab()并设置align=TRUE和linesize=.5。为了确保所有的标签都能够适当地显示在图中，我们还需要设置x轴的范围，以展示1990到2020年间的变化。
ggtree(tree, mrsd="2013-01-01") +
theme_tree2() +
geom_tiplab(align=TRUE, linesize=.5) +
xlim(1990, 2020)
msaplot是ggtree包中的一个函数，它能够将多序列比对（MSA）和系统发育树并排展示，为我们提供了一个直观的方式来查看序列的演化和基因的变异。这种视图尤其有用，当你想要探索特定基因或序列区域在不同物种中是如何变化的时候。
使用msaplot展示多序列比对和系统发育树
msaplot(p=ggtree(tree), fasta="data/flu_aasequence.fasta", window=c(150, 175))
在这个例子中，msaplot函数接受一个ggtree()生成的树对象p和一个指向FASTA格式的多序列比对文件的路径fasta。window参数用于指定你想在比对中查看的序列窗口。这里，我们只关注从第150个到第175个氨基酸残基的窗口。
尝试改变坐标系统
如果你想尝试一些有趣但可能不那么实用的视图，你可以通过在命令的末尾添加+ coord_polar(theta="y")来改变图的坐标系统，将它转换为极坐标系统：
msaplot(p=ggtree(tree), fasta="data/flu_aasequence.fasta", window=c(150, 175)) +
coord_polar(theta="y")
这将会将MSA和系统发育树的视图转换为圆形，为你提供一个不同的视角来观察序列的演化。虽然这种视图可能在实际应用中不是很常用，但它提供了一个有趣的方式来探索和展示你的数据。
更多高级注释
ggtree提供了许多高级功能来注释和定制系统发育树，包括添加额外的数据层、自定义节点和标签的样式，以及与其他生物信息学数据集进行集成。如果你对深入学习ggtree的高级注释功能感兴趣，请查看高级树注释使用说明（vignette）以获取更多信息和示例。这将为你提供更多工具和技巧，以便更有效地利用ggtree来探索和解释你的系统发育数据。
ggtree不仅允许您一次性绘制一个系统发育树，还可以同时绘制多个树，并且可以用ggplot2常规的方式来对它们进行分面展示。这为比较不同的树提供了一个方便的方法。此外，ggtree还提供了facet_plot()函数，使得可以将系统发育树与其他类型的数据面板并排展示，从而增加了更多数据的上下文理解。
绘制多个树
以下代码演示了如何生成不同尺寸的随机树，并使用facet_wrap()进行分面展示：
set.seed(42)
trees <- lapply(rep(c(10, 25, 50, 100), 3), rtree)
class(trees) <- "multiPhylo"
ggtree(trees) + facet_wrap(~.id, scale="free", ncol=4) + ggtitle("Many trees. Such phylogenetics. Wow.")
这段代码生成了10、25、50和100尖端的4棵随机树，每种尺寸重复3次，然后将它们一起绘制并分面展示。
与其他数据一起绘制树
接下来，我们将看看如何使用facet_plot()函数将系统发育树与其他面板上的自定义数据并排展示：
# 生成一个随机树
tree <- rtree(30)
# 创建基础图表
p <- ggtree(tree)
# 生成随机值数据
d1 <- data.frame(id=tree$tip.label, val=rnorm(30, sd=3))
# 添加点状图
p2 <- facet_plot(p, panel="dot", data=d1, geom=geom_point, aes(x=val), color='red3')
# 生成另一组随机值数据
d2 <- data.frame(id=tree$tip.label, value = abs(rnorm(30, mean=100, sd=50)))
# 添加柱状图
p3 <- facet_plot(p2, panel='bar', data=d2, geom=geom_segment,
aes(x=0, xend=value, y=y, yend=y), size=3, color='blue4')
# 展示全部图表并添加比例尺条
p3 + theme_tree2()
这段代码首先创建了一个30尖端的随机树，然后为每个尖端生成了一些随机值，并使用这些数据创建了点状图和柱状图。这些图表与系统发育树并排展示，提供了一个独特的视角来分析系统发育树上尖端的特性。
为了在系统发育树上叠加生物轮廓图像，ggtree提供了一种与phylopic.org资源集成的方法，后者提供了大量生物形象的免费轮廓图像。通过使用phylopic()函数，你可以将这些图像添加到你的系统发育树绘图中，从而提供更多的视觉信息。
以下是如何实现这一点的示例：
# 读取系统发育树数据
tree <- read.tree("data/tree_newick.nwk")
# 使用ggtree绘制树并添加生物轮廓图像
ggtree(tree) %>%
# 在整个图上添加一个革兰氏阴性细菌的轮廓图
phylopic("ba0a446e-18d7-4db9-9937-5adec24721b5", color="gold2", alpha = .25) %>%
# 在节点17添加Homo sapiens的轮廓图
phylopic("c089caae-43ef-4e4e-bf26-973dd4cb65c5", color="purple3", alpha = .5, node=17) %>%
# 在节点21添加狗的轮廓图
phylopic("6c9cb19d-1d8a-4215-88ba-d49cd4917a5e", color="purple3", alpha = .5, node=21)
在这个例子中，我们首先读取了Newick格式的系统发育树数据。然后，我们使用
ggtree()来绘制这棵树，并通过phylopic()函数来添加几个生物轮廓图。每个phylopic()函数调用都指定了一个轮廓图的唯一标识符（可以从phylopic.org获得），颜色和透明度，以及（对于后两个图像）要添加轮廓图的节点编号。
这种方法可以使得系统发育树的演示更加直观和生动，特别是在展示特定物种或演化关系的研究中。通过为树的不同部分叠加不同的生物轮廓图，我们可以清楚地展示系统发育树中的演化事件和分类关系。
回复tree_newick.nwk 获取树文件
https://mp.weixin.qq.com/s/zalmVPAWDjiYa-K1aOebCQ