1 ## drop.tip.R (2007-12-21)
3 ## Remove Tips in a Phylogenetic Tree
5 ## Copyright 2003-2007 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 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)))
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]
36 ## find the MRCA of the remaining tips:
38 ## This is modified since some tips were deleted!!
39 for (i in phy$edge[, 2][as.numeric(phy$edge[, 2]) > 0]) {
43 ind <- which(phy$edge[, 2] == j)
49 sn <- lapply(seq.nod, rev)
51 x <- unlist(lapply(sn, function(x) x[i]))
52 while (length(unique(x)) == 1) {
53 x <- unlist(lapply(sn, function(x) x[i]))
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
67 phy$root.edge <- newrootedge
69 if (!is.null(phy$root.edge)) phy$root.edge <- NULL
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]
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]
87 phy$tip.label <- c(phy$tip.label, new.tip.label)
89 useless.nodes <- names(which(table(phy$edge[, 1]) == 1))
91 if (!nobr) mnbr <- mean(phy$edge.length)
92 if (length(useless.nodes) == 1) n <- length(tip) else {
95 for (i in as.character(which(del))) { # it is not needed to loop through all tips!
98 while (!(j %in% useless.nodes)) {
99 ind <- which(edge.bak[, 2] == j)
101 j <- edge.bak[ind, 1]
106 n <- table(unlist(lapply(seq.nod, function(x) rev(x)[1])))
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)
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)])
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, ]
130 phy$edge.length[ind2] <- phy$edge.length[ind2] + phy$edge.length[ind1]
131 phy$edge.length <- phy$edge.length[-ind1]
135 tmp <- as.numeric(phy$edge)
136 if (!is.null(phy$node.label)) {
139 phy$node.label <- phy$node.label[-x]
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"
152 phy <- old2new.phylo(phy)
153 if (!trim.internal || subtree) {
155 phy <- if (nobr) clado.build(S) else tree.build(S)