1 ## drop.tip.R (2009-01-07)
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
25 if (node <= Ntip) stop("node number must be greater than the number of tips")
26 if (node == ROOT) return(phy)
27 phy <- reorder(phy) # insure it is in cladewise order
28 root.node <- which(phy$edge[, 2] == node)
29 start <- root.node + 1 # start of the clade looked for
30 anc <- phy$edge[root.node, 1] # the ancestor of 'node'
31 next.anc <- which(phy$edge[-(1:start), 1] == anc) # find the next occurence of 'anc'
33 keep <- if (length(next.anc)) start + 0:(next.anc[1] - 1) else start:Nedge
36 NewRootEdge <- phy$edge.length[root.node]
37 root.edge <- root.edge - 1
39 if (anc == ROOT) break
40 i <- which(phy$edge[, 2] == anc)
41 NewRootEdge <- NewRootEdge + phy$edge.length[i]
42 root.edge <- root.edge - 1
45 if (root.edge && !is.null(phy$root.edge))
46 NewRootEdge <- NewRootEdge + phy$root.edge
47 phy$root.edge <- NewRootEdge
50 phy$edge <- phy$edge[keep, ]
51 if (wbl) phy$edge.length <- phy$edge.length[keep]
52 TIPS <- phy$edge[, 2] <= Ntip
53 tip <- phy$edge[TIPS, 2]
54 phy$tip.label <- phy$tip.label[tip]
55 ## keep the ordering so no need to reorder tip.label:
56 phy$edge[TIPS, 2] <- order(tip)
57 if (!is.null(phy$node.label))
58 phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
59 Ntip <- length(phy$tip.label)
60 phy$Nnode <- dim(phy$edge)[1] - Ntip + 1L
61 ## The block below renumbers the nodes so that they conform
62 ## to the "phylo" format -- same as in root()
63 newNb <- integer(Ntip + phy$Nnode)
64 newNb[node] <- Ntip + 1L
65 sndcol <- phy$edge[, 2] > Ntip
66 ## executed from right to left, so newNb is modified before phy$edge:
67 phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
68 (Ntip + 2):(Ntip + phy$Nnode)
69 phy$edge[, 1] <- newNb[phy$edge[, 1]]
74 function(phy, tip, trim.internal = TRUE, subtree = FALSE, root.edge = 0)
76 if (class(phy) != "phylo")
77 stop('object "phy" is not of class "phylo"')
78 Ntip <- length(phy$tip.label)
79 NEWROOT <- ROOT <- Ntip + 1
81 Nedge <- dim(phy$edge)[1]
84 tr <- reorder(phy, "pruningwise")
85 N <- .C("node_depth", as.integer(Ntip), as.integer(Nnode),
86 as.integer(tr$edge[, 1]), as.integer(tr$edge[, 2]),
87 as.integer(Nedge), double(Ntip + Nnode),
88 DUP = FALSE, PACKAGE = "ape")[[6]]
90 wbl <- !is.null(phy$edge.length)
91 edge1 <- phy$edge[, 1] # local copies
92 edge2 <- phy$edge[, 2] #
93 keep <- !logical(Nedge)
94 ## find the tips to drop:
95 if (is.character(tip))
96 tip <- which(phy$tip.label %in% tip)
98 ## delete the terminal edges given by `tip':
99 keep[match(tip, edge2)] <- FALSE
102 ## delete the internal edges that do not have descendants
103 ## anymore (ie, they are in the 2nd column of `edge' but
104 ## not in the 1st one)
106 sel <- !(edge2 %in% edge1[keep]) & !trms & keep
111 ## keep the subtending edge(s):
112 subt <- edge1 %in% edge1[keep] & edge1 %in% edge1[!keep]
113 ## <FIXME> 'if (... ' needed below?
114 if (any(subt)) keep[which(subt)] <- TRUE
116 if (root.edge && wbl) {
117 degree <- tabulate(edge1[keep])
118 if (degree[ROOT] == 1) {
119 j <- integer(0) # will store the indices of the edges below the new root
121 i <- which(edge1 == NEWROOT & keep)
124 degree <- tabulate(edge1[keep])
125 if (degree[NEWROOT] > 1) break
128 if (length(j) > root.edge) j <- 1:root.edge
129 NewRootEdge <- sum(phy$edge.length[j])
130 if (length(j) < root.edge && !is.null(phy$root.edge))
131 NewRootEdge <- NewRootEdge + phy$root.edge
132 phy$root.edge <- NewRootEdge
137 if (!root.edge) phy$root.edge <- NULL
139 ## upate the tree; 1) drop the edges and tip labels
140 phy$edge <- phy$edge[keep, ]
141 if (wbl) phy$edge.length <- phy$edge.length[keep]
142 phy$tip.label <- phy$tip.label[-tip]
143 ## 2) renumber the remaining tips now
144 TIPS <- phy$edge[, 2] <= Ntip
145 ## keep the ordering so no need to reorder tip.label:
146 phy$edge[TIPS, 2] <- order(phy$edge[TIPS, 2])
147 Ntip <- length(phy$tip.label) # update Ntip
149 ## make new tip labels if necessary
150 if (subtree || !trim.internal) {
151 new.trms <- !(phy$edge[, 2] %in% phy$edge[, 1]) & phy$edge[, 2] > Ntip
152 node2tip <- phy$edge[new.trms, 2]
154 new.lab <- paste("[", N[node2tip], "_tips]", sep = "")
157 if (is.null(phy$node.label)) rep("NA", length(node2tip))
158 else phy$node.label[node2tip - Ntip]
160 ## change the #'s in the edge matrix
161 new.tip <- Ntip + 1:length(node2tip)
162 phy$edge[new.trms, 2] <- new.tip
163 phy$tip.label[new.tip] <- new.lab
164 Ntip <- length(phy$tip.label)
165 if (!is.null(phy$node.label))
166 phy$node.label <- phy$node.label[-(node2tip - Ntip)]
168 phy$Nnode <- dim(phy$edge)[1] - Ntip + 1L # 3) update phy$Nnode
170 ## The block below renumbers the nodes so that they conform
171 ## to the "phylo" format -- same as in root()
172 newNb <- integer(Ntip + phy$Nnode)
173 newNb[NEWROOT] <- Ntip + 1L
174 sndcol <- phy$edge[, 2] > Ntip
175 ## executed from right to left, so newNb is modified before phy$edge:
176 phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
177 (Ntip + 2):(Ntip + phy$Nnode)
178 phy$edge[, 1] <- newNb[phy$edge[, 1]]
180 collapse.singles(phy)