From 1090d5990d4b6f7feb10c87638f4229f53891eb7 Mon Sep 17 00:00:00 2001 From: paradis Date: Sat, 21 Nov 2009 13:56:35 +0000 Subject: [PATCH] last changes for ape 2.4-1 git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@100 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 14 +++++++++++++- DESCRIPTION | 2 +- R/ace.R | 3 ++- R/as.phylo.R | 16 ++++++++-------- R/collapse.singles.R | 4 ++-- R/read.nexus.R | 16 ++++++++++++---- R/root.R | 10 +++++----- R/vcv.phylo.R | 19 ++++++++++++++----- man/read.nexus.Rd | 18 +++++++++++++----- man/vcv.phylo.Rd | 42 +++++++++++++++++++++++++++++++----------- src/tree_build.c | 22 ++++++++-------------- 11 files changed, 109 insertions(+), 57 deletions(-) diff --git a/ChangeLog b/ChangeLog index 8b30895..61c911f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -6,6 +6,11 @@ NEW FEATURES o rtree() and rcoal() now accept a numeric vector for the 'br' argument. + o vcv() is a new generic function with methods for the classes "phylo" + and "corPhyl" so that it is possible to calculate the var-cov matrix + for "transformation models". vcv.phylo() can still be used for trees + of class "phylo"; its argument 'cor' has been renamed 'corr'. + BUG FIXES @@ -14,14 +19,21 @@ BUG FIXES o read.nexus() shuffled tip labels when the trees have no branch lengths and there is a TRANSLATE block. + o read.nexus() does not try to translate node labels if there is a + translation table in the NEXUS file. See ?read.nexus for a + clarification on this behaviour. + o plot.multiPhylo() crashed R when plotting a list of trees with - "compressed tip labels". + compressed tip labels. o write.nexus() did not translate the taxa names when asked for. o plot.phylo(type = "fan") did not rotate the tip labels correctly when the tree has branch lengths. + o ace(type = "continuous", method = "ML") now avoids sigma² being + negative (which resulted in an error). + CHANGES IN APE VERSION 2.4 diff --git a/DESCRIPTION b/DESCRIPTION index 08cc2b1..37ec283 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape Version: 2.4-1 -Date: 2009-11-10 +Date: 2009-11-21 Title: Analyses of Phylogenetics and Evolution Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Gangolf Jobb, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, Damien de Vienne Maintainer: Emmanuel Paradis diff --git a/R/ace.R b/R/ace.R index 91994be..1dd0caa 100644 --- a/R/ace.R +++ b/R/ace.R @@ -1,4 +1,4 @@ -## ace.R (2009-06-10) +## ace.R (2009-11-12) ## Ancestral Character Estimation @@ -59,6 +59,7 @@ did not match: the former were ignored in the analysis.') if (model == "BM") { tip <- phy$edge[, 2] <= nb.tip dev.BM <- function(p) { + if (p[1] < 0) return(1e100) # in case sigma² is negative x1 <- p[-1][phy$edge[, 1] - nb.tip] x2 <- numeric(length(x1)) x2[tip] <- x[phy$edge[tip, 2]] diff --git a/R/as.phylo.R b/R/as.phylo.R index 2280652..e0ee65f 100644 --- a/R/as.phylo.R +++ b/R/as.phylo.R @@ -36,27 +36,27 @@ as.phylo <- function (x, ...) as.phylo.hclust <- function(x, ...) { N <- dim(x$merge)[1] - edge <- matrix(NA, 2*N, 2) + edge <- matrix(0L, 2*N, 2) edge.length <- numeric(2*N) ## `node' gives the number of the node for the i-th row of x$merge - node <- numeric(N) - node[N] <- N + 2 - cur.nod <- N + 3 - j <- 1 + node <- integer(N) + node[N] <- N + 2L + cur.nod <- N + 3L + j <- 1L for (i in N:1) { edge[j:(j + 1), 1] <- node[i] for (l in 1:2) { - k <- j + l - 1 + k <- j + l - 1L if (x$merge[i, l] > 0) { edge[k, 2] <- node[x$merge[i, l]] <- cur.nod - cur.nod <- cur.nod + 1 + cur.nod <- cur.nod + 1L edge.length[k] <- x$height[i] - x$height[x$merge[i, l]] } else { edge[k, 2] <- -x$merge[i, l] edge.length[k] <- x$height[i] } } - j <- j + 2 + j <- j + 2L } if (is.null(x$labels)) x$labels <- as.character(1:(N + 1)) diff --git a/R/collapse.singles.R b/R/collapse.singles.R index 4d477ab..9762bb0 100644 --- a/R/collapse.singles.R +++ b/R/collapse.singles.R @@ -29,7 +29,7 @@ collapse.singles <- function(tree) prev.node <- which(xmat[, 2] == i) next.node <- which(xmat[, 1] == i) xmat[prev.node, 2] <- xmat[next.node, 2] - xmat <- xmat[xmat[, 1] != i, ] ## drop + xmat <- xmat[xmat[, 1] != i, ] # drop ## changed by EP for the new coding of "phylo" (2006-10-05): ## xmat[xmat < i] <- xmat[xmat < i] + 1 ## adjust indices xmat[xmat > i] <- xmat[xmat > i] - 1 ## adjust indices @@ -37,7 +37,7 @@ collapse.singles <- function(tree) elen[prev.node] <- elen[prev.node] + elen[next.node] ## added by Elizabeth Purdom (2008-06-19): if (!is.null(node.lab)) node.lab <- node.lab[-c(i - ntip)] - nnode <- nnode - 1 + nnode <- nnode - 1L ## end elen <- elen[-next.node] } diff --git a/R/read.nexus.R b/R/read.nexus.R index 8c9abc2..969d42c 100644 --- a/R/read.nexus.R +++ b/R/read.nexus.R @@ -1,4 +1,4 @@ -## read.nexus.R (2009-10-27) +## read.nexus.R (2009-11-21) ## Read Tree File in Nexus Format @@ -9,12 +9,20 @@ .treeBuildWithTokens <- function(x) { + ## remove potential node labels; see ?read.nexus for justification + node.label <- gsub("[:;].*$", "", strsplit(x, ")")[[1]][-1]) + has.node.labels <- FALSE + if (any(node.label != "")) { + x <- gsub(")[^:]*:", "):", x) + x <- gsub(")[^:]*;", ");", x) # if there's no root edge + has.node.labels <- TRUE + } phy <- .Call("treeBuildWithTokens", x, PACKAGE = "ape") dim(phy[[1]]) <- c(length(phy[[1]])/2, 2) - nms <- c("edge", "edge.length", "Nnode", "node.label") - if (length(phy) == 5) nms <- c(nms, "root.edge") + nms <- c("edge", "edge.length", "Nnode", "root.edge") + if (length(phy) == 3) nms <- nms[-4] names(phy) <- nms - if (!sum(phy[[4]])) phy[[4]] <- NULL + if (has.node.labels) phy$node.label <- node.label class(phy) <- "phylo" phy } diff --git a/R/root.R b/R/root.R index aadd6c1..d02044f 100644 --- a/R/root.R +++ b/R/root.R @@ -1,4 +1,4 @@ -## root.R (2009-09-09) +## root.R (2009-11-15) ## Root of Phylogenetic Trees @@ -35,11 +35,11 @@ unroot <- function(phy) ## 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 { @@ -50,12 +50,12 @@ unroot <- function(phy) ## 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 diff --git a/R/vcv.phylo.R b/R/vcv.phylo.R index 01cbff7..2473280 100644 --- a/R/vcv.phylo.R +++ b/R/vcv.phylo.R @@ -1,4 +1,4 @@ -## vcv.phylo.R (2009-07-06) +## vcv.phylo.R (2009-11-19) ## Phylogenetic Variance-Covariance or Correlation Matrix @@ -7,10 +7,10 @@ ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. -vcv.phylo <- function(phy, model = "Brownian", cor = FALSE) +vcv <- function(phy, ...) UseMethod("vcv") + +vcv.phylo <- function(phy, model = "Brownian", corr = FALSE, ...) { - if (!inherits(phy, "phylo")) - stop('object "phy" is not of class "phylo"') if (is.null(phy$edge.length)) stop("the tree has no branch lengths") @@ -45,7 +45,7 @@ vcv.phylo <- function(phy, model = "Brownian", cor = FALSE) foo(n + 1, 0, dim(phy$edge)[1]) vcv <- vcv[1:n, 1:n] - if (cor) { + if (corr) { ## This is inspired from the code of `cov2cor' (2005-09-08): M <- vcv Is <- sqrt(1/M[1 + 0:(n - 1)*(n + 1)]) @@ -55,3 +55,12 @@ vcv.phylo <- function(phy, model = "Brownian", cor = FALSE) rownames(vcv) <- colnames(vcv) <- phy$tip.label vcv } + +vcv.corPhyl <- function(phy, corr = FALSE, ...) +{ + labels <- attr(phy, "tree")$tip.label + dummy.df <- data.frame(1:length(labels), row.names = labels) + res <- nlme::corMatrix(Initialize.corPhyl(phy, dummy.df), corr = corr) + dimnames(res) <- list(labels, labels) + res +} diff --git a/man/read.nexus.Rd b/man/read.nexus.Rd index 43a0e34..a517e89 100644 --- a/man/read.nexus.Rd +++ b/man/read.nexus.Rd @@ -8,20 +8,28 @@ read.nexus(file, tree.names = NULL) \item{file}{a file name specified by either a variable of mode character, or a double-quoted string.} \item{tree.names}{if there are several trees to be read, a vector of - mode character that gives names to the individual trees; if - \code{NULL} (the default), the trees are named \code{"tree1"}, - \code{"tree2"}, ...} + mode character that gives names to the individual trees.} } \description{ This function reads one or several trees in a NEXUS file. } \details{ The present implementation tries to follow as much as possible the - NEXUS standard. Only the block ``TREES'' is read; the other data can be - read with other functions (e.g., \code{\link{read.dna}}, + NEXUS standard (but see the restriction below on TRANSLATION + tables). Only the block ``TREES'' is read; the other data can be read + with other functions (e.g., \code{\link{read.dna}}, \code{\link[utils]{read.table}}, ...). A trace of the original data is kept with the attribute \code{"origin"} (see below). + If a TRANSLATION table is present it is assumed that only the tip + labels are translated and they are all translated with integers + without gap. Consequently, if nodes have labels in the tree(s) they + are read as they are and not looked for in the translation table. The + logic behind this is that in the vast majority of cases, node labels + will be support values rather than proper taxa names. This is + consistent with \code{\link{write.nexus}} which translates only the + tip labels. + `read.nexus' tries to represent correctly trees with a badly represented root edge (i.e. with an extra pair of parentheses). For instance, the tree "((A:1,B:1):10);" will be read like "(A:1,B:1):10;" diff --git a/man/vcv.phylo.Rd b/man/vcv.phylo.Rd index 001bc68..7af6975 100644 --- a/man/vcv.phylo.Rd +++ b/man/vcv.phylo.Rd @@ -1,26 +1,31 @@ \name{vcv.phylo} +\alias{vcv} \alias{vcv.phylo} +\alias{vcv.corPhyl} \title{Phylogenetic Variance-covariance or Correlation Matrix} \usage{ -vcv.phylo(phy, model = "Brownian", cor = FALSE) +vcv(phy, ...) +\method{vcv}{phylo}(phy, model = "Brownian", corr = FALSE, ...) +\method{vcv}{corPhyl}(phy, corr = FALSE, ...) } \arguments{ - \item{phy}{an object of class \code{"phylo"}.} + \item{phy}{an object of the correct class (see above).} \item{model}{a character giving the model used to compute the - variances and covariances of the phynotype; by default - \code{"Brownian"}. Currently only the Brownian model is available.} - \item{cor}{a logical indicating whether the correlation matrix should + variances and covariances; only \code{"Brownian"} is available.} + \item{corr}{a logical indicating whether the correlation matrix should be returned (\code{TRUE}); by default the variance-covariance matrix is returned (\code{FALSE}).} + \item{\dots}{further arguments to be passed to or from other methods.} } \description{ This function computes the expected variances and covariances of a - continuous phenotype assuming it evolves under a given model - (currently only the model of Brownian motion is available). + continuous trait assuming it evolves under a given model. + + This is a generic function with methods for objects of class + \code{"phylo"} and \code{"corPhyl"}. } \value{ - a numeric matrix with the names of the tips (as given by the \code{tip.label} - of the argument \code{phy}) as colnames and rownames. + a numeric matrix with the names of the tips as colnames and rownames. } \references{ Garland, T. Jr. and Ives, A. R. (2000) Using the past to predict the @@ -28,9 +33,24 @@ vcv.phylo(phy, model = "Brownian", cor = FALSE) comparative methods. \emph{American Naturalist}, \bold{155}, 346--364. } -\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}} +\author{Emmanuel Paradis} +\note{ + Do not confuse this function with \code{\link[stats]{vcov}} which + computes the variance-covariance matrix among parameters of a fitted + model object. +} \seealso{ - \code{\link{read.tree}} to read tree files in Newick format + \code{\link{corBrownian}}, \code{\link{corMartins}}, + \code{\link{corGrafen}}, \code{\link{corPagel}}, + \code{\link{corBlomberg}} +} +\examples{ +tr <- rtree(5) +## all are the same: +vcv(tr) +vcv(corBrownian(1, tr)) +vcv(corPagel(1, tr)) +vcv(corPagel(0, tr)) } \keyword{manip} \keyword{multivariate} diff --git a/src/tree_build.c b/src/tree_build.c index a166ffc..aeb917c 100644 --- a/src/tree_build.c +++ b/src/tree_build.c @@ -1,6 +1,6 @@ -/* tree_build.c 2008-03-09 */ +/* tree_build.c 2009-11-21 */ -/* Copyright 2008 Emmanuel Paradis */ +/* Copyright 2008-2009 Emmanuel Paradis */ /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -51,16 +51,15 @@ void decode_edge(const char *x, int a, int b, int *node, double *w) l = j - 1;\ while (e[l + nedge] != curnode) l--;\ decode_edge(x, ps + 1, pt - 1, &tmpi, &tmpd);\ - nl[curnode - ntip - 1] = tmpi;\ el[l] = tmpd;\ curnode = e[l] SEXP treeBuildWithTokens(SEXP nwk) { const char *x; - int n, i, ntip = 1, nnode = 0, nedge, *e, *nl, curnode, node, j, *skeleton, nsk = 0, ps, pr, pt, tmpi, l; + int n, i, ntip = 1, nnode = 0, nedge, *e, curnode, node, j, *skeleton, nsk = 0, ps, pr, pt, tmpi, l; double *el, tmpd; - SEXP node_label, edge, edge_length, Nnode, phy; + SEXP edge, edge_length, Nnode, phy; PROTECT(nwk = coerceVector(nwk, STRSXP)); x = CHAR(STRING_ELT(nwk, 0)); @@ -79,14 +78,11 @@ SEXP treeBuildWithTokens(SEXP nwk) } nedge = ntip + nnode - 1; - PROTECT(node_label = allocVector(INTSXP, nnode)); PROTECT(Nnode = allocVector(INTSXP, 1)); PROTECT(edge = allocVector(INTSXP, nedge*2)); PROTECT(edge_length = allocVector(REALSXP, nedge)); INTEGER(Nnode)[0] = nnode; - nl = INTEGER(node_label); - memset(nl, 0, nnode*sizeof(int)); e = INTEGER(edge); el = REAL(edge_length); @@ -130,22 +126,20 @@ SEXP treeBuildWithTokens(SEXP nwk) /* is there a root edge? */ if (ps < n - 2) { - PROTECT(phy = allocVector(VECSXP, 5)); + PROTECT(phy = allocVector(VECSXP, 4)); SEXP root_edge; decode_edge(x, ps + 1, n - 2, &tmpi, &tmpd); PROTECT(root_edge = allocVector(REALSXP, 1)); - nl[0] = tmpi; REAL(root_edge)[0] = tmpd; - SET_VECTOR_ELT(phy, 4, root_edge); + SET_VECTOR_ELT(phy, 3, root_edge); UNPROTECT(1); - } else PROTECT(phy = allocVector(VECSXP, 4)); + } else PROTECT(phy = allocVector(VECSXP, 3)); SET_VECTOR_ELT(phy, 0, edge); SET_VECTOR_ELT(phy, 1, edge_length); SET_VECTOR_ELT(phy, 2, Nnode); - SET_VECTOR_ELT(phy, 3, node_label); - UNPROTECT(6); + UNPROTECT(5); return phy; } -- 2.39.2