R/drop.tip.R

   1 ## drop.tip.R (2009-03-22)
   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'
  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, root.edge = 0)
  76 {
  77     if (class(phy) != "phylo")
  78         stop('object "phy" is not of class "phylo"')
  79     phy <- reorder(phy)
  80     Ntip <- length(phy$tip.label)
  81     NEWROOT <- ROOT <- Ntip + 1
  82     Nnode <- phy$Nnode
  83     Nedge <- dim(phy$edge)[1]
  84     if (subtree) {
  85         trim.internal <- TRUE
  86         tr <- reorder(phy, "pruningwise")
  87         N <- .C("node_depth", as.integer(Ntip), as.integer(Nnode),
  88                 as.integer(tr$edge[, 1]), as.integer(tr$edge[, 2]),
  89                 as.integer(Nedge), double(Ntip + Nnode),
  90                 DUP = FALSE, PACKAGE = "ape")[[6]]
  91     }
  92     wbl <- !is.null(phy$edge.length)
  93     edge1 <- phy$edge[, 1] # local copies
  94     edge2 <- phy$edge[, 2] #
  95     keep <- !logical(Nedge)
  96     ## find the tips to drop:
  97     if (is.character(tip))
  98         tip <- which(phy$tip.label %in% tip)
  99     trms <- edge2 <= Ntip
 100     ## delete the terminal edges given by `tip':
 101     keep[match(tip, edge2)] <- FALSE
 102
 103     if (trim.internal) {
 104         ## delete the internal edges that do not have descendants
 105         ## anymore (ie, they are in the 2nd column of `edge' but
 106         ## not in the 1st one)
 107         repeat {
 108             sel <- !(edge2 %in% edge1[keep]) & !trms & keep
 109             if (!sum(sel)) break
 110             keep[sel] <- FALSE
 111         }
 112         if (subtree) {
 113             ## keep the subtending edge(s):
 114             subt <- edge1 %in% edge1[keep] & edge1 %in% edge1[!keep]
 115             ## <FIXME> 'if (... ' needed below?
 116             if (any(subt)) keep[which(subt)] <- TRUE
 117         }
 118         if (root.edge && wbl) {
 119             degree <- tabulate(edge1[keep])
 120             if (degree[ROOT] == 1) {
 121                 j <- integer(0) # will store the indices of the edges below the new root
 122                 repeat {
 123                     i <- which(edge1 == NEWROOT & keep)
 124                     j <- c(i, j)
 125                     NEWROOT <- edge2[i]
 126                     degree <- tabulate(edge1[keep])
 127                     if (degree[NEWROOT] > 1) break
 128                 }
 129                 keep[j] <- FALSE
 130                 if (length(j) > root.edge) j <- 1:root.edge
 131                 NewRootEdge <- sum(phy$edge.length[j])
 132                 if (length(j) < root.edge && !is.null(phy$root.edge))
 133                     NewRootEdge <- NewRootEdge + phy$root.edge
 134                 phy$root.edge <- NewRootEdge
 135             }
 136         }
 137     }
 138
 139     if (!root.edge) phy$root.edge <- NULL
 140
 141     ## upate the tree; 1) drop the edges and tip labels
 142     phy$edge <- phy$edge[keep, ]
 143     if (wbl) phy$edge.length <- phy$edge.length[keep]
 144     phy$tip.label <- phy$tip.label[-tip]
 145     ## 2) renumber the remaining tips now
 146     TIPS <- phy$edge[, 2] <= Ntip
 147     ## keep the ordering so no need to reorder tip.label:
 148     phy$edge[TIPS, 2] <- rank(phy$edge[TIPS, 2])
 149     ## 3) update node.label if needed
 150     if (!is.null(phy$node.label))
 151         phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
 152     Ntip <- length(phy$tip.label) # 4) update Ntip
 153
 154     ## make new tip labels if necessary
 155     if (subtree || !trim.internal) {
 156         new.trms <- !(phy$edge[, 2] %in% phy$edge[, 1]) & phy$edge[, 2] > Ntip
 157         node2tip <- phy$edge[new.trms, 2]
 158         if (subtree)
 159             new.lab <- paste("[", N[node2tip], "_tips]", sep = "")
 160         else {
 161             new.lab <-
 162               if (is.null(phy$node.label)) rep("NA", length(node2tip))
 163               else phy$node.label[node2tip - Ntip]
 164         }
 165         ## change the #'s in the edge matrix
 166         new.tip <- Ntip + 1:length(node2tip)
 167         phy$edge[new.trms, 2] <- new.tip
 168         phy$tip.label[new.tip] <- new.lab
 169         Ntip <- length(phy$tip.label)
 170         if (!is.null(phy$node.label))
 171             phy$node.label <- phy$node.label[-(node2tip - Ntip)]
 172     }
 173     phy$Nnode <- dim(phy$edge)[1] - Ntip + 1L # 3) update phy$Nnode
 174
 175     ## The block below renumbers the nodes so that they conform
 176     ## to the "phylo" format -- same as in root()
 177     newNb <- integer(Ntip + phy$Nnode)
 178     newNb[NEWROOT] <- Ntip + 1L
 179     sndcol <- phy$edge[, 2] > Ntip
 180     ## executed from right to left, so newNb is modified before phy$edge:
 181     phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
 182         (Ntip + 2):(Ntip + phy$Nnode)
 183     phy$edge[, 1] <- newNb[phy$edge[, 1]]
 184     storage.mode(phy$edge) <- "integer"
 185     collapse.singles(phy)
 186 }