+ CHANGES IN APE VERSION 2.1-2
+
+
+NEW FEATURES
+
+ o root() gains the options 'node' and 'resolve.root'
+ (FALSE by default) as well as its code being improved.
+
+
+
CHANGES IN APE VERSION 2.1-1
Package: ape
-Version: 2.1-1
-Date: 2008-02-01
+Version: 2.1-2
+Date: 2008-02-11
Title: Analyses of Phylogenetics and Evolution
Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong,
Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel,
-## identify.phylo.R (2008-01-14)
+## identify.phylo.R (2008-02-08)
## Graphical Identification of Nodes and Tips
{
cat("Click close to a node of the tree...\n")
xy <- locator(1)
- Ntip <- .last_plot.phylo$Ntip
- d <- sqrt((xy$x - .last_plot.phylo$xx)^2 +
- (xy$y - .last_plot.phylo$yy)^2)
+ Ntip <- get("last_plot.phylo$Ntip", envir = .PlotPhyloEnv)
+ d <- sqrt((xy$x - get("last_plot.phylo$xx", envir = .PlotPhyloEnv))^2 +
+ (xy$y - get("last_plot.phylo$yy", envir = .PlotPhyloEnv))^2)
NODE <- which.min(d)
res <- list()
if (NODE <= Ntip) {
-## nodelabels.R (2007-03-05)
+## nodelabels.R (2008-02-08)
## Labelling Trees
-## Copyright 2004-2007 Emmanuel Paradis, 2006 Ben Bolker, and 2006 Jim Lemon
+## Copyright 2004-2008 Emmanuel Paradis, 2006 Ben Bolker, and 2006 Jim Lemon
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,
col = "black", bg = "lightblue", ...)
{
+ xx <- get("last_plot.phylo$xx", envir = .PlotPhyloEnv)
+ yy <- get("last_plot.phylo$yy", envir = .PlotPhyloEnv)
if (missing(node))
- node <- (.last_plot.phylo$Ntip + 1):length(.last_plot.phylo$xx)
- XX <- .last_plot.phylo$xx[node]
- YY <- .last_plot.phylo$yy[node]
+ node <- (get("last_plot.phylo$Ntip",
+ envir = .PlotPhyloEnv) + 1):length(xx)
+ XX <- xx[node]
+ YY <- yy[node]
BOTHlabels(text, node, XX, YY, adj, frame, pch, thermo,
pie, piecol, col, bg, ...)
}
pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,
col = "black", bg = "yellow", ...)
{
- if (missing(tip)) tip <- 1:.last_plot.phylo$Ntip
- XX <- .last_plot.phylo$xx[tip]
- YY <- .last_plot.phylo$yy[tip]
+ if (missing(tip))
+ tip <- 1:get("last_plot.phylo$Ntip", envir = .PlotPhyloEnv)
+ XX <- get("last_plot.phylo$xx", envir = .PlotPhyloEnv)[tip]
+ YY <- get("last_plot.phylo$yy", envir = .PlotPhyloEnv)[tip]
BOTHlabels(text, tip, XX, YY, adj, frame, pch, thermo,
pie, piecol, col, bg, ...)
}
pch = NULL, thermo = NULL, pie = NULL, piecol = NULL,
col = "black", bg = "lightgreen", ...)
{
+ xx <- get("last_plot.phylo$xx", envir = .PlotPhyloEnv)
+ yy <- get("last_plot.phylo$yy", envir = .PlotPhyloEnv)
+ lastEdge <- get("last_plot.phylo$edge", envir = .PlotPhyloEnv)
if (missing(edge)) {
- sel <- 1:dim(.last_plot.phylo$edge)[1]
- subedge <- .last_plot.phylo$edge
+ sel <- 1:dim(lastEdge)[1]
+ subedge <- lastEdge
} else {
sel <- edge
- subedge <- .last_plot.phylo$edge[sel, , drop = FALSE]
+ subedge <- lastEdge[sel, , drop = FALSE]
}
- if (.last_plot.phylo$type == "phylogram") {
- if(.last_plot.phylo$direction %in% c("rightwards", "leftwards")) {
- XX <- (.last_plot.phylo$xx[subedge[, 1]] +
- .last_plot.phylo$xx[subedge[, 2]]) / 2
- YY <- .last_plot.phylo$yy[subedge[, 2]]
+ if (get("last_plot.phylo$type", envir = .PlotPhyloEnv) == "phylogram") {
+ if(get("last_plot.phylo$direction", envir = .PlotPhyloEnv)
+ %in% c("rightwards", "leftwards")) {
+ XX <- (xx[subedge[, 1]] + xx[subedge[, 2]]) / 2
+ YY <- yy[subedge[, 2]]
} else {
- XX <- .last_plot.phylo$xx[subedge[, 2]]
- YY <- (.last_plot.phylo$yy[subedge[, 1]] +
- .last_plot.phylo$yy[subedge[, 2]]) / 2
+ XX <- xx[subedge[, 2]]
+ YY <- (yy[subedge[, 1]] + yy[subedge[, 2]]) / 2
}
} else {
- XX <- (.last_plot.phylo$xx[subedge[, 1]] +
- .last_plot.phylo$xx[subedge[, 2]]) / 2
- YY <- (.last_plot.phylo$yy[subedge[, 1]] +
- .last_plot.phylo$yy[subedge[, 2]]) / 2
+ XX <- (xx[subedge[, 1]] + xx[subedge[, 2]]) / 2
+ YY <- (yy[subedge[, 1]] + yy[subedge[, 2]]) / 2
}
BOTHlabels(text, sel, XX, YY, adj, frame, pch, thermo,
pie, piecol, col, bg, ...)
-## plot.phylo.R (2008-01-12)
+## plot.phylo.R (2008-02-08)
## Plot Phylogenies
label.offset = label.offset, x.lim = x.lim, y.lim = y.lim,
direction = direction, tip.color = tip.color,
Ntip = Ntip, Nnode = Nnode)
- .last_plot.phylo <<- c(L, list(edge = xe, xx = xx, yy = yy))
+ assing("last_plot.phylo", c(L, list(edge = xe, xx = xx, yy = yy)),
+ envir = .PlotPhyloEnv)
invisible(L)
}
-## root.R (2007-12-21)
+## root.R (2008-02-12)
## Root of Phylogenetic Trees
-## Copyright 2004-2007 Emmanuel Paradis
+## Copyright 2004-2008 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
phy
}
-root <- function(phy, outgroup)
+root <- function(phy, outgroup, node = NULL, resolve.root = FALSE)
{
if (class(phy) != "phylo")
stop('object "phy" is not of class "phylo"')
- if (is.character(outgroup))
- outgroup <- which(phy$tip.label %in% outgroup)
- nb.tip <- length(phy$tip.label)
- if (length(outgroup) == nb.tip) return(phy)
+ ord <- attr(phy, "order")
+ if (!is.null(ord) && ord == "pruningwise") phy <- reorder(phy)
+ n <- length(phy$tip.label)
+ ROOT <- n + 1
+ if (!is.null(node)) {
+ if (node <= n)
+ stop("incorrect node#: should be greater than the number of taxa")
+ outgroup <- NULL
+ newroot <- node
+ } else {
+ if (is.numeric(outgroup)) {
+ if (any(outgroup > n))
+ stop("incorrect taxa#: should not be greater than the number of taxa")
+ outgroup <- sort(outgroup) # used below
+ }
+ if (is.character(outgroup))
+ outgroup <- which(phy$tip.label %in% outgroup)
+ if (length(outgroup) == n) return(phy)
- ## First check that the outgroup is monophyletic--
- ## unless there's only one tip specified of course
- ROOT <- nb.tip + 1
- if (length(outgroup) > 1) {
- msg <- "the specified outgroup is not monophyletic"
- ## If all tips in `tip' are not contiguous, then
- ## no need to go further:
- if (!all(diff(outgroup) == 1)) stop(msg)
- seq.nod <- .Call("seq_root2tip", phy$edge, nb.tip,
- 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
+ ## 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]
+ ## 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))))
+ if (length(outgroup) != length(desc)) stop(msg)
+ ## both vectors below are already sorted:
+ if (!all(outgroup == desc)) stop(msg)
+ } else newroot <- phy$edge[which(phy$edge[, 2] == outgroup), 1]
+ }
+ if (newroot == ROOT) 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)
+ INV <- logical(N)
+
+ w <- which(phy$edge[, 2] == newroot)
+ nod <- phy$edge[w, 1]
+ i <- w
+ NEXT <- 1L
+
+ ## The loop below starts from the new root and goes up in the edge
+ ## matrix reversing the edges that need to be as well as well
+ ## inverting their order. The other edges must not be changed, so
+ ## their indices are stored in `stack'.
+ ## We then bind the other edges in a straightforward way.
+
+ if (nod != ROOT) {
+ ## it is important that the 3 next lines
+ ## are inside this "if" statement
+ o[NEXT] <- w
+ NEXT <- NEXT + 1L
+ INV[w] <- TRUE
+ i <- w - 1L
+ stack <- 0L
repeat {
- x <- unique(unlist(lapply(sn, "[", i)))
- if (length(x) != 1) break
- newroot <- x
- i <- i + 1
+ if (phy$edge[i, 2] == nod) {
+ if (stack) {
+ o[NEXT:(NEXT + stack - 1L)] <- (i + 1L):(i + stack)
+ NEXT <- NEXT + stack
+ stack <- 0L
+ }
+ if (phy$edge[i, 1] == ROOT) break
+ o[NEXT] <- i
+ NEXT <- NEXT + 1L
+ INV[i] <- TRUE
+ nod <- phy$edge[i, 1]
+ } else stack <- stack + 1L
+ i <- i - 1L
}
- ## Check that all descendants of this node
- ## are included in the outgroup
- ## (1st solution... there may be something smarter)
- desc <- which(unlist(lapply(seq.nod,
- function(x) any(x %in% newroot))))
- if (length(outgroup) != length(desc)) stop(msg)
- if (!all(sort(outgroup) == sort(desc))) stop(msg)
+ }
- } else newroot <- phy$edge[which(phy$edge[, 2] == outgroup), 1]
+ ## we keep the edge leading to the old root if needed:
+ if (!fuseRoot) {
+ o[NEXT] <- i
+ INV[i] <- TRUE
+ NEXT <- NEXT + 1L
+ }
- if (newroot == ROOT) return(phy)
+ endOfOutgroup <- which(phy$edge[(w+1):N, 1] < newroot)[1] + w - 1
+ if (is.na(endOfOutgroup)) endOfOutgroup <- N
+ endOfClade <- end[end >= endOfOutgroup][1]
-### <FIXME>
-### The remaining part of the code has not been improved; this
-### does not seem obvious. This is delayed... (2006-09-23)
-### </FIXME>
-
- ## Invert all branches from the new root to the old one
- i <- which(phy$edge[, 2] == newroot)
- nod <- phy$edge[i, 1]
- while (nod != ROOT) {
- j <- which(phy$edge[, 2] == nod)
- phy$edge[i, 1] <- phy$edge[i, 2]
- phy$edge[i, 2] <- nod
- i <- j
- nod <- phy$edge[i, 1]
+ ## bind the other clades...
+ for (j in 1:Nclade) {
+ if (end[j] == endOfClade) next
+ ## do we have to fuse the two basal edges?
+ if (fuseRoot) {
+ phy$edge[start[j], 1] <- phy$edge[i, 2]
+ if (!is.null(phy$edge.length))
+ phy$edge.length[start[j]] <- phy$edge.length[start[j]] +
+ phy$edge.length[i]
+ } #else {
+ # o[NEXT] <- i#start[j]
+ # NEXT <- NEXT + 1L
+ #}
+ s <- start[j]:end[j]
+ ne <- length(s)
+ o[NEXT:(NEXT + ne - 1L)] <- s
+ NEXT <- NEXT + ne
}
- i.oroot <- which(phy$edge[, 1] == ROOT)
- ## Unroot the tree if there's a basal dichotomy...
- if (length(i.oroot) == 2) {
- j <- i.oroot[which(i.oroot != i)]
- phy$edge[j, 1] <- phy$edge[i, 2]
- phy$edge <- phy$edge[-i, ]
- 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]
+ ## possibly bind the edges below the outgroup till the end of the clade
+ if (all(endOfOutgroup != end)) {
+ j <- (endOfOutgroup + 1L):endOfClade
+ ## we must take care that the branch inversions done above
+ ## may have changed the hierarchy of clades here, so we
+ ## travel from the bottom of this series of edges
+ stack <- 0L
+ inverted <- phy$edge[INV, 1] # <- fails if ', 2]' is used
+ for (k in rev(j)) {
+ if (any(phy$edge[k, 1] == inverted)) {
+ o[NEXT] <- k
+ NEXT <- NEXT + 1L
+ if (stack){
+ o[NEXT:(NEXT + stack - 1L)] <- (k + 1L):(k + stack)
+ NEXT <- NEXT + stack
+ stack <- 0L
+ }
+ } else stack <- stack + 1L
}
- phy$edge[which(phy$edge == newroot)] <- ROOT
- } else {
- ## ... otherwise just invert the root with the newroot
- phy$edge[which(phy$edge == newroot)] <- ROOT
- phy$edge[i.oroot] <- newroot
- ## ... and invert finally! (fixed 2005-11-07)
- phy$edge[i, ] <- rev(phy$edge[i, ])
}
- if (!is.null(phy$node.label))
- ## It's important to not delete the label of the newroot
- ## to keep the positions of the other nodes
- phy$node.label[1] <- phy$node.label[newroot - nb.tip]
- ## Not needed: phy$Nnode <- phy$Nnode - 1
- read.tree(text = write.tree(phy))
+
+ ## ... and the outgroup
+ s <- (w + 1L):endOfOutgroup
+ ne <- length(s)
+ o[NEXT:(NEXT + ne - 1L)] <- s
+
+ oldNnode <- phy$Nnode
+ if (fuseRoot) {
+ phy$Nnode <- oldNnode - 1
+ N <- N - 1L
+ }
+ phy$edge[INV, ] <- phy$edge[INV, 2:1]
+ phy$edge <- phy$edge[o, ]
+ if (!is.null(phy$edge.length))
+ phy$edge.length <- phy$edge.length[o]
+
+ if (resolve.root) {
+ newnod <- oldNnode + n + 1
+ if (length(outgroup) == 1L) {
+ wh <- which(phy$edge[, 2] == outgroup)
+ phy$edge[1] <- newnod
+ phy$edge <-
+ rbind(c(newroot, newnod), phy$edge[-wh, ], phy$edge[wh, ])
+ snw <- which(phy$edge[, 1] == newroot)
+ phy$edge[snw[length(snw) - 1], 1] <- newnod
+ if (!is.null(phy$edge.length)) {
+ phy$edge.length <-
+ c(0, phy$edge.length[-wh], phy$edge.length[wh])
+ }
+ } else {
+ wh <- which(phy$edge[, 1] == newroot)
+ phy$edge[wh[-1], 1] <- newnod
+ s1 <- 1:(wh[2] - 1)
+ s2 <- wh[2]:N
+ phy$edge <-
+ rbind(phy$edge[s1, ], c(newroot, newnod), phy$edge[s2, ])
+ if (!is.null(phy$edge.length)) {
+ tmp <- phy$edge.length[1]
+ phy$edge.length[1] <- 0
+ phy$edge.length <-
+ c(phy$edge.length[s1], tmp, phy$edge.length[s2])
+ }
+ }
+ ## N <- N + 1L ... not needed
+ phy$Nnode <- phy$Nnode + 1
+ }
+
+ ## The block below renumbers the nodes so that they conform
+ ## to the "phylo" format
+ newNb <- integer(n + oldNnode)
+ 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[, 1] <- newNb[phy$edge[, 1]]
+
+ if (!is.null(phy$node.label)) {
+ 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])
+ }
+ phy
}
-## scales.R (2004-12-18)
+## scales.R (2008-02-08)
## Add a Scale Bar or Axis to a Phylogeny Plot
## add.scale.bar: add a scale bar to a phylogeny plot
## axisPhylo: add a scale axis on the side of a phylogeny plot
-## Copyright 2002-2004 Emmanuel Paradis
+## Copyright 2002-2008 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
add.scale.bar <- function(x = 0, y = 1, length = NULL, ...)
{
if (is.null(length)) {
- nb.digit <- ceiling(log10(mean(.last_plot.phylo$xx))) - 2
+ nb.digit <- ceiling(log10(mean(get("last_plot.phylo$xx",
+ envir = .PlotPhyloEnv)))) - 2
length <- eval(parse(text = paste("1e", nb.digit, sep = "")))
}
segments(x, y, x + length, y)
axisPhylo <- function(side = 1, ...)
{
- if (.last_plot.phylo$type %in% c("phylogram", "cladogram")) {
- if (.last_plot.phylo$direction %in% c("rightwards", "leftwards")) {
- x <- pretty(.last_plot.phylo$xx)
- if (.last_plot.phylo$direction == "rightwards")
- maxi <- max(.last_plot.phylo$xx)
+ type <- get("last_plot.phylo$type", envir = .PlotPhyloEnv)
+ direction <- get("last_plot.phylo$direction", envir = .PlotPhyloEnv)
+ if (type %in% c("phylogram", "cladogram")) {
+ if (direction %in% c("rightwards", "leftwards")) {
+ xx <- get("last_plot.phylo$xx", envir = .PlotPhyloEnv)
+ x <- pretty(xx)
+ if (direction == "rightwards") maxi <- max(xx)
else {
- maxi <- min(.last_plot.phylo$xx)
+ maxi <- min(xx)
x <- -x
}
} else {
- x <- pretty(.last_plot.phylo$yy)
- if (.last_plot.phylo$direction == "upwards")
- maxi <- max(.last_plot.phylo$yy)
+ yy <- get("last_plot.phylo$yy", envir = .PlotPhyloEnv)
+ x <- pretty(yy)
+ if (direction == "upwards") maxi <- max(yy)
else {
- maxi <- min(.last_plot.phylo$yy)
+ maxi <- min(yy)
x <- -x
}
}
-## zzz.R (2008-01-14)
+## zzz.R (2008-02-08)
## Library Loading
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
+.PlotPhyloEnv <- new.env()
+
.First.lib <- function(lib, pkg) {
require(nlme, quietly = TRUE)
library.dynam("ape", pkg, lib)
\alias{is.rooted}
\title{Roots Phylogenetic Trees}
\usage{
-root(phy, outgroup)
+root(phy, outgroup, node = NULL, resolve.root = FALSE)
unroot(phy)
is.rooted(phy)
}
\item{phy}{an object of class \code{"phylo"}.}
\item{outgroup}{a vector of mode numeric or character specifying the
new outgroup.}
+ \item{node}{alternatively, a node number where to root the tree.}
+ \item{resolve.root}{a logical specifying whether to resolve the new
+ root as a bifurcating node.}
}
\description{
\code{root} reroots a phylogenetic tree with respect to the specified
- outgroup.
+ outgroup or at the node specified in \code{node}.
\code{unroot} unroots a phylogenetic tree, or returns it unchanged if
it is already unrooted.
one (see examples). If \code{outgroup} is not monophyletic, the
operation fails and an error message is issued.
+ If \code{resolve.root = TRUE}, \code{root} adds a zero-length branch
+ below the MRCA of the ingroup.
+
A tree is considered rooted if either only two branches connect to the
root, or if there is a \code{root.edge} element. In all other cases,
\code{is.rooted} returns \code{FALSE}.
}
\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
\seealso{
- \code{\link{bind.tree}}, \code{\link{drop.tip}}
+ \code{\link{bind.tree}}, \code{\link{drop.tip}},
+ \code{\link{nodelabels}}, \code{\link{identify.phylo}}
}
\examples{
data(bird.orders)
### This is because the tree has been unrooted first before rerooting.
### You can delete the outgroup...
is.rooted(drop.tip(tr, "Struthioniformes"))
-### ... or resolve the basal trichotomy:
+### ... or resolve the basal trichotomy in two ways:
is.rooted(multi2di(tr))
+is.rooted(root(bird.orders, 1, r = TRUE))
### To keep the basal trichotomy but forcing the tree as rooted:
tr$root.edge <- 0
is.rooted(tr)
projects / ape.git / commitdiff