R/drop.tip.R

   1 ## drop.tip.R (2012-11-22)
   2
   3 ##   Remove Tips in a Phylogenetic Tree
   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 extract.clade <- function(phy, node, root.edge = 0, interactive = FALSE)
  11 {
  12     Ntip <- length(phy$tip.label)
  13     ROOT <- Ntip + 1
  14     Nedge <- dim(phy$edge)[1]
  15     wbl <- !is.null(phy$edge.length)
  16     if (interactive) node <- identify(phy)$nodes else {
  17         if (length(node) > 1) {
  18             node <- node[1]
  19             warning("only the first value of 'node' has been considered")
  20         }
  21         if (is.character(node)) {
  22             if (is.null(phy$node.label))
  23                 stop("the tree has no node labels")
  24             node <- which(phy$node.label %in% node) + Ntip
  25         }
  26         if (node <= Ntip)
  27             stop("node number must be greater than the number of tips")
  28     }
  29     if (node == ROOT) return(phy)
  30     phy <- reorder(phy) # insure it is in cladewise order
  31     root.node <- which(phy$edge[, 2] == node)
  32     start <- root.node + 1 # start of the clade looked for
  33     anc <- phy$edge[root.node, 1] # the ancestor of 'node'
  34     next.anc <- which(phy$edge[-(1:start), 1] <= anc) # find the next occurence of 'anc' or an 'older' node
  35
  36     keep <- if (length(next.anc)) start + 0:(next.anc[1] - 1) else start:Nedge
  37
  38     if (root.edge) {
  39         NewRootEdge <- phy$edge.length[root.node]
  40         root.edge <- root.edge - 1
  41         while (root.edge) {
  42             if (anc == ROOT) break
  43             i <- which(phy$edge[, 2] ==  anc)
  44             NewRootEdge <- NewRootEdge + phy$edge.length[i]
  45             root.edge <- root.edge - 1
  46             anc <- phy$edge[i, 1]
  47         }
  48         if (root.edge && !is.null(phy$root.edge))
  49             NewRootEdge <- NewRootEdge + phy$root.edge
  50         phy$root.edge <- NewRootEdge
  51     }
  52
  53     phy$edge <- phy$edge[keep, ]
  54     if (wbl) phy$edge.length <- phy$edge.length[keep]
  55     TIPS <- phy$edge[, 2] <= Ntip
  56     tip <- phy$edge[TIPS, 2]
  57     ## Fix by Ludovic Mallet and Mahendra Mariadassou (2011-11-21):
  58     name <- vector("character", length(tip))
  59     name[order(tip)] <- phy$tip.label[tip]
  60     phy$tip.label <- name
  61     ## End of fix
  62     ## keep the ordering so no need to reorder tip.label:
  63     phy$edge[TIPS, 2] <- order(tip)
  64     if (!is.null(phy$node.label))
  65         phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
  66     Ntip <- length(phy$tip.label)
  67     phy$Nnode <- dim(phy$edge)[1] - Ntip + 1L
  68     ## The block below renumbers the nodes so that they conform
  69     ## to the "phylo" format -- same as in root()
  70     newNb <- integer(Ntip + phy$Nnode)
  71     newNb[node] <- Ntip + 1L
  72     sndcol <- phy$edge[, 2] > Ntip
  73     ## executed from right to left, so newNb is modified before phy$edge:
  74     phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
  75         (Ntip + 2):(Ntip + phy$Nnode)
  76     phy$edge[, 1] <- newNb[phy$edge[, 1]]
  77     phy
  78 }
  79
  80 drop.tip <-
  81     function(phy, tip, trim.internal = TRUE, subtree = FALSE,
  82              root.edge = 0, rooted = is.rooted(phy), interactive = FALSE)
  83 {
  84     if (!inherits(phy, "phylo"))
  85         stop('object "phy" is not of class "phylo"')
  86     if (!length(tip)) return(phy)
  87
  88     Ntip <- length(phy$tip.label)
  89     ## find the tips to drop:
  90     if (interactive) {
  91         cat("Left-click close to the tips you want to drop; right-click when finished...\n")
  92         xy <- locator()
  93         nToDrop <- length(xy$x)
  94         tip <- integer(nToDrop)
  95         lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
  96         for (i in 1:nToDrop) {
  97             d <- sqrt((xy$x[i] - lastPP$xx)^2 + (xy$y[i] - lastPP$yy)^2)
  98             tip[i] <- which.min(d)
  99         }
 100     } else {
 101         if (is.character(tip))
 102             tip <- which(phy$tip.label %in% tip)
 103     }
 104     if (!length(tip)) return(phy)
 105     if (any(tip > Ntip))
 106         warning("some tip numbers were higher than the number of tips")
 107
 108     if (!rooted && subtree) {
 109         phy <- root(phy, (1:Ntip)[-tip][1])
 110         root.edge <- 0
 111     }
 112
 113     phy <- reorder(phy)
 114     NEWROOT <- ROOT <- Ntip + 1
 115     Nnode <- phy$Nnode
 116     Nedge <- dim(phy$edge)[1]
 117     if (subtree) {
 118         trim.internal <- TRUE
 119         tr <- reorder(phy, "pruningwise")
 120         N <- .C("node_depth", as.integer(Ntip), as.integer(Nnode),
 121                 as.integer(tr$edge[, 1]), as.integer(tr$edge[, 2]),
 122                 as.integer(Nedge), double(Ntip + Nnode),
 123                 DUP = FALSE, PACKAGE = "ape")[[6]]
 124     }
 125     wbl <- !is.null(phy$edge.length)
 126     edge1 <- phy$edge[, 1] # local copies
 127     edge2 <- phy$edge[, 2] #
 128     keep <- !logical(Nedge)
 129
 130     ## delete the terminal edges given by `tip':
 131     keep[match(tip, edge2)] <- FALSE
 132
 133     if (trim.internal) {
 134         ints <- edge2 > Ntip
 135         ## delete the internal edges that do not have anymore
 136         ## descendants (ie, they are in the 2nd col of `edge' but
 137         ## not in the 1st one)
 138         repeat {
 139             sel <- !(edge2 %in% edge1[keep]) & ints & keep
 140             if (!sum(sel)) break
 141             keep[sel] <- FALSE
 142         }
 143         if (subtree) {
 144             ## keep the subtending edge(s):
 145             subt <- edge1 %in% edge1[keep] & edge1 %in% edge1[!keep]
 146             keep[subt] <- TRUE
 147         }
 148         if (root.edge && wbl) {
 149             degree <- tabulate(edge1[keep])
 150             if (degree[ROOT] == 1) {
 151                 j <- integer(0) # will store the indices of the edges below the new root
 152                 repeat {
 153                     i <- which(edge1 == NEWROOT & keep)
 154                     j <- c(i, j)
 155                     NEWROOT <- edge2[i]
 156                     degree <- tabulate(edge1[keep])
 157                     if (degree[NEWROOT] > 1) break
 158                 }
 159                 keep[j] <- FALSE
 160                 if (length(j) > root.edge) j <- 1:root.edge
 161                 NewRootEdge <- sum(phy$edge.length[j])
 162                 if (length(j) < root.edge && !is.null(phy$root.edge))
 163                     NewRootEdge <- NewRootEdge + phy$root.edge
 164                 phy$root.edge <- NewRootEdge
 165             }
 166         }
 167     }
 168
 169     if (!root.edge) phy$root.edge <- NULL
 170
 171     ## drop the edges
 172     phy$edge <- phy$edge[keep, ]
 173     if (wbl) phy$edge.length <- phy$edge.length[keep]
 174
 175     ## find the new terminal edges (works whatever 'subtree' and 'trim.internal'):
 176     TERMS <- !(phy$edge[, 2] %in% phy$edge[, 1])
 177
 178     ## get the old No. of the nodes and tips that become tips:
 179     oldNo.ofNewTips <- phy$edge[TERMS, 2]
 180
 181     ## in case some tips are dropped but kept because of 'subtree = TRUE':
 182     if (subtree) {
 183         i <- which(tip %in% oldNo.ofNewTips)
 184         if (length(i)) {
 185             phy$tip.label[tip[i]] <- "[1_tip]"
 186             tip <- tip[-i]
 187         }
 188     }
 189
 190     n <- length(oldNo.ofNewTips) # the new number of tips in the tree
 191
 192     ## the tips may not be sorted in increasing order in the
 193     ## 2nd col of edge, so no need to reorder $tip.label
 194     phy$edge[TERMS, 2] <- rank(phy$edge[TERMS, 2])
 195     phy$tip.label <- phy$tip.label[-tip]
 196
 197     ## make new tip labels if necessary:
 198     if (subtree || !trim.internal) {
 199         ## get the numbers of the nodes that become tips:
 200         node2tip <- oldNo.ofNewTips[oldNo.ofNewTips > Ntip]
 201         new.tip.label <- if (subtree) {
 202             paste("[", N[node2tip], "_tips]", sep = "")
 203         } else {
 204             if (is.null(phy$node.label)) rep("NA", length(node2tip))
 205             else phy$node.label[node2tip - Ntip]
 206         }
 207         if (!is.null(phy$node.label))
 208             phy$node.label <- phy$node.label[-(node2tip - Ntip)]
 209         phy$tip.label <- c(phy$tip.label, new.tip.label)
 210     }
 211
 212     phy$Nnode <- dim(phy$edge)[1] - n + 1L # update phy$Nnode
 213
 214     ## The block below renumbers the nodes so that they conform
 215     ## to the "phylo" format, same as in root()
 216     newNb <- integer(Ntip + Nnode)
 217     newNb[NEWROOT] <- n + 1L
 218     sndcol <- phy$edge[, 2] > n
 219     ## executed from right to left, so newNb is modified before phy$edge:
 220     phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
 221         (n + 2):(n + phy$Nnode)
 222     phy$edge[, 1] <- newNb[phy$edge[, 1]]
 223     storage.mode(phy$edge) <- "integer"
 224     if (!is.null(phy$node.label)) # update node.label if needed
 225         phy$node.label <- phy$node.label[which(newNb > 0) - Ntip]
 226     collapse.singles(phy)
 227 }
 228