From 7af63afb401f7cc82816b0af19e21c43fab49674 Mon Sep 17 00:00:00 2001 From: paradis Date: Mon, 26 Sep 2011 08:34:50 +0000 Subject: [PATCH] new LTT plot functions + updated reference in yule.time.Rd git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@170 6e262413-ae40-0410-9e79-b911bd7a66b7 --- DESCRIPTION | 4 +- COPYING => LICENCE | 0 NEWS | 15 ++++++ R/ltt.plot.R | 123 +++++++++++++++++++++++++++++++++------------ man/ltt.plot.Rd | 103 +++++++++++++++++++++++++++---------- man/yule.time.Rd | 5 ++ 6 files changed, 188 insertions(+), 62 deletions(-) rename COPYING => LICENCE (100%) diff --git a/DESCRIPTION b/DESCRIPTION index 5501dec..6b51928 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape -Version: 2.7-3 -Date: 2011-08-08 +Version: 2.8 +Date: 2011-09-24 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, Klaus Schliep, Korbinian Strimmer, Damien de Vienne Maintainer: Emmanuel Paradis diff --git a/COPYING b/LICENCE similarity index 100% rename from COPYING rename to LICENCE diff --git a/NEWS b/NEWS index 528e408..776fd4a 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,18 @@ + CHANGES IN APE VERSION 2.8 + + +NEW FEATURES + + o ltt.plot() has two new options: backward and tol. It can now + handle non-ultrametic trees and its internal coding has been + improved. The coordinates of the plot can now be computed with + the new function ltt.plot.coords. + + o The new function ltt.coplot plots on the same scales a tree and + its LTT plot. + + + CHANGES IN APE VERSION 2.7-3 diff --git a/R/ltt.plot.R b/R/ltt.plot.R index a7eae04..17a1d31 100644 --- a/R/ltt.plot.R +++ b/R/ltt.plot.R @@ -1,45 +1,93 @@ -## ltt.plot.R (2009-05-10) +## ltt.plot.R (2011-09-24) ## Lineages Through Time Plot -## Copyright 2002-2009 Emmanuel Paradis +## Copyright 2002-2011 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. -ltt.plot <- function(phy, xlab = "Time", ylab = "N", ...) +ltt.plot.coords <- function(phy, backward = TRUE, tol = 1e-6) { - if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"') - if (!is.binary.tree(phy)) phy <- multi2di(phy) - time <- sort(branching.times(phy), decreasing = TRUE) - N <- 1:(length(time) + 1) - plot.default(-c(time, 0), N, xlab = xlab, ylab = ylab, - xaxs = "r", yaxs = "r", type = "S", ...) + if (is.ultrametric(phy, tol)) { + if (is.binary.tree(phy)) { + N <- numeric(phy$Nnode + 1) + N[] <- 1 + } else { + node.order <- tabulate(phy$edge[, 1]) + N <- node.order[-(1:length(phy$tip.label))] - 1 + } + bt <- branching.times(phy) + names(bt) <- NULL + o <- order(bt, decreasing = TRUE) + time <- c(-bt[o], 0) + if (!is.binary.tree(phy)) N <- c(1, N[o]) + } else { + if (!is.binary.tree(phy)) phy <- multi2di(phy) + n <- Ntip(phy) + m <- phy$Nnode + ROOT <- n + 1L + event <- time.event <- numeric(n + m) + + time.event[ROOT] <- 0 + phy <- reorder(phy) + + for (i in 1:nrow(phy$edge)) + time.event[phy$edge[i, 2]] <- time.event[phy$edge[i, 1]] + phy$edge.length[i] + + present <- max(time.event) + event[1:n] <- -1 + event[ROOT:(n + m)] <- 1 + + ## delete the events that are too close to present: + past.event <- present - time.event > tol + event <- event[past.event] + time.event <- time.event[past.event] + + ## reorder wrt time: + o <- order(time.event) + time.event <- time.event[o] + event <- event[o] + + time <- c(time.event - present, 0) + N <- c(1, event) + } + N <- cumsum(N) + if (!is.null(phy$root.edge)) { + time <- c(time[1] - phy$root.edge, time) + N <- c(1, N) + } + if (!backward) time <- time - time[1] + cbind(time, N) } -ltt.lines <- function(phy, ...) +ltt.plot <- function(phy, xlab = "Time", ylab = "N", + backward = TRUE, tol = 1e-6, ...) { - if (!is.binary.tree(phy)) phy <- multi2di(phy) - time <- sort(branching.times(phy), decreasing = TRUE) - N <- 1:(length(time) + 1) - lines(-c(time, 0), N, type = "S", ...) + if (!inherits(phy, "phylo")) + stop("object \"phy\" is not of class \"phylo\"") + + xy <- ltt.plot.coords(phy, backward, tol) + + plot.default(xy, xlab = xlab, ylab = ylab, xaxs = "r", + yaxs = "r", type = "S", ...) } -mltt.plot <- function(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE, - xlab = "Time", ylab = "N", log = "") +ltt.lines <- function(phy, backward = TRUE, tol = 1e-6, ...) +{ + xy <- ltt.plot.coords(phy, backward, tol) + lines(xy, type = "S", ...) +} + +mltt.plot <- + function(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE, + xlab = "Time", ylab = "N", log = "", backward = TRUE, tol = 1e-6) { - ltt.xy <- function(phy) { - if (!is.binary.tree(phy)) phy <- multi2di(phy) - x <- -c(sort(branching.times(phy), decreasing = TRUE), 0) - names(x) <- NULL - y <- 1:length(x) - cbind(x, y) - } if (inherits(phy, "phylo")) { # if a tree of class "phylo" - TREES <- list(ltt.xy(phy)) + TREES <- list(ltt.plot.coords(phy, backward, tol)) names(TREES) <- deparse(substitute(phy)) } else { # a list of trees - TREES <- lapply(phy, ltt.xy) + TREES <- lapply(phy, ltt.plot.coords, backward = backward, tol = tol) names(TREES) <- names(phy) if (is.null(names(TREES))) names(TREES) <- @@ -52,14 +100,13 @@ mltt.plot <- function(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE, nms <- mc[1:n] for (i in 1:n) { if (class(dts[[i]]) == "phylo") { - a <- list(ltt.xy(dts[[i]])) + a <- list(ltt.plot.coords(dts[[i]], backward, tol)) names(a) <- nms[i] } else { # a list of trees - a <- lapply(dts[[i]], ltt.xy) + a <- lapply(dts[[i]], ltt.plot.coords, backward = backward, tol = tol) names(a) <- names(dts[[i]]) if (is.null(names(a))) - names(a) <- - paste(deparse(substitute(phy)), "-", 1:length(a)) + names(a) <- paste(deparse(substitute(phy)), "-", 1:length(a)) } TREES <- c(TREES, a) } @@ -68,16 +115,26 @@ mltt.plot <- function(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE, xl <- c(min(unlist(lapply(TREES, function(x) min(x[, 1])))), 0) yl <- c(1, max(unlist(lapply(TREES, function(x) max(x[, 2]))))) - plot.default(1, 1, type = "n", xlim = xl, ylim = yl, xaxs = "r", + plot.default(NA, type = "n", xlim = xl, ylim = yl, xaxs = "r", yaxs = "r", xlab = xlab, ylab = ylab, log = log) lty <- if (!dlty) rep(1, n) else 1:n col <- if (!dcol) rep(1, n) else topo.colors(n) for (i in 1:n) - lines(TREES[[i]], col = col[i], lty = lty[i], type = "S") + lines(TREES[[i]], col = col[i], lty = lty[i], type = "S") if (legend) - legend(xl[1], yl[2], legend = names(TREES), - lty = lty, col = col, bty = "n") + legend(xl[1], yl[2], legend = names(TREES), + lty = lty, col = col, bty = "n") +} + +ltt.coplot <- function(phy, backward = TRUE, ...) +{ + layout(matrix(1:2, 2)) + par(mar = c(0, 3, 0.5, 0.5)) + o <- plot(phy, root.edge = TRUE, ...) + par(mar = c(3, 3, 0, 0.5)) + ltt.plot(phy, xlim = o$x.lim, backward = FALSE, xaxt = "n") + if (backward) axisPhylo() else axis(1) } diff --git a/man/ltt.plot.Rd b/man/ltt.plot.Rd index 6f09500..b40cf1e 100644 --- a/man/ltt.plot.Rd +++ b/man/ltt.plot.Rd @@ -2,24 +2,38 @@ \alias{ltt.plot} \alias{ltt.lines} \alias{mltt.plot} +\alias{ltt.coplot} +\alias{ltt.plot.coords} \title{Lineages Through Time Plot} +\description{ + These functions provide tools for plotting the numbers of lineages + through time from phylogenetic trees. +} \usage{ -ltt.plot(phy, xlab = "Time", ylab = "N", ...) -ltt.lines(phy, ...) +ltt.plot(phy, xlab = "Time", ylab = "N", + backward = TRUE, tol = 1e-6, ...) +ltt.lines(phy, backward = TRUE, tol = 1e-6, ...) mltt.plot(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE, - xlab = "Time", ylab = "N", log = "") + xlab = "Time", ylab = "N", log = "", backward = TRUE, + tol = 1e-6) +ltt.coplot(phy, backward = TRUE, ...) +ltt.plot.coords(phy, backward = TRUE, tol = 1e-6) } \arguments{ \item{phy}{an object of class \code{"phylo"}; this could be an object of class \code{"multiPhylo"} in the case of \code{mltt.plot}.} \item{xlab}{a character string (or a variable of mode character) - giving the label for the x-axis (default is "Time").} - \item{ylab}{idem for the y-axis (default is "N").} - \item{\dots}{in the cases of \code{ltt.plot()} and \code{ltt.lines()} - these are further (graphical) arguments to be passed to - \code{plot()} or \code{lines()}, respectively (see \code{Details:} - on how to transform the axes); in the case \code{mltt.plot()} these - are additional trees to be plotted (see \code{Details:}).} + giving the label for the \eqn{x}-axis (default is "Time").} + \item{ylab}{idem for the \eqn{y}-axis (default is "N").} + \item{backward}{a logical value: should the time axis be traced from + the present (the default), or from the root of the tree?} + \item{tol}{a numeric value (see details).} + \item{\dots}{in the cases of \code{ltt.plot()}, \code{ltt.lines()}, + or \code{ltt.coplot()} these are further (graphical) arguments to be + passed to \code{plot()}, \code{lines()}, or \code{plot.phylo()}, + respectively (see details on how to transform the axes); in + the case \code{mltt.plot()} these are additional trees to be plotted + (see details).} \item{dcol}{a logical specifying whether the different curves should be differentiated with colors (default is \code{TRUE}).} \item{dlty}{a logical specifying whether the different curves should @@ -31,13 +45,9 @@ mltt.plot(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE, log-transformed; must be one of the followings: \code{""}, \code{"x"}, \code{"y"}, or \code{"xy"}.} } -\description{ - These functions plot, on the current graphical device, the minimum - numbers of lineages through time from phylogenetic trees. -} \details{ \code{ltt.plot} does a simple lineages through time (LTT) - plot. Additional arguments (\code{...}) may be used to change, for + plot. Additional arguments (\code{\dots}) may be used to change, for instance, the limits on the axes (with \code{xlim} and/or \code{ylim}) or other graphical settings (\code{col} for the color, \code{lwd} for the line thickness, \code{lty} for the line type may be @@ -45,20 +55,40 @@ mltt.plot(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE, graphical parameters). The \eqn{y}-axis can be log-transformed by adding the following option: \code{log = "y"}. + The option \code{tol} is used as follows: first the most distant tip + from the root is found, then all tips whose distance to the root is + not more different from the previous one than \code{tol} are + considered to be contemporaneous with it. + + If the tree is not ultrametric, the plot is done assuming the tips, + except the most distant from the root, represent extinction events. If + a root edge is present, it is taken into account. + \code{ltt.lines} adds a LTT curve to an existing plot. Additional - arguments (\code{...}) may be used to change the settings of the added - line. Of course, the settings of the already existing LTT plot cannot - be altered this way. + arguments (\code{\dots}) may be used to change the settings of the added + line. \code{mltt.plot} does a multiple LTT plot taking as arguments one or several trees. These trees may be given as objects of class - \code{"phylo"} (single trees) or \code{"multiPhylo"} (multiple + \code{"phylo"} (single trees) and/or \code{"multiPhylo"} (multiple trees). Any number of objects may be given. This function is mainly for exploratory analyses with the advantages that the axes are set properly to view all lines, and the legend is plotted by default. For - more flexible settings of line drawings, it is probably better to - combine \code{ltt.plot()} with successive calls of \code{ltt.lines()} - (see \code{Examples:}). + more flexible settings of line drawings, it may be better to combine + \code{ltt.plot()} with successive calls of \code{ltt.lines()} (see + examples). + + \code{ltt.coplot} is meant to show how to set a tree and a LTT plots + on the same scales. All extra arguments modify only the appearance of + the tree. The code can be easily edited and tailored. +} +\value{ + \code{ltt.plot.coords} returns a two-column matrix with the time + points and the number of lineages, respectively. The \eqn{i}th value + of the second column is the number of lineages for the interval + defined by the \eqn{(i - 1)}th and the \eqn{i}th values of the first + column. These are then plotted with the option \code{type = "S"} by + the other functions. } \references{ Harvey, P. H., May, R. M. and Nee, S. (1994) Phylogenies without @@ -73,12 +103,11 @@ mltt.plot(phy, ..., dcol = TRUE, dlty = FALSE, legend = TRUE, \seealso{ \code{\link{skyline}}, \code{\link{branching.times}}, \code{\link{birthdeath}}, \code{\link{bd.ext}}, - \code{\link{yule.cov}}, \code{\link{bd.time}} + \code{\link{yule.cov}}, \code{\link{bd.time}}; \code{\link[graphics]{plot}} for the basic plotting function in R } \examples{ data(bird.families) -data(bird.orders) opar <- par(mfrow = c(2, 1)) ltt.plot(bird.families) title("Lineages Through Time Plot of the Bird Families") @@ -86,13 +115,33 @@ ltt.plot(bird.families, log = "y") title(main = "Lineages Through Time Plot of the Bird Families", sub = "(with logarithmic transformation of the y-axis)") par(opar) + ### to plot the tree and the LTT plot together +data(bird.orders) layout(matrix(1:4, 2, 2)) plot(bird.families, show.tip.label = FALSE) ltt.plot(bird.families, main = "Bird families") plot(bird.orders, show.tip.label = FALSE) ltt.plot(bird.orders, main = "Bird orders") -layout(matrix(1)) +layout(1) + +### better with ltt.coplot(): +ltt.coplot(bird.families, show.tip.label = FALSE, x.lim = 27.5) +data(chiroptera) +chiroptera <- compute.brlen(chiroptera) +ltt.coplot(chiroptera, show.tip.label = FALSE, type = "c") + +### with extinct lineages and a root edge: +omar <- par("mar") +set.seed(31) +tr <- rlineage(0.2, 0.15) +tr$root.edge <- 5 +ltt.coplot(tr, show.tip.label = FALSE, x.lim = 55) +## compare with: +## ltt.coplot(drop.fossil(tr), show.tip.label = FALSE) +layout(1) +par(mar = omar) + mltt.plot(bird.families, bird.orders) ### Generates 10 random trees with 23 tips: TR <- replicate(10, rcoal(23), FALSE) @@ -105,8 +154,8 @@ mltt.plot(TR, bird.orders) ### And now for something (not so) completely different: ltt.plot(bird.orders, lwd = 2) for (i in 1:10) ltt.lines(TR[[i]], lty = 2) -legend(-10, 5, lwd = c(2, 1), lty = c(1, 2), bty = "n", - legend = c("Bird orders", "Random trees")) +legend(-20, 10, lwd = c(2, 1), lty = c(1, 2), bty = "n", + legend = c("Bird orders", "Random (coalescent) trees")) } \keyword{hplot} \keyword{aplot} diff --git a/man/yule.time.Rd b/man/yule.time.Rd index 91cbe04..5092c2f 100644 --- a/man/yule.time.Rd +++ b/man/yule.time.Rd @@ -47,6 +47,11 @@ yule.time(phy, birth, BIRTH = NULL, root.time = 0, opti = "nlm", start = 0.01) An object of class \code{"yule"} (see \code{\link{yule}}). } \author{Emmanuel Paradis} +\references{ + Hubert, N., Paradis, E., Bruggemann, H. & Planes, S. (2011) Community + assembly and diversification in Indo-Pacific coral reef + fishes. \emph{Ecology and Evolution} (in press). +} \seealso{ \code{\link{branching.times}}, \code{\link{ltt.plot}}, \code{\link{birthdeath}}, \code{\link{yule}}, \code{\link{yule.cov}}, -- 2.39.2