1 ## drop.tip.R (2012-11-29)
3 ## Remove Tips in a Phylogenetic Tree
5 ## Copyright 2003-2012 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 extract.clade <- function(phy, node, root.edge = 0, interactive = FALSE)
12 Ntip <- length(phy$tip.label)
14 Nedge <- dim(phy$edge)[1]
15 wbl <- !is.null(phy$edge.length)
16 if (interactive) node <- identify(phy)$nodes else {
17 if (length(node) > 1) {
19 warning("only the first value of 'node' has been considered")
21 if (is.character(node)) {
22 if (is.null(phy$node.label))
23 stop("the tree has no node labels")
24 node <- which(phy$node.label %in% node) + Ntip
27 stop("node number must be greater than the number of tips")
29 if (node == ROOT) return(phy)
30 phy <- reorder(phy) # insure it is in cladewise order
31 root.node <- which(phy$edge[, 2] == node)
32 start <- root.node + 1 # start of the clade looked for
33 anc <- phy$edge[root.node, 1] # the ancestor of 'node'
34 next.anc <- which(phy$edge[-(1:start), 1] <= anc) # find the next occurence of 'anc' or an 'older' node
36 keep <- if (length(next.anc)) start + 0:(next.anc[1] - 1) else start:Nedge
39 NewRootEdge <- phy$edge.length[root.node]
40 root.edge <- root.edge - 1
42 if (anc == ROOT) break
43 i <- which(phy$edge[, 2] == anc)
44 NewRootEdge <- NewRootEdge + phy$edge.length[i]
45 root.edge <- root.edge - 1
48 if (root.edge && !is.null(phy$root.edge))
49 NewRootEdge <- NewRootEdge + phy$root.edge
50 phy$root.edge <- NewRootEdge
53 phy$edge <- phy$edge[keep, ]
54 if (wbl) phy$edge.length <- phy$edge.length[keep]
55 TIPS <- phy$edge[, 2] <= Ntip
56 tip <- phy$edge[TIPS, 2]
57 ## Fix by Ludovic Mallet and Mahendra Mariadassou (2011-11-21):
58 name <- vector("character", length(tip))
59 name[order(tip)] <- phy$tip.label[tip]
62 ## keep the ordering so no need to reorder tip.label:
63 phy$edge[TIPS, 2] <- order(tip)
64 if (!is.null(phy$node.label))
65 phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
66 Ntip <- length(phy$tip.label)
67 phy$Nnode <- dim(phy$edge)[1] - Ntip + 1L
68 ## The block below renumbers the nodes so that they conform
69 ## to the "phylo" format -- same as in root()
70 newNb <- integer(Ntip + phy$Nnode)
71 newNb[node] <- Ntip + 1L
72 sndcol <- phy$edge[, 2] > Ntip
73 ## executed from right to left, so newNb is modified before phy$edge:
74 phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
75 (Ntip + 2):(Ntip + phy$Nnode)
76 phy$edge[, 1] <- newNb[phy$edge[, 1]]
81 function(phy, tip, trim.internal = TRUE, subtree = FALSE,
82 root.edge = 0, rooted = is.rooted(phy), interactive = FALSE)
84 if (!inherits(phy, "phylo"))
85 stop('object "phy" is not of class "phylo"')
86 if (!length(tip)) return(phy)
88 Ntip <- length(phy$tip.label)
89 ## find the tips to drop:
91 cat("Left-click close to the tips you want to drop; right-click when finished...\n")
93 nToDrop <- length(xy$x)
94 tip <- integer(nToDrop)
95 lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
96 for (i in 1:nToDrop) {
97 d <- sqrt((xy$x[i] - lastPP$xx)^2 + (xy$y[i] - lastPP$yy)^2)
98 tip[i] <- which.min(d)
101 if (is.character(tip))
102 tip <- which(phy$tip.label %in% tip)
104 if (!length(tip)) return(phy)
106 warning("some tip numbers were higher than the number of tips")
108 if (!rooted && subtree) {
109 phy <- root(phy, (1:Ntip)[-tip][1])
114 NEWROOT <- ROOT <- Ntip + 1
116 Nedge <- dim(phy$edge)[1]
118 trim.internal <- TRUE
119 tr <- reorder(phy, "pruningwise")
120 N <- .C("node_depth", as.integer(Ntip), as.integer(Nnode),
121 as.integer(tr$edge[, 1]), as.integer(tr$edge[, 2]),
122 as.integer(Nedge), double(Ntip + Nnode),
123 DUP = FALSE, PACKAGE = "ape")[[6]]
125 wbl <- !is.null(phy$edge.length)
126 edge1 <- phy$edge[, 1] # local copies
127 edge2 <- phy$edge[, 2] #
128 keep <- !logical(Nedge)
130 ## delete the terminal edges given by `tip':
131 keep[match(tip, edge2)] <- FALSE
135 ## delete the internal edges that do not have anymore
136 ## descendants (ie, they are in the 2nd col of `edge' but
137 ## not in the 1st one)
139 sel <- !(edge2 %in% edge1[keep]) & ints & keep
144 ## keep the subtending edge(s):
145 subt <- edge1 %in% edge1[keep] & edge1 %in% edge1[!keep]
148 if (root.edge && wbl) {
149 degree <- tabulate(edge1[keep])
150 if (degree[ROOT] == 1) {
151 j <- integer(0) # will store the indices of the edges below the new root
153 i <- which(edge1 == NEWROOT & keep)
156 degree <- tabulate(edge1[keep])
157 if (degree[NEWROOT] > 1) break
160 if (length(j) > root.edge) j <- 1:root.edge
161 NewRootEdge <- sum(phy$edge.length[j])
162 if (length(j) < root.edge && !is.null(phy$root.edge))
163 NewRootEdge <- NewRootEdge + phy$root.edge
164 phy$root.edge <- NewRootEdge
169 if (!root.edge) phy$root.edge <- NULL
172 phy$edge <- phy$edge[keep, ]
173 if (wbl) phy$edge.length <- phy$edge.length[keep]
175 ## find the new terminal edges (works whatever 'subtree' and 'trim.internal'):
176 TERMS <- !(phy$edge[, 2] %in% phy$edge[, 1])
178 ## get the old No. of the nodes and tips that become tips:
179 oldNo.ofNewTips <- phy$edge[TERMS, 2]
181 ## in case some tips are dropped but kept because of 'subtree = TRUE':
183 i <- which(tip %in% oldNo.ofNewTips)
185 phy$tip.label[tip[i]] <- "[1_tip]"
190 n <- length(oldNo.ofNewTips) # the new number of tips in the tree
192 ## the tips may not be sorted in increasing order in the
193 ## 2nd col of edge, so no need to reorder $tip.label
194 phy$edge[TERMS, 2] <- rank(phy$edge[TERMS, 2])
195 phy$tip.label <- phy$tip.label[-tip]
197 ## make new tip labels if necessary:
198 if (subtree || !trim.internal) {
199 ## get the numbers of the nodes that become tips:
200 node2tip <- oldNo.ofNewTips[oldNo.ofNewTips > Ntip]
201 new.tip.label <- if (subtree) {
202 paste("[", N[node2tip], "_tips]", sep = "")
204 if (is.null(phy$node.label)) rep("NA", length(node2tip))
205 else phy$node.label[node2tip - Ntip]
207 # if (!is.null(phy$node.label))
208 # phy$node.label <- phy$node.label[-(node2tip - Ntip)]
209 phy$tip.label <- c(phy$tip.label, new.tip.label)
212 phy$Nnode <- dim(phy$edge)[1] - n + 1L # update phy$Nnode
214 ## The block below renumbers the nodes so that they conform
215 ## to the "phylo" format, same as in root()
216 newNb <- integer(Ntip + Nnode)
217 newNb[NEWROOT] <- n + 1L
218 sndcol <- phy$edge[, 2] > n
219 ## executed from right to left, so newNb is modified before phy$edge:
220 phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
221 (n + 2):(n + phy$Nnode)
222 phy$edge[, 1] <- newNb[phy$edge[, 1]]
223 storage.mode(phy$edge) <- "integer"
224 if (!is.null(phy$node.label)) # update node.label if needed
225 phy$node.label <- phy$node.label[which(newNb > 0) - Ntip]
226 collapse.singles(phy)