R/drop.tip.R

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