1 ## drop.tip.R (2009-09-09)
3 ## Remove Tips in a Phylogenetic Tree
5 ## Copyright 2003-2009 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)
12 Ntip <- length(phy$tip.label)
14 Nedge <- dim(phy$edge)[1]
15 wbl <- !is.null(phy$edge.length)
16 if (length(node) > 1) {
18 warning("only the first value of 'node' has been considered")
20 if (is.character(node)) {
21 if (is.null(phy$node.label))
22 stop("the tree has no node labels")
23 node <- which(phy$node.label %in% node) + Ntip
26 stop("node number must be greater than the number of tips")
27 if (node == ROOT) return(phy)
28 phy <- reorder(phy) # insure it is in cladewise order
29 root.node <- which(phy$edge[, 2] == node)
30 start <- root.node + 1 # start of the clade looked for
31 anc <- phy$edge[root.node, 1] # the ancestor of 'node'
32 next.anc <- which(phy$edge[-(1:start), 1] <= anc) # find the next occurence of 'anc' or an 'older' node
34 keep <- if (length(next.anc)) start + 0:(next.anc[1] - 1) else start:Nedge
37 NewRootEdge <- phy$edge.length[root.node]
38 root.edge <- root.edge - 1
40 if (anc == ROOT) break
41 i <- which(phy$edge[, 2] == anc)
42 NewRootEdge <- NewRootEdge + phy$edge.length[i]
43 root.edge <- root.edge - 1
46 if (root.edge && !is.null(phy$root.edge))
47 NewRootEdge <- NewRootEdge + phy$root.edge
48 phy$root.edge <- NewRootEdge
51 phy$edge <- phy$edge[keep, ]
52 if (wbl) phy$edge.length <- phy$edge.length[keep]
53 TIPS <- phy$edge[, 2] <= Ntip
54 tip <- phy$edge[TIPS, 2]
55 phy$tip.label <- phy$tip.label[tip]
56 ## keep the ordering so no need to reorder tip.label:
57 phy$edge[TIPS, 2] <- order(tip)
58 if (!is.null(phy$node.label))
59 phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
60 Ntip <- length(phy$tip.label)
61 phy$Nnode <- dim(phy$edge)[1] - Ntip + 1L
62 ## The block below renumbers the nodes so that they conform
63 ## to the "phylo" format -- same as in root()
64 newNb <- integer(Ntip + phy$Nnode)
65 newNb[node] <- Ntip + 1L
66 sndcol <- phy$edge[, 2] > Ntip
67 ## executed from right to left, so newNb is modified before phy$edge:
68 phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
69 (Ntip + 2):(Ntip + phy$Nnode)
70 phy$edge[, 1] <- newNb[phy$edge[, 1]]
75 function(phy, tip, trim.internal = TRUE, subtree = FALSE,
76 root.edge = 0, rooted = is.rooted(phy))
78 if (!inherits(phy, "phylo"))
79 stop('object "phy" is not of class "phylo"')
81 Ntip <- length(phy$tip.label)
82 ## find the tips to drop:
83 if (is.character(tip))
84 tip <- which(phy$tip.label %in% tip)
86 if (!rooted && subtree) {
87 phy <- root(phy, (1:Ntip)[-tip][1])
92 NEWROOT <- ROOT <- Ntip + 1
94 Nedge <- dim(phy$edge)[1]
97 tr <- reorder(phy, "pruningwise")
98 N <- .C("node_depth", as.integer(Ntip), as.integer(Nnode),
99 as.integer(tr$edge[, 1]), as.integer(tr$edge[, 2]),
100 as.integer(Nedge), double(Ntip + Nnode),
101 DUP = FALSE, PACKAGE = "ape")[[6]]
103 wbl <- !is.null(phy$edge.length)
104 edge1 <- phy$edge[, 1] # local copies
105 edge2 <- phy$edge[, 2] #
106 keep <- !logical(Nedge)
108 ## find the tips to drop:
109 if (is.character(tip))
110 tip <- which(phy$tip.label %in% tip)
112 if (!rooted && subtree) {
113 phy <- root(phy, (1:Ntip)[-tip][1])
117 ## delete the terminal edges given by `tip':
118 keep[match(tip, edge2)] <- FALSE
122 ## delete the internal edges that do not have anymore
123 ## descendants (ie, they are in the 2nd col of `edge' but
124 ## not in the 1st one)
126 sel <- !(edge2 %in% edge1[keep]) & ints & keep
131 ## keep the subtending edge(s):
132 subt <- edge1 %in% edge1[keep] & edge1 %in% edge1[!keep]
135 if (root.edge && wbl) {
136 degree <- tabulate(edge1[keep])
137 if (degree[ROOT] == 1) {
138 j <- integer(0) # will store the indices of the edges below the new root
140 i <- which(edge1 == NEWROOT & keep)
143 degree <- tabulate(edge1[keep])
144 if (degree[NEWROOT] > 1) break
147 if (length(j) > root.edge) j <- 1:root.edge
148 NewRootEdge <- sum(phy$edge.length[j])
149 if (length(j) < root.edge && !is.null(phy$root.edge))
150 NewRootEdge <- NewRootEdge + phy$root.edge
151 phy$root.edge <- NewRootEdge
156 if (!root.edge) phy$root.edge <- NULL
159 phy$edge <- phy$edge[keep, ]
160 if (wbl) phy$edge.length <- phy$edge.length[keep]
162 ## find the new terminal edges (works whatever 'subtree' and 'trim.internal'):
163 TERMS <- !(phy$edge[, 2] %in% phy$edge[, 1])
165 ## get the old No. of the nodes and tips that become tips:
166 oldNo.ofNewTips <- phy$edge[TERMS, 2]
168 n <- length(oldNo.ofNewTips) # the new number of tips in the tree
170 ## the tips may not be sorted in increasing order of their
171 ## in the 2nd col of edge, so no need to reorder $tip.label
172 phy$edge[TERMS, 2] <- rank(phy$edge[TERMS, 2])
174 ## make new tip labels if necessary:
175 if (subtree || !trim.internal) {
176 ## get the logical indices of the tips that are kept within 'oldNo.ofNewTips':
177 tips.kept <- oldNo.ofNewTips <= Ntip & !(oldNo.ofNewTips %in% tip)
178 new.tip.label <- character(n)
179 new.tip.label[tips.kept] <- phy$tip.label[-tip]
180 ## get the numbers of the nodes that become tips:
181 node2tip <- oldNo.ofNewTips[!tips.kept]
182 new.tip.label[!tips.kept] <- if (subtree) {
183 paste("[", N[node2tip], "_tips]", sep = "")
185 if (is.null(phy$node.label)) rep("NA", length(node2tip))
186 else phy$node.label[node2tip - Ntip]
188 if (!is.null(phy$node.label))
189 phy$node.label <- phy$node.label[-(node2tip - Ntip)]
190 phy$tip.label <- new.tip.label
191 } else phy$tip.label <- phy$tip.label[-tip]
193 ## update node.label if needed:
194 if (!is.null(phy$node.label))
195 phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
197 phy$Nnode <- dim(phy$edge)[1] - n + 1L # update phy$Nnode
199 ## The block below renumbers the nodes so that they conform
200 ## to the "phylo" format -- same as in root()
201 newNb <- integer(n + phy$Nnode)
202 newNb[NEWROOT] <- n + 1L
203 sndcol <- phy$edge[, 2] > n
204 ## executed from right to left, so newNb is modified before phy$edge:
205 phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
206 (n + 2):(n + phy$Nnode)
207 phy$edge[, 1] <- newNb[phy$edge[, 1]]
208 storage.mode(phy$edge) <- "integer"
209 collapse.singles(phy)