X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Froot.R;h=8526b8f3d4f74add15240b88eca14a8267785496;hb=cc266bd204f991845994007b428ba9e20f92963c;hp=d02044fa627360bafaedf5055e8bfe3dce3d7d8c;hpb=1090d5990d4b6f7feb10c87638f4229f53891eb7;p=ape.git diff --git a/R/root.R b/R/root.R index d02044f..8526b8f 100644 --- a/R/root.R +++ b/R/root.R @@ -1,8 +1,8 @@ -## root.R (2009-11-15) +## 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. @@ -10,19 +10,19 @@ 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"') + 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 @@ -61,13 +61,18 @@ unroot <- function(phy) 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") @@ -86,34 +91,21 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE) ## 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) @@ -247,7 +239,7 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE) 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] @@ -283,7 +275,7 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE) } } ## 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