From 767a9ed6bc4444aac3dc1a26a91fc3986e99665c Mon Sep 17 00:00:00 2001 From: paradis Date: Tue, 7 Jul 2009 06:35:45 +0000 Subject: [PATCH] several bug fixes while in JKT git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@81 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 20 ++++++++++++++++++++ DESCRIPTION | 4 ++-- R/all.equal.phylo.R | 34 +++++++++++++++++++++------------- R/chronopl.R | 39 ++++++++++++++++++--------------------- R/dist.topo.R | 16 +++++++++++----- R/drop.tip.R | 4 ++-- R/root.R | 25 ++++++++++++++----------- R/vcv.phylo.R | 3 ++- Thanks | 4 ++-- 9 files changed, 92 insertions(+), 57 deletions(-) diff --git a/ChangeLog b/ChangeLog index 98602f9..dbe6fe0 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,23 @@ + CHANGES IN APE VERSION 2.3-2 + + +BUG FIXES + + o all.equal.phylo() did not compare unrooted trees correctly. + + o dist.topo(... method = "PH85") did not treat unrooted trees + correctly (thanks to Tim Wallstrom for the fix). + + o root() sometimes failed to test for the monophyly of the + outgroup correctly. + + o extract.clade() sometimes included too many edges. + + o vcv.phylo() did not work correctly when the tree is in + "pruningwise" order. + + + CHANGES IN APE VERSION 2.3-1 diff --git a/DESCRIPTION b/DESCRIPTION index 4a0986b..0ff3acf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape -Version: 2.3-1 -Date: 2009-06-23 +Version: 2.3-2 +Date: 2009-07-06 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/all.equal.phylo.R b/R/all.equal.phylo.R index cd05437..79533d5 100644 --- a/R/all.equal.phylo.R +++ b/R/all.equal.phylo.R @@ -1,4 +1,4 @@ -## all.equal.phylo.R (2006-09-12) +## all.equal.phylo.R (2009-07-05) ## ## Global Comparison of two Phylogenies @@ -8,18 +8,18 @@ ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. -### Recherche de la correspondance entre deux arbres -### Parcours en profondeur et en parallèle des deux arbres (current et target) -### current, target: les deux arbres à comparer -### use.edge.length: faut-il comparer les longueurs de branches ? -### use.tip.label: faut-il comparer les étiquettes de feuilles ou seulement la -### topologie des deux arbres ? -### index.return: si TRUE, retourner la matrice de correspondance entre noeuds -### et feuilles, une matrice à deux colonnes (current et target) avec pour -### chaque ligne des paires d'identifiants de noeuds/feuilles, tels qu'ils -### apparaissent dans l'attribut 'edge' des objets phylo -### tolerance, scale: paramètres de comparaison des longueurs de branches -### (voir 'all.equal') +## Recherche de la correspondance entre deux arbres +## Parcours en profondeur et en parallèle des deux arbres (current et target) +## current, target: les deux arbres à comparer +## use.edge.length: faut-il comparer les longueurs de branches ? +## use.tip.label: faut-il comparer les étiquettes de feuilles ou seulement la +## topologie des deux arbres ? +## index.return: si TRUE, retourner la matrice de correspondance entre noeuds +## et feuilles, une matrice à deux colonnes (current et target) avec pour +## chaque ligne des paires d'identifiants de noeuds/feuilles, tels qu'ils +## apparaissent dans l'attribut 'edge' des objets phylo +## tolerance, scale: paramètres de comparaison des longueurs de branches +## (voir 'all.equal') all.equal.phylo <- function(target, current, use.edge.length = TRUE, use.tip.label = TRUE, @@ -66,6 +66,14 @@ all.equal.phylo <- function(target, current, root1 <- Ntip1 + 1 root2 <- Ntip2 + 1 if (root1 != root2) return(FALSE) + ## Fix by EP so that unrooted trees are correctly compared: + if (!is.rooted(target) && !is.rooted(current)) { + outg <- target$tip.label[1] + if (! outg %in% current$tip.label) return(FALSE) + target <- root(target, outg) + current <- root(current, outg) + } + ## End result <- same.node(root1, root2) if (!isTRUE(index.return)) return(!is.null(result)) if (is.null(result)) return(result) diff --git a/R/chronopl.R b/R/chronopl.R index 6a0bd89..570fa26 100644 --- a/R/chronopl.R +++ b/R/chronopl.R @@ -1,15 +1,16 @@ -## chronopl.R (2008-11-04) +## chronopl.R (2009-07-06) ## Molecular Dating With Penalized Likelihood -## Copyright 2005-2008 Emmanuel Paradis +## Copyright 2005-2009 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. -chronopl <- function(phy, lambda, age.min = 1, age.max = NULL, - node = "root", S = 1, tol = 1e-8, - CV = FALSE, eval.max = 500, iter.max = 500, ...) +chronopl <- + function(phy, lambda, age.min = 1, age.max = NULL, + node = "root", S = 1, tol = 1e-8, + CV = FALSE, eval.max = 500, iter.max = 500, ...) { n <- length(phy$tip.label) ROOT <- n + 1 @@ -84,27 +85,23 @@ chronopl <- function(phy, lambda, age.min = 1, age.max = NULL, real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]] while (any(real.edge.length <= 0)) { for (i in EDGES) { - if (real.edge.length[i] <= 0) { - if (e[i, 1] %in% node) { - ini.time[e[i, 2]] <- - ini.time[e[, 2]] - 2*real.edge.length[i] - next - } - if (e[i, 2] %in% node) { - ini.time[e[i, 1]] <- - ini.time[e[, 1]] + 2*real.edge.length[i] - next - } - ini.time[e[i, 2]] <- - ini.time[e[, 2]] - real.edge.length[i] - ini.time[e[i, 1]] <- - ini.time[e[, 1]] + real.edge.length[i] + if (real.edge.length[i] > 0) next + if (e[i, 1] %in% node) { + ini.time[e[i, 2]] <- ini.time[e[1, 2]] - 2 * real.edge.length[i] + next } + if (e[i, 2] %in% node) { + ini.time[e[i, 1]] <- ini.time[e[1, 1]] + 2 * real.edge.length[i] + next + } + browser() + ini.time[e[i, 2]] <- ini.time[e[1, 2]] - real.edge.length[i] + ini.time[e[i, 1]] <- ini.time[e[1, 1]] + real.edge.length[i] } real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]] + print(min(real.edge.length)) } } - ## `unknown.ages' will contain the index of the nodes of unknown age: unknown.ages <- n + 1:m diff --git a/R/dist.topo.R b/R/dist.topo.R index 6741dcc..1248c49 100644 --- a/R/dist.topo.R +++ b/R/dist.topo.R @@ -1,4 +1,4 @@ -## dist.topo.R (2009-05-10) +## dist.topo.R (2009-07-06) ## Topological Distances, Tree Bipartitions, ## Consensus Trees, and Bootstrapping Phylogenies @@ -15,21 +15,27 @@ dist.topo <- function(x, y, method = "PH85") n <- length(x$tip.label) bp1 <- .Call("bipartition", x$edge, n, x$Nnode, PACKAGE = "ape") bp1 <- lapply(bp1, function(xx) sort(x$tip.label[xx])) - bp2 <- .Call("bipartition", y$edge, n, y$Nnode, PACKAGE = "ape") - bp2 <- lapply(bp2, function(xx) sort(y$tip.label[xx])) + ## fix by Tim Wallstrom: + bp2.tmp <- .Call("bipartition", y$edge, n, y$Nnode, PACKAGE = "ape") + bp2 <- lapply(bp2.tmp, function(xx) sort(y$tip.label[xx])) + bp2.comp <- lapply(bp2.tmp, function(xx) setdiff(1:n, xx)) + bp2.comp <- lapply(bp2.comp, function(xx) sort(y$tip.label[xx])) + ## End q1 <- length(bp1) q2 <- length(bp2) if (method == "PH85") { p <- 0 for (i in 1:q1) { for (j in 1:q2) { - if (identical(all.equal(bp1[[i]], bp2[[j]]), TRUE)) { + if (identical(bp1[[i]], bp2[[j]]) | + identical(bp1[[i]], bp2.comp[[j]])) { p <- p + 1 break } } } - dT <- if (q1 == q2) 2*(q1 - p) else 2*(min(q1, q2) - p) + abs(q1 - q2) + dT <- q1 + q2 - 2 * p # same than: + ##dT <- if (q1 == q2) 2*(q1 - p) else 2*(min(q1, q2) - p) + abs(q1 - q2) } if (method == "BHV01") { dT <- 0 diff --git a/R/drop.tip.R b/R/drop.tip.R index 6e497cb..acdbe68 100644 --- a/R/drop.tip.R +++ b/R/drop.tip.R @@ -1,4 +1,4 @@ -## drop.tip.R (2009-06-23) +## drop.tip.R (2009-07-06) ## Remove Tips in a Phylogenetic Tree @@ -29,7 +29,7 @@ extract.clade <- function(phy, node, root.edge = 0) root.node <- which(phy$edge[, 2] == node) start <- root.node + 1 # start of the clade looked for anc <- phy$edge[root.node, 1] # the ancestor of 'node' - next.anc <- which(phy$edge[-(1:start), 1] == anc) # find the next occurence of 'anc' + next.anc <- which(phy$edge[-(1:start), 1] <= anc) # find the next occurence of 'anc' or an 'older' node keep <- if (length(next.anc)) start + 0:(next.anc[1] - 1) else start:Nedge diff --git a/R/root.R b/R/root.R index 832ab0f..3abdf65 100644 --- a/R/root.R +++ b/R/root.R @@ -1,4 +1,4 @@ -## root.R (2009-05-10) +## root.R (2009-07-06) ## Root of Phylogenetic Trees @@ -11,11 +11,10 @@ is.rooted <- function(phy) { if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"') - if (!is.null(phy$root.edge)) return(TRUE) + if (!is.null(phy$root.edge)) TRUE else if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2) - return(FALSE) - else return(TRUE) + FALSE else TRUE } unroot <- function(phy) @@ -66,8 +65,7 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE) { if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"') - ord <- attr(phy, "order") - if (!is.null(ord) && ord == "pruningwise") phy <- reorder(phy) + phy <- reorder(phy) n <- length(phy$tip.label) ROOT <- n + 1 if (!is.null(node)) { @@ -88,7 +86,6 @@ 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) { - msg <- "the specified outgroup is not monophyletic" seq.nod <- .Call("seq_root2tip", phy$edge, n, phy$Nnode, PACKAGE = "ape") sn <- seq.nod[outgroup] @@ -108,9 +105,15 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE) ## (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] } N <- Nedge(phy) @@ -128,7 +131,7 @@ root <- function(phy, outgroup, node = NULL, resolve.root = FALSE) if (!is.null(phy$edge.length)) phy$edge.length <- c(phy$edge.length[a], 0, phy$edge.length[b]) - phy$Nnode <- phy$Nnode + 1 + phy$Nnode <- phy$Nnode + 1L ## node renumbering (see comments below) newNb <- integer(n + oldNnode) newNb[newroot] <- n + 1L diff --git a/R/vcv.phylo.R b/R/vcv.phylo.R index a51c6e4..01cbff7 100644 --- a/R/vcv.phylo.R +++ b/R/vcv.phylo.R @@ -1,4 +1,4 @@ -## vcv.phylo.R (2009-05-10) +## vcv.phylo.R (2009-07-06) ## Phylogenetic Variance-Covariance or Correlation Matrix @@ -37,6 +37,7 @@ vcv.phylo <- function(phy, model = "Brownian", cor = FALSE) } } + phy <- reorder(phy) n <- length(phy$tip.label) n.node <- phy$Nnode N <- n.node + n diff --git a/Thanks b/Thanks index 157d4cc..abbeca8 100644 --- a/Thanks +++ b/Thanks @@ -6,8 +6,8 @@ comments, or bug reports: thanks to all of you! Significant bug fixes were provided by Cécile Ané, James Bullard, Éric Durand, Olivier François, Rich FitzJohn Bret Larget, Nick Matzke, -Elizabeth Purdom, Dan Rabosky, Klaus Schliep, Li-San Wang, Yan Wong, -and Peter Wragg. Contact me if I forgot someone. +Elizabeth Purdom, Dan Rabosky, Klaus Schliep, Tim Wallstrom, Li-San +Wang, Yan Wong, and Peter Wragg. Contact me if I forgot someone. Kurt Hornik, of the R Core Team, helped in several occasions to fix some problems and bugs. -- 2.39.2