R/bind.tree.R

   1 ## bind.tree.R (2012-02-13)
   2
   3 ##    Bind Trees
   4
   5 ## Copyright 2003-2012 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 `+.phylo` <- function(x, y)
  11 {
  12     p <- if (is.null(x$root.edge)) 0 else x$root.edge
  13     bind.tree(x, y, position = p)
  14 }
  15
  16 bind.tree <- function(x, y, where = "root", position = 0, interactive = FALSE)
  17 {
  18     nx <- length(x$tip.label)
  19     mx <- x$Nnode
  20     ROOTx <- nx + 1L
  21     ny <- length(y$tip.label)
  22     my <- y$Nnode
  23
  24     if (interactive) {
  25         lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
  26         if (lastPP$type != "phylogram" || lastPP$direction != "rightwards")
  27             stop("you must plot tree 'x' as a 'rightward phylogram'")
  28         cat("Click where you want to graft tree 'y'...\n")
  29         xy <- locator(1)
  30         d <- abs(xy$y - lastPP$yy)
  31         d[lastPP$xx - xy$x < 0] <- Inf
  32         where <- which.min(d)
  33         position <- lastPP$xx[where] - xy$x
  34         if (position < 0) position <- 0
  35         cat("The following parameters are used:\n")
  36         cat("  where =", where, " position =", position, "\n")
  37     } else {
  38         if (where == 0 || where == "root") where <- ROOTx
  39         if (position < 0) position <- 0
  40         if (where > nx + mx)
  41             stop("argument 'where' out of range for tree 'x'")
  42     }
  43
  44     ## check whether both trees have branch lengths:
  45     switch(is.null(x$edge.length) + is.null(y$edge.length) + 1L,
  46            wbl <- TRUE, {
  47                x$edge.length <- y$edge.length <- NULL
  48                wbl <- FALSE
  49                warning("one tree has no branch lengths, they have been ignored")
  50            },
  51            wbl <- FALSE)
  52
  53     yHasNoRootEdge <- is.null(y$root.edge)
  54     xHasNoRootEdge <- is.null(x$root.edge)
  55
  56     ## find the row of 'where' before renumbering
  57     if (where == ROOTx) case <- 1 else {
  58         i <- which(x$edge[, 2] == where)
  59         case <- if (where <= nx) 2 else 3
  60     }
  61     ## case = 1 -> y is bound on the root of x
  62     ## case = 2 -> y is bound on a tip of x
  63     ## case = 3 -> y is bound on a node of x
  64
  65     ## check that 'position' is correct
  66     if (position && wbl) {
  67 ### New in ape 3.0-1: this makes possible binding 'y' below
  68 ### a node of 'x' thus creating a new node in 'x'
  69 ###        if (!wbl)
  70 ###            stop("'position' is non-null but trees have no branch lengths")
  71         if (case == 1) {
  72             if (xHasNoRootEdge)
  73                 stop("tree 'x' has no root edge")
  74             if (position > x$root.edge)
  75                 stop("'position' is larger than x's root edge")
  76         } else {
  77             if (x$edge.length[i] < position)
  78                 stop("'position' is larger than the branch length")
  79         }
  80     }
  81
  82     ## the special case of substituting two tips:
  83     if (case == 2 && ny == 1 && !position) {
  84         x$tip.label[x$edge[i, 2]] <- y$tip.label
  85         if (wbl)
  86             x$edge.length[i] <- x$edge.length[i] + y$edge.length
  87         return(x)
  88     }
  89
  90     x <- reorder(x)
  91     y <- reorder(y)
  92
  93 ### because in all situations internal nodes need to be
  94 ### renumbered, they are changed to negatives first, and
  95 ### nodes eventually added will be numbered sequentially
  96
  97     nodes <- x$edge > nx
  98     x$edge[nodes] <- -(x$edge[nodes] - nx) # -1, ..., -mx
  99     nodes <- y$edge > ny
 100     y$edge[nodes] <- -(y$edge[nodes] - ny + mx) # -(mx+1), ..., -(mx+my)
 101     ROOT <- -1L # may change later
 102     next.node <- -(mx + my) - 1L
 103
 104     ## renumber now the tips in y:
 105     new.nx <- if (where <= nx && !position) nx - 1L else nx
 106     y$edge[!nodes] <- y$edge[!nodes] + new.nx
 107
 108     ## if 'y' as a root edge, use it:
 109     if (!yHasNoRootEdge) {
 110         y$edge <- rbind(c(0, y$edge[1]), y$edge)
 111         ##                ^ will be filled later
 112         next.node <- next.node - 1L
 113         if (wbl) y$edge.length <- c(y$root.edge, y$edge.length)
 114     }
 115
 116     switch(case, { # case = 1
 117         if (position) {
 118             x$root.edge <- x$root.edge - position
 119             x$edge <- rbind(c(next.node, x$edge[1]), x$edge)
 120             ROOT <- next.node
 121             if (wbl) x$edge.length <- c(position, x$edge.length)
 122         }
 123         if (yHasNoRootEdge) {
 124             j <- which(y$edge[, 1] == y$edge[1])
 125             y$edge[j, 1] <- ROOT
 126         } else y$edge[1] <- ROOT
 127         x$edge <- rbind(x$edge, y$edge)
 128         if (wbl)
 129             x$edge.length <- c(x$edge.length, y$edge.length)
 130     }, { # case = 2
 131         if (position) {
 132             x$edge[i, 2] <- next.node
 133             x$edge <- rbind(x$edge[1:i, ], c(next.node, where), x$edge[-(1:i), ])
 134             if (wbl) {
 135                 x$edge.length[i] <- x$edge.length[i] - position
 136                 x$edge.length <- c(x$edge.length[1:i], position, x$edge.length[-(1:i)])
 137             }
 138             i <- i + 1L
 139             if (yHasNoRootEdge) {
 140                 j <- which(y$edge[, 1] == y$edge[1])
 141                 y$edge[j, 1] <- x$edge[i, 1]
 142             } else y$edge[1] <- x$edge[i, 1]
 143         } else {
 144             if (yHasNoRootEdge) x$edge[i, 2] <- y$edge[1]
 145             else {
 146                 ## the root edge of y is fused with the terminal edge of x
 147                 if (wbl) y$edge.length[1] <- y$edge.length[1] + x$edge.length[i]
 148                 y$edge[1] <- x$edge[i, 1]
 149                 ## delete i-th edge in x:
 150                 x$edge <- x$edge[-i, ]
 151                 if (wbl) x$edge.length <- x$edge.length[-i]
 152                 i <- i - 1L
 153             }
 154             x$tip.label <- x$tip.label[-where]
 155             ## renumber the tips that need to:
 156             ii <- which(x$edge[, 2] > where & x$edge[, 2] <= nx)
 157             x$edge[ii, 2] <- x$edge[ii, 2] - 1L
 158         }
 159         x$edge <- rbind(x$edge[1:i, ], y$edge, x$edge[-(1:i), ])
 160         if (wbl)
 161             x$edge.length <- c(x$edge.length[1:i], y$edge.length, x$edge.length[-(1:i)])
 162     }, { # case = 3
 163         if (position) {
 164             if (yHasNoRootEdge) {
 165                 j <- which(y$edge[, 1] == y$edge[1])
 166                 y$edge[j, 1] <- next.node
 167             } else y$edge[1] <- next.node
 168             x$edge <- rbind(x$edge[1:i, ], c(next.node, x$edge[i, 2]), x$edge[-(1:i), ])
 169             x$edge[i, 2] <- next.node
 170             if (wbl) {
 171                 x$edge.length[i] <- x$edge.length[i] - position
 172                 x$edge.length <- c(x$edge.length[1:i], position, x$edge.length[-(1:i)])
 173             }
 174             i <- i + 1L
 175         } else {
 176             if (yHasNoRootEdge) {
 177                 j <- which(y$edge[, 1] == y$edge[1])
 178                 y$edge[j, 1] <- x$edge[i, 2]
 179             } else y$edge[1] <- x$edge[i, 2]
 180         }
 181         x$edge <- rbind(x$edge[1:i, ], y$edge, x$edge[-(1:i), ])
 182         if (wbl)
 183             x$edge.length <- c(x$edge.length[1:i], y$edge.length, x$edge.length[-(1:i)])
 184     })
 185
 186     x$tip.label <- c(x$tip.label, y$tip.label)
 187
 188     if (is.null(x$node.label)) {
 189         if (!is.null(y$node.label))
 190             x$node.label <- c(rep(NA, mx), y$node.label)
 191     } else {
 192         x$node.label <-
 193             if (is.null(y$node.label)) c(x$node.label, rep(NA, my))
 194             else c(x$node.label, y$node.label)
 195     }
 196
 197     n <- length(x$tip.label)
 198     x$Nnode <- dim(x$edge)[1] + 1L - n
 199
 200     ## update the node labels before renumbering (this adds NA for
 201     ## the added nodes, and drops the label for those deleted)
 202     if (!is.null(x$node.label))
 203         x$node.label <- x$node.label[sort(-unique(x$edge[, 1]))]
 204
 205     ## renumber nodes:
 206     newNb <- integer(x$Nnode)
 207     newNb[-ROOT] <- n + 1L
 208     sndcol <- x$edge[, 2] < 0
 209     ## executed from right to left, so newNb is modified before x$edge:
 210     x$edge[sndcol, 2] <- newNb[-x$edge[sndcol, 2]] <- n + 2:x$Nnode
 211     x$edge[, 1] <- newNb[-x$edge[, 1]]
 212
 213     if (!is.null(x$node.label))
 214         x$node.label <- x$node.label[order(newNb[newNb > 0])]
 215
 216     x
 217 }