R/drop.tip.R

   1 ## drop.tip.R (2008-04-17)
   2
   3 ##   Remove Tips in a Phylogenetic Tree
   4
   5 ## Copyright 2003-2008 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 drop.tip <- function(phy, tip, trim.internal = TRUE, subtree = FALSE,
  11                      root.edge = 0)
  12 {
  13     if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
  14     phy <- new2old.phylo(phy)
  15     if (subtree) {
  16         trim.internal <- TRUE
  17         edge.bak <- phy$edge
  18     }
  19     tmp <- as.numeric(phy$edge)
  20     nb.tip <- max(tmp)
  21     ## fix by Yan Wong:
  22     nodes <- setdiff(tmp,1:nb.tip) #not sure if this also needs sorting into order
  23     ## end
  24     nobr <- is.null(phy$edge.length)
  25     if (is.numeric(tip)) tip <- phy$tip.label[tip]
  26     ## find the tips to drop...:
  27     del <- phy$tip.label %in% tip
  28     ## ... and the corresponding terminal branches:
  29     ind <- which(phy$edge[, 2] %in% as.character(which(del)))
  30     ## drop them...:
  31     phy$edge <- phy$edge[-ind, ]
  32     ## ... and the lengths if applies:
  33     if (!nobr) phy$edge.length <- phy$edge.length[-ind]
  34     ## drop the tip labels:
  35     phy$tip.label <- phy$tip.label[!del]
  36     if (trim.internal) {
  37         if (root.edge) {
  38             ## find the MRCA of the remaining tips:
  39             seq.nod <- list()
  40             ## This is modified since some tips were deleted!!
  41             for (i in phy$edge[, 2][as.numeric(phy$edge[, 2]) > 0]) {
  42                 vec <- i
  43                 j <- i
  44                 while (j != "-1") {
  45                     ind <- which(phy$edge[, 2] == j)
  46                     j <- phy$edge[ind, 1]
  47                     vec <- c(vec, j)
  48                 }
  49                 seq.nod[[i]] <- vec
  50             }
  51             sn <- lapply(seq.nod, rev)
  52             i <- 1
  53             x <- unlist(lapply(sn, function(x) x[i]))
  54             while (length(unique(x)) == 1) {
  55                 x <- unlist(lapply(sn, function(x) x[i]))
  56                 i <-  i + 1
  57             }
  58             MRCA <- sn[[1]][i - 2]
  59             newrootedge <- if (is.null(phy$root.edge)) 0 else phy$root.edge
  60             for (i in 1:root.edge) {
  61                 ind <- which(phy$edge[, 2] == MRCA)
  62                 newrootedge <- newrootedge + phy$edge.length[ind]
  63                 MRCA <- phy$edge[ind, 1]
  64                 if (MRCA == "-1" && i < root.edge) {
  65                     newrootedge <- newrootedge
  66                     break
  67                 }
  68             }
  69             phy$root.edge <- newrootedge
  70         } else {
  71             if (!is.null(phy$root.edge)) phy$root.edge <- NULL
  72         }
  73         while (!all(phy$edge[, 2][as.numeric(phy$edge[, 2]) < 0] %in% phy$edge[, 1])) {
  74             temp <- phy$edge[, 2][as.numeric(phy$edge[, 2]) < 0]
  75             k <- temp %in% phy$edge[, 1]
  76             ind <- phy$edge[, 2] %in% temp[!k]
  77             phy$edge <- phy$edge[!ind, ]
  78             if (!nobr) phy$edge.length <- phy$edge.length[!ind]
  79         }
  80     } else {
  81         ## fix by Yan Wong:
  82         k <- nodes %in% phy$edge[, 1] #nodes that have descendants
  83         ind <- phy$edge[, 2] %in% nodes[!k]
  84         phy$edge[which(ind), 2] <- as.character(nb.tip + (1:sum(ind)))
  85         if (is.null(phy$node.label)) new.tip.label <- rep("NA", sum(ind))
  86         else new.tip.label <- phy$node.label[!k]
  87         phy$tip.label <- c(phy$tip.label, new.tip.label)
  88         #N.B. phy$node.label can be left: it is altered later
  89         ## end
  90     }
  91     useless.nodes <- names(which(table(phy$edge[, 1]) == 1))
  92     if (subtree) {
  93         if (!nobr) mnbr <- mean(phy$edge.length)
  94         if (length(useless.nodes) == 1) n <- length(tip) else {
  95             seq.nod <- list()
  96             wh <- numeric(0)
  97             for (i in as.character(which(del))) { # it is not needed to loop through all tips!
  98                 vec <- i
  99                 j <- i
 100                 while (!(j %in% useless.nodes)) {
 101                     ind <- which(edge.bak[, 2] == j)
 102                     wh <- c(wh, ind)
 103                     j <- edge.bak[ind, 1]
 104                     vec <- c(vec, j)
 105                 }
 106                 seq.nod[[i]] <- vec
 107             }
 108             n <- table(unlist(lapply(seq.nod, function(x) rev(x)[1])))
 109         }
 110         new.lab <- paste("[", n, "_tips]", sep = "")
 111         for (i in 1:length(useless.nodes)) {
 112             wh <- which(phy$edge[, 1] == useless.nodes[i])
 113             phy$tip.label <- c(phy$tip.label, new.lab[i])
 114             if (wh == dim(phy$edge)[1]) {
 115                 phy$edge <- rbind(phy$edge, c(useless.nodes[i], as.character(nb.tip + i)))
 116                 if (!nobr) phy$edge.length <- c(phy$edge.length, mnbr)
 117             } else {
 118                 phy$edge <- rbind(phy$edge[1:wh, ],
 119                                   c(useless.nodes[i], as.character(nb.tip + i)),
 120                                   phy$edge[(wh + 1):dim(phy$edge)[1], ])
 121                 if (!nobr) phy$edge.length <- c(phy$edge.length[1:wh], mnbr,
 122                                                 phy$edge.length[(wh + 1):(dim(phy$edge)[1] - 1)])
 123             }
 124         }
 125     } else {
 126         for (i in useless.nodes) {
 127             ind1 <- which(phy$edge[, 1] == i)
 128             ind2 <- which(phy$edge[, 2] == i)
 129             phy$edge[ind2, 2] <- phy$edge[ind1, 2]
 130             phy$edge <- phy$edge[-ind1, ]
 131             if (!nobr) {
 132                 phy$edge.length[ind2] <- phy$edge.length[ind2] + phy$edge.length[ind1]
 133                 phy$edge.length <- phy$edge.length[-ind1]
 134             }
 135         }
 136     }
 137     tmp <- as.numeric(phy$edge)
 138     if (!is.null(phy$node.label)) {
 139         x <- unique(tmp)
 140         x <- x[x < 0]
 141         phy$node.label <- phy$node.label[-x]
 142     }
 143     n <- length(tmp)
 144     nodes <- tmp < 0
 145     ind.nodes <- (1:n)[nodes]
 146     ind.tips <- (1:n)[!nodes]
 147     new.nodes <- -as.numeric(factor(-tmp[nodes]))
 148     new.tips <- as.numeric(factor(tmp[!nodes]))
 149     tmp[ind.nodes] <- new.nodes
 150     tmp[ind.tips] <- new.tips
 151     dim(tmp) <- c(n / 2, 2)
 152     mode(tmp) <- "character"
 153     phy$edge <- tmp
 154     phy <- old2new.phylo(phy)
 155     if (!trim.internal || subtree) {
 156         S <- write.tree(phy)
 157         phy <- if (nobr) clado.build(S) else tree.build(S)
 158     }
 159     phy
 160 }