1 ## drop.tip.R (2008-04-17)
3 ## Remove Tips in a Phylogenetic Tree
5 ## Copyright 2003-2008 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 drop.tip <- function(phy, tip, trim.internal = TRUE, subtree = FALSE,
13 if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
14 phy <- new2old.phylo(phy)
19 tmp <- as.numeric(phy$edge)
22 nodes <- setdiff(tmp,1:nb.tip) #not sure if this also needs sorting into order
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)))
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]
38 ## find the MRCA of the remaining tips:
40 ## This is modified since some tips were deleted!!
41 for (i in phy$edge[, 2][as.numeric(phy$edge[, 2]) > 0]) {
45 ind <- which(phy$edge[, 2] == j)
51 sn <- lapply(seq.nod, rev)
53 x <- unlist(lapply(sn, function(x) x[i]))
54 while (length(unique(x)) == 1) {
55 x <- unlist(lapply(sn, function(x) x[i]))
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
69 phy$root.edge <- newrootedge
71 if (!is.null(phy$root.edge)) phy$root.edge <- NULL
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]
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
91 useless.nodes <- names(which(table(phy$edge[, 1]) == 1))
93 if (!nobr) mnbr <- mean(phy$edge.length)
94 if (length(useless.nodes) == 1) n <- length(tip) else {
97 for (i in as.character(which(del))) { # it is not needed to loop through all tips!
100 while (!(j %in% useless.nodes)) {
101 ind <- which(edge.bak[, 2] == j)
103 j <- edge.bak[ind, 1]
108 n <- table(unlist(lapply(seq.nod, function(x) rev(x)[1])))
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)
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)])
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, ]
132 phy$edge.length[ind2] <- phy$edge.length[ind2] + phy$edge.length[ind1]
133 phy$edge.length <- phy$edge.length[-ind1]
137 tmp <- as.numeric(phy$edge)
138 if (!is.null(phy$node.label)) {
141 phy$node.label <- phy$node.label[-x]
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"
154 phy <- old2new.phylo(phy)
155 if (!trim.internal || subtree) {
157 phy <- if (nobr) clado.build(S) else tree.build(S)