[Phylobase-commits] r611 - in pkg: R inst/unitTests man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Sep 2 09:10:15 CEST 2009
Author: regetz
Date: 2009-09-02 09:10:14 +0200 (Wed, 02 Sep 2009)
New Revision: 611
Added:
pkg/tests/trees.RData
Modified:
pkg/R/prune.R
pkg/inst/unitTests/runit.subset.R
pkg/man/prune-methods.Rd
Log:
Added new phylo4 prune method, replacing the ape:drop.tip wrapper for
both prune and subset for both phylo4 and phylo4d objects. Partially
updated documentation and tests.
Modified: pkg/R/prune.R
===================================================================
--- pkg/R/prune.R 2009-08-31 21:54:27 UTC (rev 610)
+++ pkg/R/prune.R 2009-09-02 07:10:14 UTC (rev 611)
@@ -30,23 +30,129 @@
}
setMethod("prune","phylo4",
- function(phy, tip, trim.internal = TRUE, subtree = FALSE,
- ...) {
- DropTip(phy,tip,trim.internal, subtree)
- })
+ function(phy, tip, trim.internal = TRUE, subtree = FALSE, ...) {
+ if (subtree) {
+ #stop("subtree option is not currently supported for phylo4")
+ # do this for now (at least to allow examples to pass check)
+ warning("subtree option is not currently supported for phylo4")
+ return(DropTip(phy,tip,trim.internal, subtree))
+ }
+
+ makeEdgeNames <- function(edge) {
+ paste(edge[,1], edge[,2], sep="-")
+ }
+
+ ## drop tips and obsolete internal nodes from edge matrix
+ tip.drop <- getNode(phy, tip, missing="fail")
+ tip.keep <- setdiff(seq_len(nTips(phy)), tip.drop)
+ nodes <- seq_len(nrow(phy at edge))
+ node.keep <- logical(nrow(phy at edge))
+ node.keep[tip.keep] <- TRUE
+ if (trim.internal) {
+ if (phy at order == "postorder") {
+ edge.post <- phy at edge
+ } else {
+ edge.post <- reorder(phy, "postorder")@edge
+ }
+ for (i in seq_along(edge.post[,2])) {
+ if (node.keep[edge.post[i,2]]) {
+ node.keep[edge.post[i,1]] <- TRUE
+ }
+ }
+ } else {
+ node.keep[nTips(phy) + seq_len(nNodes(phy))] <- TRUE
+ }
+ edge.new <- phy at edge[phy at edge[,2] %in% nodes[node.keep], ]
+
+ ## remove singletons
+ edge.length.new <- phy at edge.length
+ edge.label.new <- phy at edge.label
+ singletons <- which(tabulate(na.omit(edge.new[,1]))==1)
+ while (length(singletons)>0) {
+ sing.node <- singletons[1]
+
+ ## update edge matrix
+ edges.drop <- which(edge.new==sing.node, arr.ind=TRUE)[,"row"]
+ sing.edges <- edge.new[edges.drop,]
+ edge.new[edges.drop[2], ] <- c(sing.edges[2,1], sing.edges[1,2])
+ edge.new <- edge.new[-edges.drop[1], ]
+
+ ## update edge lengths and edge labels
+ edge.names.drop <- makeEdgeNames(sing.edges)
+ edge.name.new <- paste(sing.edges[2,1], sing.edges[1,2], sep="-")
+ edge.length.new[edge.name.new] <-
+ sum(edge.length.new[edge.names.drop])
+ edge.length.new <- edge.length.new[-match(edge.names.drop,
+ names(edge.length.new))]
+ edge.label.new[edge.name.new] <- NA
+ edge.label.new <- edge.label.new[-match(edge.names.drop,
+ names(edge.label.new))]
+
+ singletons <- which(tabulate(na.omit(edge.new[,1]))==1)
+ }
+
+ ## remove dropped elements from tip.label and node.label
+ tip.label.new <- phy at tip.label[names(phy at tip.label) %in% edge.new]
+ node.label.new <- phy at node.label[names(phy at node.label) %in% edge.new]
+
+ ## subset and order edge.length and edge.label with respect to edge
+ edge.names <- makeEdgeNames(edge.new)
+ edge.length.new <- edge.length.new[edge.names]
+ edge.label.new <- edge.label.new[edge.names]
+
+ if (!trim.internal) {
+ ## make sure now-terminal internal nodes are treated as tips
+ tip.now <- setdiff(edge.new[,2], edge.new[,1])
+ tip.add <- tip.now[tip.now>nTips(phy)]
+ if (length(tip.add)>0) {
+ ind <- match(tip.add, names(node.label.new))
+
+ ## node renumbering workaround to satisfy plot method
+ newid <- sapply(tip.add, function(x) descendants(phy, x)[1])
+ names(node.label.new)[ind] <- newid
+ edge.new[match(tip.add, edge.new)] <- newid
+ tip.now[match(tip.add, tip.now)] <- newid
+
+ tip.label.new <- c(tip.label.new, node.label.new[ind])
+ node.label.new <- node.label.new[-ind]
+ isTip <- edge.new %in% tip.now
+ edge.new[isTip] <- match(edge.new[isTip],
+ sort(unique.default(edge.new[isTip])))
+ }
+ }
+
+ ## renumber nodes in the edge matrix
+ edge.new[] <- match(edge.new, sort(unique.default(edge.new)))
+
+ ## update corresponding element names in the other slots
+ edge.names <- makeEdgeNames(edge.new)
+ names(edge.length.new) <- edge.names
+ names(edge.label.new) <- edge.names
+ tip.label.new <- tip.label.new[order(as.numeric(names(tip.label.new)))]
+ names(tip.label.new) <- seq_along(tip.label.new)
+ names(node.label.new) <- seq_along(node.label.new) + length(tip.label.new)
+
+ ## create and return new phylo4 object
+ ## NOTE: a faster but looser approach would be to replace the phy
+ ## slots with their new values (including Nnode) and return phy
+ phylo4(x=edge.new, edge.length = edge.length.new, tip.label =
+ tip.label.new, node.label = node.label.new, edge.label =
+ edge.label.new, annote=phy at annote)
+})
+
## trace("prune", browser, signature = "phylo4d")
## untrace("prune", signature = "phylo4d")
setMethod("prune", "phylo4d", function(phy, tip, trim.internal=TRUE,
subtree=FALSE, ...) {
tree <- extractTree(phy)
- phytr <- DropTip(tree, tip, trim.internal, subtree)
+ phytr <- prune(tree, tip, trim.internal, subtree)
## create temporary phylo4 object with unique labels
tmpLbl <- .genlab("n", nTips(phy)+nNodes(phy))
tmpPhy <- tree
labels(tmpPhy, "all") <- tmpLbl
- tmpPhytr <- DropTip(tmpPhy, getNode(phy, tip), trim.internal, subtree)
+ tmpPhytr <- prune(tmpPhy, getNode(phy, tip), trim.internal, subtree)
## get node numbers to keep
oldLbl <- labels(tmpPhy, "all")
Modified: pkg/inst/unitTests/runit.subset.R
===================================================================
--- pkg/inst/unitTests/runit.subset.R 2009-08-31 21:54:27 UTC (rev 610)
+++ pkg/inst/unitTests/runit.subset.R 2009-09-02 07:10:14 UTC (rev 611)
@@ -1,6 +1,9 @@
#
# --- Test subset.R ---
#
+
+# load test comparison objects
+load("trees.RData")
# Create sample tree for testing (ape::phylo object)
tr <- read.tree(text="(((t1:0.2,(t2:0.1,t3:0.1):0.15):0.5,t4:0.7):0.2,t5:1):0.4;")
@@ -21,8 +24,7 @@
# subset 4 tips
checkEquals(tr.sub4, subset(tr, tips.include=c(1, 2, 4, 5)))
checkEquals(tr.sub4, subset(tr, tips.exclude=3))
- checkEquals(tr.sub4, subset(tr, tips.include=c("t1", "t2", "t4",
- "t5")))
+ checkEquals(tr.sub4, subset(tr, tips.include=c("t1", "t2", "t4", "t5")))
checkEquals(tr.sub4, subset(tr, tips.exclude="t3"))
# check variants that should all return the original object
checkEquals(tr, subset(tr))
@@ -34,9 +36,6 @@
}
test.subset.phylo4 <- function() {
- phy <- phylo4(tr)
- phy.sub2 <- phylo4(tr.sub2)
- phy.sub4 <- phylo4(tr.sub4)
# subset 2 tips
checkEquals(phy.sub2, subset(phy, tips.include=c(2, 5)))
checkEquals(phy.sub2, subset(phy, tips.exclude=c(1, 3, 4)))
@@ -63,11 +62,6 @@
}
test.subset.phylo4d <- function() {
- phyd <- phylo4d(tr, data.frame(x=1:5, row.names=paste("t", 1:5, sep="")))
- phyd.sub2 <- phylo4d(tr.sub2, data.frame(x=c(2,5),
- row.names=paste("t", c(2,5), sep="")))
- phyd.sub4 <- phylo4d(tr.sub4, data.frame(x=c(1,2,4,5),
- row.names=paste("t", c(1,2,4,5), sep="")))
# subset 2 tips
checkEquals(phyd.sub2, subset(phyd, tips.include=c(2, 5)))
checkEquals(phyd.sub2, subset(phyd, tips.exclude=c(1, 3, 4)))
Modified: pkg/man/prune-methods.Rd
===================================================================
--- pkg/man/prune-methods.Rd 2009-08-31 21:54:27 UTC (rev 610)
+++ pkg/man/prune-methods.Rd 2009-09-02 07:10:14 UTC (rev 611)
@@ -21,9 +21,7 @@
\item{phy = "phylo"}{drop tips}
}
}
-\note{At the moment, this simply wraps \code{ape::drop.tip}.
- Renamed from \code{drop.tip} to \code{prune} in order to
- avoid conflicts with \code{ape}.}
+\note{The phylo method simply wraps \code{ape::drop.tip}.}
\usage{
\S4method{prune}{phylo4}(phy, tip, trim.internal = TRUE, subtree = FALSE, \dots)
\S4method{prune}{phylo4d}(phy, tip, trim.internal = TRUE, subtree = FALSE, \dots)
@@ -53,53 +51,54 @@
If \code{subtree = TRUE}, the returned tree has one or several
terminal branches indicating how many tips have been removed (with a
label \code{"[x_tips]"}). This is done for as many monophyletic groups
- that have been deleted
+ that have been deleted. This option is only implemented for the
+ phylo method, and thus if used on phylo4 or phylo4d objects, they will
+ be coerced to phylo
Note that \code{subtree = TRUE} implies \code{trim.internal = TRUE}
- To understand how the option \code{root.edge} works, see the examples
- below
+ Handling of root edge lengths differs between \code{phylobase} and
+ \code{ape}. For phylo4/phylo4d objects, the new root edge length
+ will be the total distance back through the original root edge (or NA
+ if the original edge had length NA)
+
+ A \code{root.edge} argument can be passed to the phylo method. To
+ understand how this option works, see the examples below
}
\value{
an object of class \code{"phylo4"}
}
-\author{Emmanuel Paradis \email{Emmanuel.Paradis at mpl.ird.fr}
-(original ape version)}
+\author{
+ Jim Regetz \email{regetz at nceas.ucsb.edu} (phylo4 method),
+ Emmanuel Paradis \email{Emmanuel.Paradis at mpl.ird.fr
+ (original ape version used in phylo method)}
+}
\examples{
-require(ape)
-data(bird.families)
-tip <- c(
-"Eopsaltriidae", "Acanthisittidae", "Pittidae", "Eurylaimidae",
-"Philepittidae", "Tyrannidae", "Thamnophilidae", "Furnariidae",
-"Formicariidae", "Conopophagidae", "Rhinocryptidae", "Climacteridae",
-"Menuridae", "Ptilonorhynchidae", "Maluridae", "Meliphagidae",
-"Pardalotidae", "Petroicidae", "Irenidae", "Orthonychidae",
-"Pomatostomidae", "Laniidae", "Vireonidae", "Corvidae",
-"Callaeatidae", "Picathartidae", "Bombycillidae", "Cinclidae",
-"Muscicapidae", "Sturnidae", "Sittidae", "Certhiidae",
-"Paridae", "Aegithalidae", "Hirundinidae", "Regulidae",
-"Pycnonotidae", "Hypocoliidae", "Cisticolidae", "Zosteropidae",
-"Sylviidae", "Alaudidae", "Nectariniidae", "Melanocharitidae",
-"Paramythiidae","Passeridae", "Fringillidae")
-plot(prune(bird.families, tip),cex=0.5)
-plot(prune(bird.families, tip, trim.internal = FALSE),cex=0.5)
-data(bird.orders)
-plot(prune(bird.orders, 6:23, subtree = TRUE))
-plot(prune(bird.orders, c(1:5, 20:23), subtree = TRUE))
+data(geospiza)
-### Examples of the use of `root.edge'
-tr <- as(ape::read.tree(text = "(A:1,(B:1,(C:1,(D:1,E:1):1):1):1):1;"),"phylo4")
-prune(tr, c("A", "B"), root.edge = 0) # = (C:1,(D:1,E:1):1);
-prune(tr, c("A", "B"), root.edge = 1) # = (C:1,(D:1,E:1):1):1;
-prune(tr, c("A", "B"), root.edge = 2) # = (C:1,(D:1,E:1):1):2;
-prune(tr, c("A", "B"), root.edge = 3) # = (C:1,(D:1,E:1):1):3;
+## Pruning phylo4 objects
+geotree <- extractTree(geospiza)
+tip <- c("difficilis", "fortis", "fuliginosa", "fusca", "olivacea",
+ "pallida", "parvulus", "scandens")
+plot(prune(geotree, tip))
+plot(prune(geotree, tip, trim.internal = FALSE))
+plot(prune(geotree, 1:11, subtree = TRUE))
-## Dropping tips on phylo4d objects
+## Pruning phylo4d objects
r1 <- as(rcoal(5), "phylo4")
td <- data.frame(a=1:5, row.names=paste("t", 1:5, sep=""))
nd <- data.frame(a=6:9, row.names=nodeId(r1, "internal"))
r2 <- phylo4d(r1, tip.data=td, node.data=nd)
r3 <- prune(r2, "t1")
+
+## Pruning phylo objects, including the use of `root.edge'
+require(ape)
+tr <- as(ape::read.tree(text = "(A:1,(B:1,(C:1,(D:1,E:1):1):1):1):1;"),"phylo4")
+prune(tr, c("A", "B"), root.edge = 0) # = (C:1,(D:1,E:1):1);
+prune(tr, c("A", "B"), root.edge = 1) # = (C:1,(D:1,E:1):1):1;
+prune(tr, c("A", "B"), root.edge = 2) # = (C:1,(D:1,E:1):1):2;
+prune(tr, c("A", "B"), root.edge = 3) # = (C:1,(D:1,E:1):1):3;
+
}
\keyword{manip}
\keyword{methods}
Added: pkg/tests/trees.RData
===================================================================
(Binary files differ)
Property changes on: pkg/tests/trees.RData
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
More information about the Phylobase-commits
mailing list