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, а не с сюжетной филогенией.Поэтому вызов tree_layout() после plot.phylo(type='p') против plot.phylo(type='u') приведет к тому же объекту.
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.Все, что нужно, это взять yy координаты, упорядоченные из первого столбца edge, из вывода get("last_plot.phylo",envir=.PlotPhyloEnv).Теперь все работает нормально с корневыми и нерутированными деревьями.
Обновленная функция:
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)
