禁止商业或二改转载，仅供自学使用，侵权必究，如需截取部分内容请后台联系作者!
教程
教程内容
介绍
构成要素
有根树与无根树
构建方法
应用领域
说明的问题
加载R包
数据下载
导入数据
数据预处理
系统发育树可视化
准备画图数据
1.构建基础系统发育树 p1
2.添加条形图 p2
3.添加热图 p3
4.添加第二个热图 p4
保存PDF
总结
系统信息
介绍
物种系统发育树（Phylogenetic tree），也称为进化树或系统进化树，是一种以树状分支图形来表示各物种或基因之间的亲缘关系的图表。它利用生物的形态特征、分子序列（如DNA、RNA或蛋白质序列）等数据，通过数理统计算法来计算生物之间的进化关系，从而构建出一个反映物种进化历史的拓扑结构。
构成要素
系统发育树由节点（node）和进化分支（branch）组成：
节点：表示一个分类学单元，如属、种群、个体等。分支末端的节点对应一个基因或者生物体，代表实际观察到的最终分类。
进化分支：定义了分类单元（祖先与后代）之间的关系，一个分支只能连接两个相邻的节点。
分支长度：表示该分支在进化过程中的变化程度。标有分支长度的进化分支称为标度枝（scaled branch）。校正后的标度树（scaled tree）常常用年代表示，这样的树通常根据某一或部分基因的理论分析而得出。
有根树与无根树
系统发育树可以是有根的（rooted），也可以是无根的（unrooted）：
有根树：有一个明确的根节点，表示所有物种的共同祖先。这种树可以清晰地显示物种的进化方向。
无根树：没有明确的根节点，仅表示物种之间的亲缘关系，不显示进化方向。
构建方法
构建系统发育树的方法主要有以下几种：
距离法（Distance-based methods）：基于物种之间的进化距离，通过计算和比较序列之间的差异来构建树。常用的方法包括邻接法（Neighbor-Joining, NJ）和最小进化法（Minimum Evolution, ME）。
最大简约法（Maximum Parsimony, MP）：寻找需要最少进化步骤的树，即假设进化过程中变化的次数最少。
最大似然法（Maximum Likelihood, ML）：基于统计模型，寻找最有可能产生观测数据的树。
贝叶斯推断法（Bayesian Inference, BI）：利用贝叶斯统计方法，通过计算后验概率来构建树。
应用领域
系统发育树在多个领域有广泛应用：
生物分类学：帮助确定物种的分类地位，构建生物的分类系统。
进化生物学：研究物种的起源、分化和进化历史，了解生物的进化过程。
生态学：研究物种的地理分布和生态适应性，了解物种在不同环境中的进化和适应机制。
分子生物学：研究基因家族的进化，了解基因的功能和进化历史。
医学：研究病原体的进化和传播，指导药物研发和疾病防控。
说明的问题
系统发育树主要说明以下问题：
亲缘关系：显示不同物种或基因之间的亲缘关系，帮助理解它们的共同祖先和进化路径。
进化历史：反映物种的进化历程，包括分化时间、进化速率和关键进化事件。
分类地位：确定物种在生物分类系统中的位置，为生物分类提供科学依据。
功能预测：通过比较基因或蛋白质的进化关系，预测其功能和结构特征。
生态适应：研究物种在不同生态环境中的适应性进化，了解生态系统的演变过程。
加载R包
安装ggtree的时候注意R包的依赖问题，耐心等待安装。
library(ape)
library(ggnewscale)
library(ggplot2)
library(ggtree)
library(ggtreeExtra)
library(RColorBrewer)
library(tidyverse)
# devtools::install_github("kevinwolz/hisafer")
library(hisafer)
数据下载
R语言系统发育树数据集下载链接：
百度网盘链接: 百度网盘下载链接
提取码：R语言绘制微生物物种系统发育树（phylogenetic tree）
导入数据
本次加载数据有四种数据类型，它们分别是
树结构文件：pMAGs_bact_gtdtk_midroot.tree；
物种信息：pMAGS_tax.tsv；
物种出现情况：pMAGS_Presence_Datasets.txt；
物种加和：pMAGs_cov_sum.txt
BacTree <- read.tree("pMAGs_bact_gtdtk_midroot.tree")
dat <- read_tsv("pMAGS_tax.tsv")
dataset <- read_tsv("pMAGS_Presence_Datasets.txt")
covmax <- read_tsv("pMAGs_cov_sum.txt")
数据预处理
修改物种phylum水平的名称；
选择相对丰度前16的物种；
过滤物种类型
dat$p_c <- if_else(dat$p == "p__Proteobacteria", dat$c, dat$p)
dat$p_c <- gsub(".__", "", dat$p_c, perl = T)
dat$p_c <- gsub("_.$", "", dat$p_c, perl = T)
bacDat <- dat %>%
filter(d == "d__Bacteria") %>%
mutate(Abundance = 0.2)
mostAbundantTax <- bacDat %>%
group_by(p_c) %>%
summarise(total = n()) %>%
arrange(desc(total)) %>%
slice(1:16)
list <- mostAbundantTax$p_c
list <- c(list, "Nitrospirota")
bacDat <- bacDat %>%
mutate(p_c = if_else(p_c %in% list, p_c, "Other"))
bacDatset <- dataset %>%
filter(MAGs %in% bacDat$MAGs)
bacDatset$Busi <-
gsub("Busi", "Busi et al., Nat Com, 2022", bacDatset$Busi)
bacDatset$ENSEMBLE <-
gsub("ENSEMBLE", "Michoud et al, L&O, 2023", bacDatset$ENSEMBLE)
bacDatset$Tibet <-
gsub("Tibet", "Tibetan Glacier Genome and Gene", bacDatset$Tibet)
bacDatset$Tara <- gsub("Tara", "Tara Oceans", bacDatset$Tara)
bacDatset <- as.data.frame(bacDatset)
rownames(bacDatset) <- bacDatset$MAGs
bacDatset$MAGs <- NULL
bactcov <- covmax %>%
filter(MAGs %in% bacDat$MAGs) %>%
select(MAGs, count)
bactcov <- as.data.frame(bactcov)
rownames(bactcov) <- bactcov$MAGs
bactcov$MAGs <- NULL
系统发育树可视化
构建和美化了一个微生物物种的系统发育树，并添加了多个图层来展示不同的数据。
p1：基础的圆形系统发育树，节点颜色根据 group 变量着色。
p2：在 p1 的基础上添加了条形图，展示每个节点的 MAGs 值。
p3：在 p2 的基础上添加了热图，展示 bacDatset 数据。
p4：在 p3 的基础上添加了第二个热图，展示 bactcov 数据。
准备画图数据
标准化物种的频次；
提取物种的树结构；
设置物种的颜色；
bactcov$count <- log10(bactcov$count)
a <- split(bacDat$MAGs, bacDat$p_c)
tree <- groupOTU(BacTree, a)
getPaletteBact <- colorRampPalette(brewer.pal(9, "Set1"))
bactColor <- getPaletteBact(length(unique(bacDat$p_c)) + 1)
bactColor[1] <- "black"
1. 构建基础系统发育树 p1
ggtree(tree, layout = "circular", aes(color = group))：使用 ggtree 函数创建一个圆形布局的系统发育树，节点颜色根据 group 变量进行着色。
geom_tree()：添加树的线条。
theme_tree()：应用 ggtree 的默认主题。
geom_treescale(width = 0.1)：添加一个比例尺，宽度为 0.1。
scale_color_manual(values = bactColor, na.value = "transparent", guide = "none")：自定义节点颜色，缺失值透明，不显示图例。
theme(legend.position = "right")：将图例位置设置在右侧。
p1 <-
ggtree(tree,
layout = "circular",
aes(color = group)
) + #
geom_tree() +
theme_tree() +
geom_treescale(width = 0.1) +
scale_color_manual(
values = bactColor,
na.value = "transparent",
guide = "none"
) +
# geom_text2(aes(subset=!isTip, label=node), hjust=-.3)+
theme(legend.position = "right")
p1
2. 添加条形图 p2
new_scale_colour() 和 new_scale_fill()：重置颜色和填充尺度，以便添加新的图层。
geom_fruit：在树的每个节点上添加条形图，数据来自 bacDat，条形图的宽度为 0.01，高度为 MAGs，填充颜色根据 p_c 变量。
scale_fill_manual(values = bactColor[-1])：自定义条形图的填充颜色。
labs(fill = "Taxa")：设置填充颜色的图例标签为 "Taxa"。
p2 <- p1 +
new_scale_colour() +
new_scale_fill() +
geom_fruit(
data = bacDat,
pwidth = 0.01,
geom = geom_bar,
mapping = aes(y = MAGs, fill = p_c, x = 1),
# orientation="y",
stat = "identity",
) +
scale_fill_manual(values = bactColor[-1]) +
labs(fill = "Taxa") +
new_scale_colour() +
new_scale_fill()
p2
3. 添加热图 p3
gheatmap：在 p2 的基础上添加一个热图，数据来自 bacDatset，热图宽度为 0.2，偏移量为 0.1，不显示列名，颜色默认。
scale_fill_discrete(na.translate = FALSE)：处理缺失值，不翻译为颜色。
labs(fill = "Datasets")：设置填充颜色的图例标签为 "Datasets"。
p3 <-
gheatmap(
p2,
bacDatset,
width = 0.2,
offset = 0.1,
# offset=8, width=0.6,
colnames = FALSE,
color = NULL
) +
scale_fill_discrete(na.translate = F) +
labs(fill = "Datasets") +
new_scale_colour() +
new_scale_fill()
p3
4. 添加第二个热图 p4
gheatmap：在 p3 的基础上添加另一个热图，数据来自 bactcov，热图宽度为 0.05，偏移量为 0.6，不显示列名，颜色默认。
scale_fill_viridis_c()：使用 viridis 色板填充热图。
labs(fill = "Normalized log10\nabundance")：设置填充颜色的图例标签为 "Normalized log10 abundance"。
p4 <- gheatmap(
p3,
bactcov,
offset = 0.6,
width = 0.05,
colnames = FALSE,
color = NULL
) +
scale_fill_viridis_c() +
labs(fill = "Normalized log10\nabundance")
p4
结果：微生物物种系统发育树图，它由以下部分组成：
最外层：物种的相对丰度值；
次外层：物种来自于的数据集分布；
次里层：物种的phylum水平；
最里层：不同分类物种的系统发育树结构。
保存PDF
将图片以PDF格式保存到本地
ggsave_fitmax <- function(
filename,
plot,
maxheight = 7,
maxwidth = maxheight,
units = "in", ...) {
if (is.null(plot)) {
return(FALSE)
}
dims <- get_dims(
ggobj = plot,
maxheight = maxheight,
maxwidth = maxwidth,
units = units
)
ggplot2::ggsave(
filename = filename,
plot = plot,
height = dims$height,
width = dims$width,
units = units, ...
)
}
get_dims <- function(
ggobj,
maxheight,
maxwidth = maxheight,
units = "in", ...) {
# Internal helper function:
# Treat all null units in a unit object as if they were inches.
# This is a bad idea in gneral, but I use it here as a workaround.
# Extracting unit names from non-atomic unit objects is a pain,
# so questions like "which rows of this table layout have null heights?"
# are hard to answer. To work around it, I exploit an (undocumented!)
# quirk: When calculating the size of a table layout inside a Grid plot,
# convertUnit(...) treats null units as zero.
# Therefore
# (convertHeight(grob_height, "in", valueOnly=TRUE)
# - convertHeight(null_as_if_inch(grob_height), "in", valueOnly=TRUE))
# does the inverse of convertUnit: It gives the sum of all *null* heights
# in the object, treating *fixed* units as zero.
#
# Warning: I repeat, this approach ONLY makes any sense if
# convertUnit(unit(1, "null"), "in", "x", valueOnly=T) == 0
# is true. Please check that it is before calling this code.
.null_as_if_inch <- function(u) {
stopifnot(packageVersion("grid") < "4.0")
if (!grid::is.unit(u)) {
return(u)
}
if (is.atomic(u)) {
if ("null" %in% attr(u, "unit")) {
d <- attr(u, "data")
u <- unit(
x = as.vector(u),
units = gsub("null", "in", attr(u, "unit")),
data = d
)
}
return(u)
}
if (inherits(u, "unit.arithmetic")) {
l <- .null_as_if_inch(u$arg1)
r <- .null_as_if_inch(u$arg2)
if (is.null(r)) {
args <- list(l)
} else {
args <- list(l, r)
}
return(do.call(u$fname, args))
}
if (inherits(u, "unit.list")) {
return(do.call(grid::unit.c, lapply(u, .null_as_if_inch)))
}
return(u)
}
if (inherits(ggobj, "ggplot") && !isTRUE(ggobj$respect) &&
is.null(ggobj$theme$aspect.ratio) && is.null(ggobj$coordinates$ratio) &&
is.null(ggplot2::theme_get()$aspect.ratio)) {
return(list(height = maxheight, width = maxwidth))
}
tmpf <- tempfile(pattern = "dispos-a-plot", fileext = ".png")
png(
filename = tmpf,
height = maxheight,
width = maxwidth,
units = units,
res = 120, ...
)
on.exit({
dev.off()
unlink(tmpf)
})
if (inherits(ggobj, "ggplot")) {
g <- ggplot2::ggplotGrob(ggobj)
} else if (inherits(ggobj, "gtable")) {
g <- ggobj
} else {
stop("Don't know how to get sizes for object of class ", deparse(class(ggobj)))
}
stopifnot(grid::convertUnit(grid::unit(1, "null"), "in", "x", valueOnly = TRUE) == 0)
known_ht <- sum(grid::convertHeight(g$heights, units, valueOnly = TRUE))
known_wd <- sum(grid::convertWidth(g$widths, units, valueOnly = TRUE))
free_ht <- maxheight - known_ht
free_wd <- maxwidth - known_wd
if (packageVersion("grid") >= "4.0.0") {
null_rowhts <- as.numeric(g$heights[grid::unitType(g$heights) == "null"])
null_colwds <- as.numeric(g$widths[grid::unitType(g$widths) == "null"])
panel_asps <- (
matrix(null_rowhts, ncol = 1)
%*% matrix(1 / null_colwds, nrow = 1))
} else {
all_null_rowhts <- (
grid::convertHeight(.null_as_if_inch(g$heights), "in", valueOnly = TRUE)
- grid::convertHeight(g$heights, "in", valueOnly = TRUE))
all_null_colwds <- (
grid::convertWidth(.null_as_if_inch(g$widths), "in", valueOnly = TRUE)
- grid::convertWidth(g$widths, "in", valueOnly = TRUE))
null_rowhts <- all_null_rowhts[all_null_rowhts > 0]
null_colwds <- all_null_colwds[all_null_colwds > 0]
panel_asps <- (matrix(null_rowhts, ncol = 1) %*% matrix(1 / null_colwds, nrow = 1))
}
panel_asps <- matrix(null_rowhts, ncol = 1) %*% matrix(1 / null_colwds, nrow = 1)
max_rowhts <- free_ht / sum(null_rowhts) * null_rowhts
max_colwds <- free_wd / sum(null_colwds) * null_colwds
rowhts_if_maxwd <- max_colwds[1] * panel_asps[, 1]
colwds_if_maxht <- max_rowhts[1] / panel_asps[1, ]
height <- min(maxheight, known_ht + sum(rowhts_if_maxwd))
width <- min(maxwidth, known_wd + sum(colwds_if_maxht))
return(list(height = height, width = width))
}
ggsave_fitmax("bacterialTree.pdf",
p4,
maxheight = 15
)
总结
本教程将指导你如何使用ggtree等一系列包在R语言环境中构建微生物物种的系统发育树。为了帮助读者更好地理解和应用，本教程提供了完整的数据和代码示例。
系统信息
R version 4.3.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.2
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Asia/Shanghai
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] hisafer_1.5.1 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[6] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 tidyverse_2.0.0
[11] RColorBrewer_1.1-3 ggtreeExtra_1.12.0 ggtree_3.8.2 ggplot2_3.5.1 ggnewscale_0.4.10
[16] ape_5.7-1
loaded via a namespace (and not attached):
[1] yulab.utils_0.1.4 utf8_1.2.4 generics_0.1.3 ggplotify_0.1.2 stringi_1.8.3
[6] lattice_0.22-6 hms_1.1.3 digest_0.6.35 magrittr_2.0.3 timechange_0.3.0
[11] grid_4.3.3 fastmap_1.1.1 jsonlite_1.8.8 fansi_1.0.6 aplot_0.2.2
[16] scales_1.3.0 lazyeval_0.2.2 cli_3.6.3 rlang_1.1.4 munsell_0.5.0
[21] tidytree_0.4.6 withr_3.0.2 cachem_1.0.8 tools_4.3.3 parallel_4.3.3
[26] tzdb_0.4.0 memoise_2.0.1 colorspace_2.1-0 vctrs_0.6.5 R6_2.5.1
[31] gridGraphics_0.5-1 lifecycle_1.0.4 fs_1.6.3 ggfun_0.1.4 treeio_1.26.0
[36] pkgconfig_2.0.3 pillar_1.9.0 gtable_0.3.4 glue_1.8.0 Rcpp_1.0.12
[41] tidyselect_1.2.1 rstudioapi_0.16.0 farver_2.1.1 nlme_3.1-164 patchwork_1.3.0
[46] compiler_4.3.3