X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Fdrop.tip.R;h=5b635f680f4459a579242dfc4fdf2c549d9b5f22;hb=46bc16a7531d871fccb6ea7c2a31d980aac21bd1;hp=69c957629e6533598488d0f7b3c4ede773faede7;hpb=0e8dcca3046844b0b36b4f5855e6c10319db45a2;p=ape.git diff --git a/R/drop.tip.R b/R/drop.tip.R index 69c9576..5b635f6 100644 --- a/R/drop.tip.R +++ b/R/drop.tip.R @@ -1,35 +1,37 @@ -## drop.tip.R (2009-03-04) +## drop.tip.R (2011-11-21) ## Remove Tips in a Phylogenetic Tree -## Copyright 2003-2009 Emmanuel Paradis +## Copyright 2003-2011 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. -extract.clade <- function(phy, node, root.edge = 0) +extract.clade <- function(phy, node, root.edge = 0, interactive = FALSE) { Ntip <- length(phy$tip.label) ROOT <- Ntip + 1 Nedge <- dim(phy$edge)[1] wbl <- !is.null(phy$edge.length) - if (length(node) > 1) { - node <- node[1] - warning("only the first value of 'node' has been considered") - } - if (is.character(node)) { - if (is.null(phy$node.label)) - stop("the tree has no node labels") - node <- which(phy$node.label %in% node) + Ntip + if (interactive) node <- identify(phy)$nodes else { + if (length(node) > 1) { + node <- node[1] + warning("only the first value of 'node' has been considered") + } + if (is.character(node)) { + if (is.null(phy$node.label)) + stop("the tree has no node labels") + node <- which(phy$node.label %in% node) + Ntip + } + if (node <= Ntip) + stop("node number must be greater than the number of tips") } - if (node <= Ntip) - stop("node number must be greater than the number of tips") if (node == ROOT) return(phy) phy <- reorder(phy) # insure it is in cladewise order root.node <- which(phy$edge[, 2] == node) start <- root.node + 1 # start of the clade looked for anc <- phy$edge[root.node, 1] # the ancestor of 'node' - next.anc <- which(phy$edge[-(1:start), 1] == anc) # find the next occurence of 'anc' + next.anc <- which(phy$edge[-(1:start), 1] <= anc) # find the next occurence of 'anc' or an 'older' node keep <- if (length(next.anc)) start + 0:(next.anc[1] - 1) else start:Nedge @@ -52,7 +54,11 @@ extract.clade <- function(phy, node, root.edge = 0) if (wbl) phy$edge.length <- phy$edge.length[keep] TIPS <- phy$edge[, 2] <= Ntip tip <- phy$edge[TIPS, 2] - phy$tip.label <- phy$tip.label[tip] + ## Fix by Ludovic Mallet and Mahendra Mariadassou (2011-11-21): + name <- vector("character", length(tip)) + name[order(tip)] <- phy$tip.label[tip] + phy$tip.label <- name + ## End of fix ## keep the ordering so no need to reorder tip.label: phy$edge[TIPS, 2] <- order(tip) if (!is.null(phy$node.label)) @@ -72,11 +78,38 @@ extract.clade <- function(phy, node, root.edge = 0) } drop.tip <- - function(phy, tip, trim.internal = TRUE, subtree = FALSE, root.edge = 0) + function(phy, tip, trim.internal = TRUE, subtree = FALSE, + root.edge = 0, rooted = is.rooted(phy), interactive = FALSE) { - if (class(phy) != "phylo") + if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"') + if (!length(tip)) return(phy) + Ntip <- length(phy$tip.label) + ## find the tips to drop: + if (interactive) { + cat("Left-click close to the tips you want to drop; right-click when finished...\n") + xy <- locator() + nToDrop <- length(xy$x) + tip <- integer(nToDrop) + lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv) + for (i in 1:nToDrop) { + d <- sqrt((xy$x[i] - lastPP$xx)^2 + (xy$y[i] - lastPP$yy)^2) + tip[i] <- which.min(d) + } + } else { + if (is.character(tip)) + tip <- which(phy$tip.label %in% tip) + } + if (any(tip > Ntip)) + warning("some tip numbers were higher than the number of tips") + + if (!rooted && subtree) { + phy <- root(phy, (1:Ntip)[-tip][1]) + root.edge <- 0 + } + + phy <- reorder(phy) NEWROOT <- ROOT <- Ntip + 1 Nnode <- phy$Nnode Nedge <- dim(phy$edge)[1] @@ -92,27 +125,24 @@ drop.tip <- edge1 <- phy$edge[, 1] # local copies edge2 <- phy$edge[, 2] # keep <- !logical(Nedge) - ## find the tips to drop: - if (is.character(tip)) - tip <- which(phy$tip.label %in% tip) - trms <- edge2 <= Ntip + ## delete the terminal edges given by `tip': keep[match(tip, edge2)] <- FALSE if (trim.internal) { - ## delete the internal edges that do not have descendants - ## anymore (ie, they are in the 2nd column of `edge' but + ints <- edge2 > Ntip + ## delete the internal edges that do not have anymore + ## descendants (ie, they are in the 2nd col of `edge' but ## not in the 1st one) repeat { - sel <- !(edge2 %in% edge1[keep]) & !trms & keep + sel <- !(edge2 %in% edge1[keep]) & ints & keep if (!sum(sel)) break keep[sel] <- FALSE } if (subtree) { ## keep the subtending edge(s): subt <- edge1 %in% edge1[keep] & edge1 %in% edge1[!keep] - ## 'if (... ' needed below? - if (any(subt)) keep[which(subt)] <- TRUE + keep[subt] <- TRUE } if (root.edge && wbl) { degree <- tabulate(edge1[keep]) @@ -137,46 +167,63 @@ drop.tip <- if (!root.edge) phy$root.edge <- NULL - ## upate the tree; 1) drop the edges and tip labels + ## drop the edges phy$edge <- phy$edge[keep, ] if (wbl) phy$edge.length <- phy$edge.length[keep] + + ## find the new terminal edges (works whatever 'subtree' and 'trim.internal'): + TERMS <- !(phy$edge[, 2] %in% phy$edge[, 1]) + + ## get the old No. of the nodes and tips that become tips: + oldNo.ofNewTips <- phy$edge[TERMS, 2] + + ## in case some tips are dropped but kept because of 'subtree = TRUE': + if (subtree) { + i <- which(tip %in% oldNo.ofNewTips) + if (length(i)) { + phy$tip.label[tip[i]] <- "[1_tip]" + tip <- tip[-i] + } + } + + n <- length(oldNo.ofNewTips) # the new number of tips in the tree + + ## the tips may not be sorted in increasing order in the + ## 2nd col of edge, so no need to reorder $tip.label + phy$edge[TERMS, 2] <- rank(phy$edge[TERMS, 2]) phy$tip.label <- phy$tip.label[-tip] - ## 2) renumber the remaining tips now - TIPS <- phy$edge[, 2] <= Ntip - ## keep the ordering so no need to reorder tip.label: - phy$edge[TIPS, 2] <- rank(phy$edge[TIPS, 2]) - Ntip <- length(phy$tip.label) # update Ntip - ## make new tip labels if necessary + ## make new tip labels if necessary: if (subtree || !trim.internal) { - new.trms <- !(phy$edge[, 2] %in% phy$edge[, 1]) & phy$edge[, 2] > Ntip - node2tip <- phy$edge[new.trms, 2] - if (subtree) - new.lab <- paste("[", N[node2tip], "_tips]", sep = "") - else { - new.lab <- - if (is.null(phy$node.label)) rep("NA", length(node2tip)) - else phy$node.label[node2tip - Ntip] + ## get the numbers of the nodes that become tips: + node2tip <- oldNo.ofNewTips[oldNo.ofNewTips > Ntip] + new.tip.label <- if (subtree) { + paste("[", N[node2tip], "_tips]", sep = "") + } else { + if (is.null(phy$node.label)) rep("NA", length(node2tip)) + else phy$node.label[node2tip - Ntip] } - ## change the #'s in the edge matrix - new.tip <- Ntip + 1:length(node2tip) - phy$edge[new.trms, 2] <- new.tip - phy$tip.label[new.tip] <- new.lab - Ntip <- length(phy$tip.label) if (!is.null(phy$node.label)) phy$node.label <- phy$node.label[-(node2tip - Ntip)] + phy$tip.label <- c(phy$tip.label, new.tip.label) } - phy$Nnode <- dim(phy$edge)[1] - Ntip + 1L # 3) update phy$Nnode + + ## update node.label if needed: + if (!is.null(phy$node.label)) + phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip] + + phy$Nnode <- dim(phy$edge)[1] - n + 1L # update phy$Nnode ## The block below renumbers the nodes so that they conform ## to the "phylo" format -- same as in root() - newNb <- integer(Ntip + phy$Nnode) - newNb[NEWROOT] <- Ntip + 1L - sndcol <- phy$edge[, 2] > Ntip + newNb <- integer(n + phy$Nnode) + newNb[NEWROOT] <- n + 1L + sndcol <- phy$edge[, 2] > n ## executed from right to left, so newNb is modified before phy$edge: phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <- - (Ntip + 2):(Ntip + phy$Nnode) + (n + 2):(n + phy$Nnode) phy$edge[, 1] <- newNb[phy$edge[, 1]] storage.mode(phy$edge) <- "integer" collapse.singles(phy) } +