1

有没有办法用phyloseq::tree_layout()函数获取无根树的结构？

使用tree_layout()将为您提供构成下面绘制的树的节点和段的坐标。然后，您可以轻松地重绘该树。

library(ape)
library(phyloseq) # bioconductor
#
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")


par(mfrow=c(2,1))

#original tree
plot.phylo(tree.owls,type = 'p',main = 'plot.phylo')

#redraw structure with tree_layout object
tree.ly <- tree_layout(tree.owls)
plot(1,type='n',axes=FALSE, xlim = c(0,max(tree.ly$edgeDT$xright)),ylim = c(0,max(tree.ly$edgeDT$y)),main = 'tree_layout')
segments(x0=tree.ly$edgeDT$xleft,y0=tree.ly$edgeDT$y,x1=tree.ly$edgeDT$xright,y1=tree.ly$edgeDT$y)
segments(x0=tree.ly$vertDT$x,y0=tree.ly$vertDT$vmin,x1=tree.ly$vertDT$x,y1=tree.ly$vertDT$vmax)

但是如果你想重绘这棵树怎么办：plot.phylo(tree.owls,type = 'u'). 你会怎么做？

4

1 回答 1

1 
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')vs.之后调用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由 的第一列排序的坐标。现在，有根树和无根树都可以正常工作。edgeget("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)

在此处输入图像描述

于 2019-09-20T21:40:41.067 回答