library(ape)
library(phyloseq)
cat(
"(((Strix_aluco:4.2,Asio_otus:4.2):3.1,",
"Athene_noctua:7.3):6.3,Tyto_alba:13.5);",
file = "ex.tre",
sep = "\n"
)
tree.owls <- read.tree("ex.tre")
tree_layout() 与phylo 对象一起使用，而不是与绘制的系统发育一起使用。所以在plot.phylo(type='p') 与plot.phylo(type='u') 之后调用tree_layout() 将产生相同的对象。
plot.phylo(tree.owls, type = 'p')
a <- tree_layout(tree.owls)
plot.phylo(tree.owls, type = 'u')
b <- tree_layout(tree.owls)
identical(a, b)
[1] TRUE
要获得已经绘制的系统发育的坐标，我们可以像这样使用last_plot.phylo：
plot.phylo(tree.owls, type = 'p')
coords_phylogram <- get("last_plot.phylo",envir=.PlotPhyloEnv)
plot.phylo(tree.owls, type = 'u')
coords_unrooted <- get("last_plot.phylo",envir=.PlotPhyloEnv)
identical(coords_phylogram, coords_unrooted)
[1] FALSE
我们可以看到系统图和无根网络之间的坐标现在不同了。
get("last_plot.phylo",envir=.PlotPhyloEnv) 的输出是一个命名列表，因此我们可以轻松提取特定元素。例如，如果我们想在节点和提示的坐标处添加一些注解，我们可以这样做：
points(x = coords_unrooted$xx, y=coords_unrooted$yy, pch=21, bg="yellow", cex=3)
编辑：创建一个布局数据框（如tree_layout）并绘制一个无根系统发育
预赛。 phyloseq 不再需要，因为我们正在手动提取布局。
library(ape)
library(dplyr)
tr <- "(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3,Tyto_alba:13.5);"
tree.owls <- read.tree(text = tr)
创建tree_layouts() 的edgeDT 和vertDT 的函数。要绘制不仅水平或垂直的线段（就像您在分支图或系统图中的那样），我们需要另一列用于 y0 坐标，因为这些现在与 y 坐标不同。通过这种方式，我们可以绘制倾斜的系统图和无根树。
重要提示：此函数计算并存储在edgeDT 中的边长不正确（在某些情况下）！这是因为它们是通过减去坐标来计算的，如果边缘以一定角度绘制（如在倾斜的系统发育图或无根网络中），这显然不起作用。这对于绘图无关紧要，因为绘图功能不使用edge.lengths。但请记住。
layout_from_plot <- function(tree, type, drop_root=FALSE, use_vert=FALSE) {
plot.phylo(x = tree, type = type)
title("Plot to get coordinates")
coords <- get("last_plot.phylo",envir=.PlotPhyloEnv)
edgeDT <- tibble(
xright = coords$xx[coords$edge[,2]],
xleft = coords$xx[coords$edge[,1]],
y = coords$yy[coords$edge[,2]],
edge.length = xright-xleft,
V1 = coords$edge[,1],
V2 = coords$edge[,2]
)
edgeDT <- edgeDT %>%
arrange(V2) %>%
mutate(OTU = c(tree$tip.label, rep(NA_character_, coords$Nnode - 1))) %>%
select(V1, V2, edge.length, OTU, xleft, xright, y)
if (!use_vert) {
edgeDT <- mutate(edgeDT, y0=y[V1-1])
} else {
edgeDT <- mutate(edgeDT, y0=y)
}
# this is a bit hacky, but it does work.
# double check that OTUs are in the right positions
# and branch lengths are correct.
# ideally, you would unroot the tree, then plot the network
# to extract coordinates. See below.
if (is.rooted(phy = tree)) {
if (drop_root) {
edgeDT <- edgeDT %>%
group_by(y0) %>%
mutate(xleft = ifelse(y == y0, xright, xleft)) %>%
mutate(xright = ifelse(y == y0, lead(xright), xright)) %>%
mutate(y = ifelse(y == y0, lead(y), y)) %>%
distinct(y, y0, .keep_all = TRUE)
}
}
vertDT <- edgeDT %>%
group_by(V1) %>%
mutate(vmin=min(y), vmax=max(y)) %>%
mutate(x=xleft[which(y==min(y))]) %>%
select(V1, x, vmin, vmax) %>%
distinct()
return(list("edgeDT"=edgeDT, "vertDT"=vertDT))
}
根据提取的树布局绘制不同类型树的函数。
plot_from_layout <- function(tree_ly, plot_vert=FALSE) {
plot(
1,
type = 'n',
axes = TRUE,
xlim = c(0, max(tree_ly$edgeDT$xright)),
ylim = c(0, max(c(tree_ly$edgeDT$y, tree_ly$edgeDT$y0))),
main = 'tree_layout'
)
segments(
x0 = tree_ly$edgeDT$xleft,
y0 = tree_ly$edgeDT$y0,
x1 = tree_ly$edgeDT$xright,
y1 = tree_ly$edgeDT$y
)
if (plot_vert) {
segments(
x0 = tree_ly$vertDT$x,
y0 = tree_ly$vertDT$vmin,
x1 = tree_ly$vertDT$x,
y1 = tree_ly$vertDT$vmax
)
}
}
测试。
par(mfrow=c(1,2))
layout_from_plot(
tree = tree.owls,
type = "p",
drop_root = FALSE,
use_vert = TRUE
) %>%
plot_from_layout(tree_ly = ., plot_vert = TRUE)
layout_from_plot(
tree = tree.owls,
type = "c",
drop_root = FALSE,
use_vert = FALSE
) %>%
plot_from_layout(tree_ly = ., plot_vert = FALSE)
layout_from_plot(
tree = tree.owls,
type = "u",
drop_root = TRUE,
use_vert = FALSE
) %>%
plot_from_layout(tree_ly = ., plot_vert = FALSE)
最好的方法不是手动移除边缘（如上），而是先将树解根，然后绘制它。为此，您需要修改提取坐标的方式。主要问题是根节点如何分割边缘以及移除根节点后将剩余的添加到哪个分支。默认情况下，这个额外的边缘似乎被添加到网络的一个内部边缘，但考虑到那里的系统图和分支长度，这似乎并不正确。见下文。
owls.unrooted <- unroot(tree.owls)
layout_from_plot(
tree = owls.unrooted,
type = "u",
drop_root = FALSE,
use_vert = FALSE
) %>%
plot_from_layout(tree_ly = ., plot_vert = FALSE)
编辑 2：更新 layout_from_plot 函数以正确处理 y0 坐标并使用无根树而不是手动删除边缘。
我意识到为y0 坐标创建列过于复杂。所需要做的就是从get("last_plot.phylo",envir=.PlotPhyloEnv) 的输出中获取edge 第一列中的yy 坐标。现在，rooted 和 unrooted 树都可以正常工作了。
更新功能：
layout_from_plot <- function(tree, type, drop_root=FALSE, use_vert=FALSE) {
if (drop_root) {
tree <- unroot(tree)
}
plot.phylo(x = tree, type = type)
title("Plot to get coordinates")
coords <- get("last_plot.phylo",envir=.PlotPhyloEnv)
edgeDT <- tibble(
xright = coords$xx[coords$edge[,2]],
xleft = coords$xx[coords$edge[,1]],
y = coords$yy[coords$edge[,2]],
y0 = coords$yy[coords$edge[,1]],
edge.length = xright-xleft,
V1 = coords$edge[,1],
V2 = coords$edge[,2]
)
edgeDT <- edgeDT %>%
arrange(V2) %>%
mutate(OTU = c(tree$tip.label, rep(NA_character_, coords$Nnode - 1))) %>%
select(V1, V2, edge.length, OTU, xleft, xright, y, y0)
if (use_vert) {
edgeDT <- mutate(edgeDT, y0=y)
}
vertDT <- edgeDT %>%
group_by(V1) %>%
mutate(vmin=min(y), vmax=max(y)) %>%
mutate(x=xleft[which(y==min(y))]) %>%
select(V1, x, vmin, vmax) %>%
distinct()
return(list("edgeDT"=edgeDT, "vertDT"=vertDT))
}
附加测试代码：
layout_from_plot(tree = bird.orders, type = "c", use_vert = FALSE, drop_root = TRUE) %>%
plot_from_layout(tree_ly = ., plot_vert = FALSE)
layout_from_plot(tree = bird.orders, type = "u", drop_root = TRUE, use_vert = FALSE) %>%
plot_from_layout(tree_ly = ., plot_vert = FALSE)