-## root.R (2008-02-12)
+## root.R (2010-02-11)
## Root of Phylogenetic Trees
-## Copyright 2004-2008 Emmanuel Paradis
+## Copyright 2004-2010 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
is.rooted <- function(phy)
{
- if (!"phylo" %in% class(phy))
- stop('object "phy" is not of class "phylo"')
- if (!is.null(phy$root.edge)) return(TRUE)
+ if (!inherits(phy, "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)
- return(FALSE)
- else return(TRUE)
+ if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2)
+ FALSE else TRUE
}
unroot <- function(phy)
{
- if (class(phy) != "phylo")
- stop('object "phy" is not of class "phylo"')
+ 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("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
## nodes by adding 1, except the root (this remains the
## origin of the tree).
nb.tip <- length(phy$tip.label)
- ROOT <- nb.tip + 1
+ ROOT <- nb.tip + 1L
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) {
+ if (phy$edge[EDGEROOT[1], 2] == ROOT + 1L) {
j <- EDGEROOT[2]
i <- EDGEROOT[1]
} else {
## cladewise order.
phy$edge <- phy$edge[-i, ]
nodes <- phy$edge > ROOT # renumber all nodes except the root
- phy$edge[nodes] <- phy$edge[nodes] - 1
+ phy$edge[nodes] <- phy$edge[nodes] - 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$Nnode <- phy$Nnode - 1
+ phy$Nnode <- phy$Nnode - 1L
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 (class(phy) != "phylo")
- stop('object "phy" is not of class "phylo"')
- ord <- attr(phy, "order")
- if (!is.null(ord) && ord == "pruningwise") phy <- reorder(phy)
+ if (!inherits(phy, "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) {
- msg <- "the specified outgroup is not monophyletic"
seq.nod <- .Call("seq_root2tip", phy$edge, n,
phy$Nnode, PACKAGE = "ape")
sn <- seq.nod[outgroup]
## (below is slightly faster than calling "bipartition")
desc <- which(unlist(lapply(seq.nod,
function(x) any(x %in% newroot))))
- if (length(outgroup) != length(desc)) stop(msg)
- ## both vectors below are already sorted:
- if (!all(outgroup == desc)) stop(msg)
+ msg <- "the specified outgroup is not monophyletic"
+ 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)
+ }
} else newroot <- phy$edge[which(phy$edge[, 2] == outgroup), 1]
}
- if (newroot == ROOT) return(phy)
+ N <- Nedge(phy)
+ oldNnode <- phy$Nnode
+ if (newroot == ROOT) {
+ if (resolve.root) {
+ snw <- which(phy$edge[, 1] == newroot)
+ if (length(snw) > 2) {
+ a <- snw[1]:(snw[2] - 1)
+ b <- snw[2]:N
+ newnod <- oldNnode + n + 1
+ phy$edge[snw[-1], 1] <- newnod
+ 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$Nnode <- phy$Nnode + 1L
+ ## node renumbering (see comments below)
+ newNb <- integer(n + oldNnode)
+ newNb[newroot] <- n + 1L
+ sndcol <- phy$edge[, 2] > n
+ phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
+ (n + 2):(n + phy$Nnode)
+ phy$edge[, 1] <- newNb[phy$edge[, 1]]
+ }
+ }
+ return(phy)
+ }
phy$root.edge <- NULL # just in case...
Nclade <- tabulate(phy$edge[, 1])[ROOT] # degree of the root node
## if only 2 edges connect to the root, we have to fuse them:
fuseRoot <- Nclade == 2
- N <- Nedge(phy)
start <- which(phy$edge[, 1] == ROOT)
end <- c(start[-1] - 1, N)
o <- integer(N)
ne <- length(s)
o[NEXT:(NEXT + ne - 1L)] <- s
- oldNnode <- phy$Nnode
if (fuseRoot) {
phy$Nnode <- oldNnode - 1
N <- N - 1L
phy$edge[, 1] <- newNb[phy$edge[, 1]]
if (!is.null(phy$node.label)) {
+ #browser()
newNb <- newNb[-(1:n)]
if (fuseRoot) {
newNb <- newNb[-1]
phy$node.label <- phy$node.label[-1]
}
phy$node.label <- phy$node.label[order(newNb)]
- if (resolve.root)
- phy$node.label <- c(phy$node.label[1], NA, phy$node.label[-1])
+ if (resolve.root) {
+ phy$node.label[is.na(phy$node.label)] <- phy$node.label[1]
+ phy$node.label[1] <- NA
+ ##phy$node.label <- c(phy$node.label[1], NA, phy$node.label[-1])
+ ##phy$node.label <- c("NA", phy$node.label)
+ }
}
phy
}