-## root.R (2009-09-09)
+## root.R (2011-08-05)
## Root of Phylogenetic Trees
-## Copyright 2004-2009 Emmanuel Paradis
+## Copyright 2004-2011 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
is.rooted <- function(phy)
{
if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
+ stop('object "phy" is not of class "phylo"')
if (!is.null(phy$root.edge)) TRUE
else
- if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2)
- FALSE else TRUE
+ if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2)
+ FALSE else TRUE
}
unroot <- function(phy)
{
if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
- if (dim(phy$edge)[1] < 3)
- stop("cannot unroot a tree with two edges.")
+ stop('object "phy" is not of class "phylo"')
+ N <- dim(phy$edge)[1]
+ if (N < 3)
+ stop("cannot unroot a tree with less than three edges.")
+
## delete FIRST the root.edge (in case this is sufficient to
## unroot the tree, i.e. there is a multichotomy at the root)
if (!is.null(phy$root.edge)) phy$root.edge <- NULL
if (!is.rooted(phy)) return(phy)
- ## We remove one of the edges coming from the root, and
- ## eventually adding the branch length to the other one
- ## also coming from the root.
- ## In all cases, the node deleted is the 2nd one (numbered
- ## nb.tip+2 in `edge'), so we simply need to renumber the
- ## nodes by adding 1, except the root (this remains the
- ## origin of the tree).
- nb.tip <- length(phy$tip.label)
- ROOT <- nb.tip + 1
- EDGEROOT <- which(phy$edge[, 1] == ROOT)
- ## j: the target where to stick the edge
- ## i: the edge to delete
- if (phy$edge[EDGEROOT[1], 2] == ROOT + 1) {
- j <- EDGEROOT[2]
- i <- EDGEROOT[1]
+
+ n <- length(phy$tip.label)
+ ROOT <- n + 1L
+
+### EDGEROOT[1]: the edge to delete
+### EDGEROOT[2]: the target where to stick the edge to delete
+
+### If the tree is in pruningwise order, then the last two edges
+### are those connected to the root; the node situated in
+### phy$edge[N - 2L, 1L] will be the new root...
+
+ ophy <- attr(phy, "order")
+ if (!is.null(ophy) && ophy == "pruningwise") {
+ NEWROOT <- phy$edge[N - 2L, 1L]
+ EDGEROOT <- c(N, N - 1L)
} else {
- j <- EDGEROOT[1]
- i <- EDGEROOT[2]
+
+### ... otherwise, we remove one of the edges coming from
+### the root, and eventually adding the branch length to
+### the other one also coming from the root.
+### In all cases, the node deleted is the 2nd one (numbered
+### nb.tip+2 in 'edge'), so we simply need to renumber the
+### nodes by adding 1, except the root (this remains the
+### origin of the tree).
+
+ EDGEROOT <- which(phy$edge[, 1L] == ROOT)
+ NEWROOT <- ROOT + 1L
}
- ## This should work whether the tree is in pruningwise or
- ## cladewise order.
- phy$edge <- phy$edge[-i, ]
- nodes <- phy$edge > ROOT # renumber all nodes except the root
- phy$edge[nodes] <- phy$edge[nodes] - 1
+
+ ## make sure EDGEROOT is ordered as described above:
+ if (phy$edge[EDGEROOT[1L], 2L] != NEWROOT)
+ EDGEROOT <- EDGEROOT[2:1]
+
+ phy$edge <- phy$edge[-EDGEROOT[1L], ]
+
+ s <- phy$edge == NEWROOT # renumber the new root
+ phy$edge[s] <- ROOT
+
+ s <- phy$edge > NEWROOT # renumber all nodes greater than the new root
+ phy$edge[s] <- phy$edge[s] - 1L
+
if (!is.null(phy$edge.length)) {
- phy$edge.length[j] <- phy$edge.length[j] + phy$edge.length[i]
- phy$edge.length <- phy$edge.length[-i]
+ phy$edge.length[EDGEROOT[2L]] <-
+ phy$edge.length[EDGEROOT[2L]] + phy$edge.length[EDGEROOT[1L]]
+ phy$edge.length <- phy$edge.length[-EDGEROOT[1L]]
+ }
+
+ phy$Nnode <- phy$Nnode - 1L
+
+ if (!is.null(phy$node.label)) {
+ if (NEWROOT == n + 2L)
+ phy$node.label <- phy$node.label[-1]
+ else {
+ lbs <- phy$node.label
+ tmp <- lbs[NEWROOT - n]
+ lbs <- lbs[-c(1, NEWROOT)]
+ phy$node.label <- c(tmp, lbs)
+ }
}
- phy$Nnode <- phy$Nnode - 1
- if (!is.null(phy$node.label))
- phy$node.label <- phy$node.label[-2]
phy
}
-root <- function(phy, outgroup, node = NULL, resolve.root = FALSE)
+root <- function(phy, outgroup, node = NULL,
+ resolve.root = FALSE, interactive = FALSE)
{
if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
+ stop('object "phy" is not of class "phylo"')
phy <- reorder(phy)
n <- length(phy$tip.label)
- ROOT <- n + 1
+ ROOT <- n + 1L
+ if (interactive) {
+ node <- identify(phy)$nodes
+ cat("You have set resolve.root =", resolve.root, "\n")
+ }
if (!is.null(node)) {
if (node <= n)
stop("incorrect node#: should be greater than the number of taxa")
## First check that the outgroup is monophyletic--
## unless there's only one tip specified of course
if (length(outgroup) > 1) {
- seq.nod <- .Call("seq_root2tip", phy$edge, n,
- phy$Nnode, PACKAGE = "ape")
- sn <- seq.nod[outgroup]
- ## We go from the root to the tips: the sequence of nodes
- ## is identical until the MRCA:
- newroot <- ROOT
- i <- 2 # we start at the 2nd position since the root
- # of the tree is a common ancestor to all tips
- repeat {
- x <- unique(unlist(lapply(sn, "[", i)))
- if (length(x) != 1) break
- newroot <- x
- i <- i + 1
- }
- ## Check that all descendants of this node
- ## are included in the outgroup.
- ## (below is slightly faster than calling "bipartition")
- desc <- which(unlist(lapply(seq.nod,
- function(x) any(x %in% newroot))))
- msg <- "the specified outgroup is not monophyletic"
+ pp <- prop.part(phy)
ingroup <- (1:n)[-outgroup]
- ## 'outgroup' and 'desc' are already sorted:
- if (newroot != ROOT) {
- if (!identical(outgroup, desc) && !identical(ingroup, desc))
- stop(msg)
- } else { # otherwise check monophyly of the ingroup
- if (!is.monophyletic(phy, ingroup)) stop(msg)
+ newroot <- 0L
+ for (i in 2:phy$Nnode) {
+ if (identical(pp[[i]], ingroup)) {
+ newroot <- i + n
+ break
+ }
+ if (identical(pp[[i]], outgroup)) {
+ newroot <- phy$edge[which(phy$edge[, 2] == i + n), 1]
+ break
+ }
}
+ if (!newroot)
+ stop("the specified outgroup is not monophyletic")
} else newroot <- phy$edge[which(phy$edge[, 2] == outgroup), 1]
}
N <- Nedge(phy)
phy$edge <- rbind(phy$edge[a, ], c(ROOT, newnod),
phy$edge[b, ])
if (!is.null(phy$edge.length))
- phy$edge.length <-
- c(phy$edge.length[a], 0, phy$edge.length[b])
+ phy$edge.length <-
+ c(phy$edge.length[a], 0, phy$edge.length[b])
phy$Nnode <- phy$Nnode + 1L
## node renumbering (see comments below)
newNb <- integer(n + oldNnode)
o[NEXT:(NEXT + ne - 1L)] <- s
if (fuseRoot) {
- phy$Nnode <- oldNnode - 1
+ phy$Nnode <- oldNnode - 1L
N <- N - 1L
}
phy$edge[INV, ] <- phy$edge[INV, 2:1]
}
}
## N <- N + 1L ... not needed
- phy$Nnode <- phy$Nnode + 1
+ phy$Nnode <- phy$Nnode + 1L
}
## The block below renumbers the nodes so that they conform
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]] <-
- (n + 2):(n + phy$Nnode)
+ phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <- n + 2:phy$Nnode
phy$edge[, 1] <- newNb[phy$edge[, 1]]
if (!is.null(phy$node.label)) {
projects / ape.git / blobdiff