R/drop.tip.R

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