From: paradis Date: Wed, 19 Dec 2012 07:35:43 +0000 (+0000) Subject: a few changes.... X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=a2fc961dffe9f9b7994ed880e68c03b2334dc341 a few changes.... git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@201 6e262413-ae40-0410-9e79-b911bd7a66b7 --- diff --git a/DESCRIPTION b/DESCRIPTION index a81d208..c2302fd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape Version: 3.0-7 -Date: 2012-11-28 +Date: 2012-12-19 Title: Analyses of Phylogenetics and Evolution Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Andrei-Alin Popescu, Klaus Schliep, Korbinian Strimmer, Damien de Vienne Maintainer: Emmanuel Paradis diff --git a/NEWS b/NEWS index 9d57433..f6acbff 100644 --- a/NEWS +++ b/NEWS @@ -6,13 +6,28 @@ NEW FEATURES o The new function 'where' searches patterns in DNA sequences. o pic() gains an option 'rescaled.tree = FALSE' to return the tree - with its branch lengths rescaled for the PICs calculation. + with its branch lengths rescaled for the PIC calculation. + + o clustal(), muscle(), and tcoffee() gain an option + 'original.ordering = TRUE' to ease the comparisons of + alignments. + + o plot.phylo() has a new option, open.angle, used when plotting + circular trees. BUG FIXES o drop.tip() shuffled node labels on some trees. + o axisPhylo() now works correctly with circular trees, and gives a + sensible error message when type = "r" or "u". + + +OTHER CHANGES + + o .compressTipLabel() is 10 times faster thanks to Joseph Brown. + CHANGES IN APE VERSION 3.0-6 diff --git a/R/DNA.R b/R/DNA.R index 1d0ce79..31ef4bb 100644 --- a/R/DNA.R +++ b/R/DNA.R @@ -481,7 +481,7 @@ where <- function(x, pattern) s <- as.integer(length(x)) if (s < p) stop("sequence shorter than the pattern") ans <- .C("where", x, pat, s, p, integer(s), 0L, - DUP = FALSE, NAOK = TRUE, PACKAGE = "apex") + DUP = FALSE, NAOK = TRUE, PACKAGE = "ape") n <- ans[[6]] if (n) ans[[5]][seq_len(n)] - p + 2L else integer() } diff --git a/R/clustal.R b/R/clustal.R index 656ecb2..c777945 100644 --- a/R/clustal.R +++ b/R/clustal.R @@ -1,4 +1,4 @@ -## clustal.R (2012-01-12) +## clustal.R (2012-11-28) ## Multiple Sequence Alignment with External Applications @@ -9,7 +9,7 @@ clustal <- function(x, pw.gapopen = 10, pw.gapext = 0.1, gapopen = 10, gapext = 0.2, exec = NULL, - MoreArgs = "", quiet = TRUE) + MoreArgs = "", quiet = TRUE, original.ordering = TRUE) { os <- Sys.info()[1] if (is.null(exec)) { @@ -32,10 +32,12 @@ clustal <- function(x, pw.gapopen = 10, pw.gapext = 0.1, opts <- paste(prefix, suffix, sep = "=", collapse = " ") opts <- paste(opts, MoreArgs) system(paste(exec, opts), ignore.stdout = quiet) - read.dna(outf, "clustal") + res <- read.dna(outf, "clustal") + if (original.ordering) res <- res[labels(x), ] + res } -muscle <- function(x, exec = "muscle", MoreArgs = "", quiet = TRUE) +muscle <- function(x, exec = "muscle", MoreArgs = "", quiet = TRUE, original.ordering = TRUE) { if (missing(x)) { system(exec) @@ -50,10 +52,12 @@ muscle <- function(x, exec = "muscle", MoreArgs = "", quiet = TRUE) if (quiet) opts <- paste(opts, "-quiet") opts <- paste(opts, MoreArgs) system(paste(exec, opts)) - read.dna(outf, "fasta") + res <- read.dna(outf, "fasta") + if (original.ordering) res <- res[labels(x), ] + res } -tcoffee <- function(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE) +tcoffee <- function(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE, original.ordering = TRUE) { if (missing(x)) { system(exec) @@ -68,5 +72,7 @@ tcoffee <- function(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE) opts <- paste(inf, MoreArgs) if (quiet) opts <- paste(opts, "-quiet=nothing") system(paste(exec, opts)) - read.dna("input_tcoffee.aln", "clustal") + res <- read.dna("input_tcoffee.aln", "clustal") + if (original.ordering) res <- res[labels(x), ] + res } diff --git a/R/cophyloplot.R b/R/cophyloplot.R index 85f091f..29e5b3c 100644 --- a/R/cophyloplot.R +++ b/R/cophyloplot.R @@ -42,11 +42,11 @@ cophyloplot <- } plotCophylo2 <- - function (x, y, assoc = assoc, use.edge.length = use.edge.length, - space = space, length.line = length.line, gap = gap, - type = type, return = return, col = col, lwd=lwd, lty=lty, - show.tip.label = show.tip.label, - font = font, ...) + function(x, y, assoc = assoc, use.edge.length = use.edge.length, + space = space, length.line = length.line, gap = gap, + type = type, return = return, col = col, lwd=lwd, lty=lty, + show.tip.label = show.tip.label, + font = font, ...) { res <- list() ###choice of the minimum space between the trees diff --git a/R/dist.topo.R b/R/dist.topo.R index 5777482..230f009 100644 --- a/R/dist.topo.R +++ b/R/dist.topo.R @@ -1,4 +1,4 @@ -## dist.topo.R (2012-03-13) +## dist.topo.R (2012-12-12) ## Topological Distances, Tree Bipartitions, ## Consensus Trees, and Bootstrapping Phylogenies @@ -68,32 +68,28 @@ dist.topo <- function(x, y, method = "PH85") ## 'x' is a list of objects of class "phylo" possibly with no class if (!is.null(attr(x, "TipLabel"))) return(x) ref <- x[[1]]$tip.label - if (any(table(ref) != 1)) - stop("some tip labels are duplicated in tree no. 1") n <- length(ref) - Ntree <- length(x) - if (Ntree > 1) { - for (i in 2:Ntree) { - label <- x[[i]]$tip.label - if (!identical(label, ref)) { - if (length(label) != length(ref)) - stop(paste("tree no.", i, "has a different number of tips")) - ilab <- match(label, ref) - ## can use tabulate here because 'ilab' contains integers - if (any(is.na(ilab))) - stop(paste("tree no.", i, "has different tip labels")) -### the test below does not seem useful anymore -### if (any(tabulate(ilab) > 1)) -### stop(paste("some tip labels are duplicated in tree no.", i)) -### - ie <- match(1:n, x[[i]]$edge[, 2]) - x[[i]]$edge[ie, 2] <- ilab - } - x[[i]]$tip.label <- NULL + if (length(unique(ref)) != n) + stop("some tip labels are duplicated in tree no. 1") + + ## serious improvement by Joseph W. Brown! + relabel <- function (y) { + label <- y$tip.label + if (!identical(label, ref)) { + if (length(label) != length(ref)) + stop(paste("tree ", y, "has a different number of tips")) + ilab <- match(label, ref) + if (any(is.na(ilab))) + stop(paste("tree ", y, "has different tip labels")) + ie <- match(1:n, y$edge[, 2]) + y$edge[ie, 2] <- ilab } + y$tip.label <- NULL + y } - x[[1]]$tip.label <- NULL + x <- lapply(x, relabel) attr(x, "TipLabel") <- ref + class(x) <- "multiPhylo" x } diff --git a/R/drop.tip.R b/R/drop.tip.R index f1b413f..113c448 100644 --- a/R/drop.tip.R +++ b/R/drop.tip.R @@ -1,4 +1,4 @@ -## drop.tip.R (2012-11-22) +## drop.tip.R (2012-11-29) ## Remove Tips in a Phylogenetic Tree @@ -204,8 +204,8 @@ drop.tip <- if (is.null(phy$node.label)) rep("NA", length(node2tip)) else phy$node.label[node2tip - Ntip] } - if (!is.null(phy$node.label)) - phy$node.label <- phy$node.label[-(node2tip - Ntip)] +# if (!is.null(phy$node.label)) +# phy$node.label <- phy$node.label[-(node2tip - Ntip)] phy$tip.label <- c(phy$tip.label, new.tip.label) } diff --git a/R/plot.phylo.R b/R/plot.phylo.R index 4c23c6f..96d9aa4 100644 --- a/R/plot.phylo.R +++ b/R/plot.phylo.R @@ -1,4 +1,4 @@ -## plot.phylo.R (2012-10-20) +## plot.phylo.R (2012-12-19) ## Plot Phylogenies @@ -15,7 +15,8 @@ plot.phylo <- adj = NULL, srt = 0, no.margin = FALSE, root.edge = FALSE, label.offset = 0, underscore = FALSE, x.lim = NULL, y.lim = NULL, direction = "rightwards", lab4ut = "horizontal", - tip.color = "black", plot = TRUE, rotate.tree = 0, ...) + tip.color = "black", plot = TRUE, rotate.tree = 0, + open.angle = 0, ...) { Ntip <- length(x$tip.label) if (Ntip < 2) { @@ -122,13 +123,15 @@ plot.phylo <- xx <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, z$edge.length) } } else { - rotate.tree <- 2 * pi * rotate.tree/360 + twopi <- 2 * pi + rotate.tree <- twopi * rotate.tree/360 switch(type, "fan" = { ## if the tips are not in the same order in tip.label ## and in edge[, 2], we must reorder the angles: we ## use `xx' to store temporarily the angles TIPS <- x$edge[which(x$edge[, 2] <= Ntip), 2] - xx <- seq(0, 2*pi*(1 - 1/Ntip), 2*pi/Ntip) + xx <- seq(0, twopi * (1 - 1/Ntip) - twopi * open.angle/360, + length.out = Ntip) theta <- double(Ntip) theta[TIPS] <- xx theta <- c(theta, numeric(Nnode)) @@ -140,8 +143,8 @@ plot.phylo <- r <- 1/r } theta <- theta + rotate.tree - xx <- r*cos(theta) - yy <- r*sin(theta) + xx <- r * cos(theta) + yy <- r * sin(theta) }, "unrooted" = { nb.sp <- .nodeDepth(Ntip, Nnode, z$edge, Nedge) XY <- if (use.edge.length) @@ -157,7 +160,7 @@ plot.phylo <- ## radius: X <- 1 - X/Ntip ## angle (1st compute the angles for the tips): - yy <- c((1:Ntip)*2*pi/Ntip, rep(0, Nnode)) + yy <- c((1:Ntip)*twopi/Ntip, rep(0, Nnode)) Y <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy) xx <- X * cos(Y + rotate.tree) yy <- X * sin(Y + rotate.tree) @@ -375,7 +378,6 @@ if (plot) { if (type %in% c("fan", "radial")) { xx.tips <- xx[1:Ntip] yy.tips <- yy[1:Ntip] - ## using atan2 considerably facilitates things compared to acos... angle <- atan2(yy.tips, xx.tips) # in radians if (label.offset) { xx.tips <- xx.tips + label.offset * cos(angle) diff --git a/R/scales.R b/R/scales.R index 7176d1f..c2e9d61 100644 --- a/R/scales.R +++ b/R/scales.R @@ -1,4 +1,4 @@ -## scales.R (2011-05-31) +## scales.R (2012-12-19) ## Add a Scale Bar or Axis to a Phylogeny Plot @@ -76,7 +76,14 @@ add.scale.bar <- function(x, y, length = NULL, ask = FALSE, axisPhylo <- function(side = 1, ...) { lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv) - if (lastPP$type %in% c("phylogram", "cladogram")) { + type <- lastPP$type + + if (type == "unrooted") + stop("axisPhylo() not available for unrooted plots; try add.scale.bar()") + if (type == "radial") + stop("axisPhylo() not meaningful for this type of plot") + + if (type %in% c("phylogram", "cladogram")) { if (lastPP$direction %in% c("rightwards", "leftwards")) { x <- pretty(lastPP$xx) if (lastPP$direction == "rightwards") maxi <- max(lastPP$xx) @@ -92,6 +99,33 @@ axisPhylo <- function(side = 1, ...) x <- -x } } + axis(side = side, at = c(maxi - x), labels = abs(x), ...) + } else { # type == "fan" + n <- lastPP$Ntip + xx <- lastPP$xx[1:n]; yy <- lastPP$yy[1:n] + r0 <- max(sqrt(xx^2 + yy^2)) + firstandlast <- c(1, n) + theta0 <- mean(atan2(yy[firstandlast], xx[firstandlast])) + x0 <- r0 * cos(theta0); y0 <- r0 * sin(theta0) + inc <- diff(pretty(c(0, r0))[1:2]) + srt <- 360*theta0/(2*pi) + coef <- -1 + if (abs(srt) > 90) { + srt <- srt + 180 + coef <- 1 + } + len <- 0.025 * r0 # the length of tick marks + r <- r0 + while (r > 1e-8) { + x <- r * cos(theta0); y <- r * sin(theta0) + if (len/r < 1) { + ra <- sqrt(len^2 + r^2); thetaa <- theta0 + coef * asin(len/r) + xa <- ra * cos(thetaa); ya <- ra * sin(thetaa) + segments(xa, ya, x, y) + text(xa, ya, r0 - r, srt = srt, adj = c(0.5, 1.1), ...) + } + r <- r - inc + } + segments(x, y, x0, y0) } - axis(side = side, at = c(maxi - x), labels = abs(x), ...) } diff --git a/man/clustal.Rd b/man/clustal.Rd index 0e00198..afb8d73 100644 --- a/man/clustal.Rd +++ b/man/clustal.Rd @@ -10,9 +10,11 @@ \usage{ clustal(x, pw.gapopen = 10, pw.gapext = 0.1, gapopen = 10, gapext = 0.2, exec = NULL, - MoreArgs = "", quiet = TRUE) -muscle(x, exec = "muscle", MoreArgs = "", quiet = TRUE) -tcoffee(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE) + MoreArgs = "", quiet = TRUE, original.ordering = TRUE) +muscle(x, exec = "muscle", MoreArgs = "", quiet = TRUE, + original.ordering = TRUE) +tcoffee(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE, + original.ordering = TRUE) } \arguments{ \item{x}{an object of class \code{"DNAbin"}.} @@ -25,6 +27,8 @@ tcoffee(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE) \item{MoreArgs}{a character string giving additional options.} \item{quiet}{a logical: the default is to not print on \R's console the messages from the external program.} + \item{original.ordering}{a logical specifying whether to return the + aligned sequences in the same order than in \code{x}.} } \details{ \code{clustal} tries to guess the name of the executable program diff --git a/man/plot.phylo.Rd b/man/plot.phylo.Rd index ee89287..6d08946 100644 --- a/man/plot.phylo.Rd +++ b/man/plot.phylo.Rd @@ -10,7 +10,7 @@ root.edge = FALSE, label.offset = 0, underscore = FALSE, x.lim = NULL, y.lim = NULL, direction = "rightwards", lab4ut = "horizontal", tip.color = "black", plot = TRUE, - rotate.tree = 0, ...) + rotate.tree = 0, open.angle = 0, ...) \method{plot}{multiPhylo}(x, layout = 1, ...) } \arguments{ @@ -100,6 +100,9 @@ \item{layout}{the number of trees to be plotted simultaneously.} \item{\dots}{further arguments to be passed to \code{plot} or to \code{plot.phylo}.} + \item{open.angle}{if \code{type = "f"}, the angle in degrees left + blank. Use a non-zero value if you want to call + \code{\link{axisPhylo}} after the tree is plotted.} } \description{ These functions plot phylogenetic trees on the current graphical diff --git a/man/where.Rd b/man/where.Rd index faddeba..09e0ab8 100644 --- a/man/where.Rd +++ b/man/where.Rd @@ -1,4 +1,4 @@ -\names{where} +\name{where} \alias{where} \title{Find Patterns in DNA Sequences} \description{