R/AddTip.R

In TreeTools: Create, Modify and Analyse Phylogenetic Trees

#' Add a tip to a phylogenetic tree
#'
#' `AddTip()` adds a tip to a phylogenetic tree at a specified location.
#'
#' `AddTip()` extends \code{\link{bind.tree}}, which cannot handle
#'   single-taxon trees.
#'
#' @template treeParam
#' @param where The node or tip that should form the sister taxon to the new
#' node.  To add a new tip at the root, use `where = 0`.  By default, the
#' new tip is added to a random edge.
#' @param label Character string providing the label to apply to the new tip.
#' @param nodeLabel Character string providing a label to apply to the newly
#' created node, if `tree$node.label` is specified.
#' @param edgeLength Numeric specifying length of new edge. If `NULL`,
#' defaults to `lengthBelow`.
# Notice added in v1.10.0.9001, 2024-02-20:
#' This will become the default behaviour in a
#' future release; please manually specify the desired behaviour in your code.
#' @param lengthBelow Numeric specifying length below neighbour at which to
#' graft new edge. Values greater than the length of the edge will result
#' in negative edge lengths. If `NULL`, the default, the new tip will be added
#' at the midpoint of the broken edge. If inserting at the root (`where = 0`),
#' a new edge of length `lengthBelow` will be inserted.
#' @param nTip,nNode,rootNode Optional integer vectors specifying number of tips
#' and nodes in `tree`, and index of root node.
#' Not checked for correctness: specifying values here yields a marginal speed
#' increase at the cost of code safety.
#'
#' @return `AddTip()` returns a tree of class `phylo` with an additional tip
#' at the desired location.
#'
#' @template MRS
#'
#' @seealso Add one tree to another: \code{\link{bind.tree}()}
#'
#' @examples
#' tree <- BalancedTree(10)
#' 
#' # Add a leaf below an internal node
#' plot(tree)
#' ape::nodelabels()
#' node <- 15
#' ape::nodelabels(bg = ifelse(NodeNumbers(tree) == node, "green", "grey"))
#'
#' plot(AddTip(tree, 15, "NEW_TIP"))
#' 
#' # Add edge lengths for an ultrametric tree
#' tree$edge.length <- rep(c(rep(1, 5), 2, 1, 2, 2), 2)
#' 
#' # Add a leaf to an external edge
#' leaf <- 5
#' plot(tree)
#' ape::tiplabels(bg = ifelse(seq_len(NTip(tree)) == leaf, "green", "grey"))
#' 
#' plot(AddTip(tree, 5, "NEW_TIP", edgeLength = NULL))
#' 
#' @keywords tree
#' @family tree manipulation
#' @export
AddTip <- function(tree,
                   where = sample.int(tree[["Nnode"]] * 2 + 2L, size = 1) - 1L,
                   label = "New tip",
                   nodeLabel = "",
                   edgeLength = 0,
                   lengthBelow = NULL,
                   nTip = NTip(tree),
                   nNode = tree[["Nnode"]],
                   rootNode = RootNode(tree)
) {
  newTipNumber <- nTip + 1L
  treeEdge <- tree[["edge"]]
  edgeLengths <- tree[["edge.length"]]
  lengths <- !is.null(edgeLengths)
  
  if (is.character(where)) {
    tmp <- match(where, TipLabels(tree))
    if (is.na(tmp)) {
      stop("No tip labelled '", where, "'")
    }
    where <- tmp
  }
  ## find the row of "where" before renumbering
  if (where < 1L || where == rootNode) {
    case <- 1L
  } else {
    insertionEdge <- which(treeEdge[, 2] == where)
    case <- if (where <= nTip) 2L else 3L
  }
  # case = 1 -> y is bound on the root of x
  # case = 2 -> y is bound on a tip of x
  # case = 3 -> y is bound on a node of x
  
  # Because in all situations internal nodes need to be
  # renumbered, they are changed to negatives first, and
  # nodes eventually added will be numbered sequentially
  nodes <- treeEdge > nTip
  treeEdge[nodes] <- nTip - treeEdge[nodes]  # -1, ..., -nTip
  nextNode <- -nNode - 1L
  rootNode <- nTip - rootNode
  
  switch(case, { # case = 1 -> y is bound on the root of x
    treeEdge <- rbind(c(nextNode, treeEdge[1]), treeEdge, c(nextNode, newTipNumber))
    if (lengths) {
      if (is.null(lengthBelow)) {
        lengthBelow <- 0
      }
      edgeLengths <- c(lengthBelow, edgeLengths, 
                       if(is.null(edgeLength)) lengthBelow else edgeLength)
    }
    rootNode <- nextNode
  }, { # case = 2 -> y is bound on a tip of x
    beforeInsertion <- seq_len(insertionEdge)
    treeEdge[insertionEdge, 2] <- nextNode
    treeEdge <- rbind(treeEdge[beforeInsertion, ],
                      c(nextNode, where),
                      c(nextNode, newTipNumber),
                      treeEdge[-beforeInsertion, ])
    if (lengths) {
      if (is.null(lengthBelow)) {
        lengthBelow <- edgeLengths[insertionEdge] / 2L
      }
      edgeLengths <- c(edgeLengths[beforeInsertion[-insertionEdge]],
                       edgeLengths[insertionEdge] - lengthBelow,
                       lengthBelow,
                       if(is.null(edgeLength)) lengthBelow else edgeLength,
                       edgeLengths[-beforeInsertion])
    }
  }, { # case = 3 -> y is bound on a node of x
    beforeInsertion <- seq_len(insertionEdge)
    
    treeEdge <- rbind(treeEdge[beforeInsertion, ],
                      c(nextNode, newTipNumber),
                      c(nextNode, treeEdge[insertionEdge, 2]),
                      treeEdge[-beforeInsertion, ])
    treeEdge[insertionEdge, 2] <- nextNode
    
    if (lengths) {
      if (is.null(lengthBelow)) {
        lengthBelow <- edgeLengths[insertionEdge] / 2L
      }
      edgeLengths <- c(edgeLengths[beforeInsertion[-insertionEdge]],
                       edgeLengths[insertionEdge] - lengthBelow,
                       if(is.null(edgeLength)) lengthBelow else edgeLength,
                       lengthBelow,
                       edgeLengths[-beforeInsertion])
    }
    
  }
  )
  tree[["tip.label"]] <- c(tree[["tip.label"]], label)
  
  nNode <- nNode + 1L
  tree[["Nnode"]] <- nNode
  
  ## renumber nodes:
  newNumbering <- integer(nNode)
  newNumbering[-rootNode] <- newTipNumber + 1L
  childNodes <- treeEdge[, 2] < 0L
  
  ## executed from right to left, so newNb is modified before x$edge:
  treeEdge[childNodes, 2] <-
    newNumbering[-treeEdge[childNodes, 2]] <-
    newTipNumber + 2:nNode
  treeEdge[, 1] <- newNumbering[-treeEdge[, 1]]
  
  tree[["edge"]] <- treeEdge
  if (lengths) {
    tree[["edge.length"]] <- edgeLengths
  }
  
  nodeLabels <- tree[["node.label"]]
  if (!is.null(nodeLabels)) {
    newLabels <- character(nNode)
    newLabels[newNumbering - newTipNumber] <- c(nodeLabels, "")
    tree[["node.label"]] <- newLabels
  }
  
  # Return:
  tree
}

#' @rdname AddTip
#'
#' @details `AddTipEverywhere()` adds a tip to each edge in turn.
#'
#' @param includeRoot Logical; if `TRUE`, each position adjacent
#' to the root edge is considered to represent distinct edges; if `FALSE`,
#' they are treated as a single edge.
#' @return `AddTipEverywhere()` returns a list of class `multiPhylo` containing
#' the trees produced by adding `label` to each edge of `tree` in turn.
#'
#' @examples
#' # Set up multi-panel plot
#' oldPar <- par(mfrow = c(2, 4), mar = rep(0.3, 4), cex = 0.9)
#'
#' # Add leaf to each edge on a tree in turn
#' backbone <- BalancedTree(4)
#' # Treating the position of the root as instructive:
#' additions <- AddTipEverywhere(backbone, includeRoot = TRUE)
#' xx <- lapply(additions, plot)
#'
#' par(mfrow = c(2, 3))
#' # Don't treat root edges as distinct:
#' additions <- AddTipEverywhere(backbone, includeRoot = FALSE)
#' xx <- lapply(additions, plot)
#'
#' # Restore original plotting parameters
#' par(oldPar)
#'
#' @importFrom ape is.rooted
#' @export
AddTipEverywhere <- function(tree, label = "New tip", includeRoot = FALSE) {
  nTip <- NTip(tree)
  if (nTip == 0L) return(list(SingleTaxonTree(label)))
  if (nTip == 1L) return(list(StarTree(c(tree[["tip.label"]], label))))
  whichNodes <- seq_len(nTip + tree[["Nnode"]])
  edge <- tree[["edge"]]
  root <- RootNode(edge)
  if (!includeRoot) {
    parent <- edge[, 1]
    child <- edge[, 2]
    rootChildren <- child[parent == root]
    
    whichNodes <- if (length(rootChildren) == 2L && nTip > 2L) {
      rootChildrenNodes <- rootChildren[rootChildren > nTip]
      whichNodes[-c(root, rootChildrenNodes[1])]
    } else {
      whichNodes[-root]
    }
  }
  lapply(whichNodes, AddTip, tree = tree, label = label, nTip = nTip,
         rootNode = root)
}