new speciesTree()

git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@133 6e262413-ae40-0410-9e79-b911bd7a66b7

diff --git a/ChangeLog b/ChangeLog
index b4a7e8a377c17ed2308d424dfcf3f608363c71e8..16201fb94ce6683b5f4aad5cdb2e2f644f1bac6c 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,6 +1,13 @@
                CHANGES IN APE VERSION 2.6-1
 
 
+NEW FEATURES
+
+    o The new speciesTree calculates the species tree from a set of gene
+      trees. Two methods are currently available: maximum tree and
+      shallowest divergence tree.
+
+
 BUG FIXES
 
     o as.list.DNAbin() did not work correctly with vectors.

diff --git a/DESCRIPTION b/DESCRIPTION
index e165f27fbecb9593f8f2ab6427abb0394c943d79..27ce5f30651b33a89927f550a3c42dcfdf0e5730 100644 (file)
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ape
 Version: 2.6-1
-Date: 2010-10-19
+Date: 2010-10-20
 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 <Emmanuel.Paradis@ird.fr>

diff --git a/R/speciesTree.R b/R/speciesTree.R
new file mode 100644 (file)
index 0000000..ecaa63c
--- /dev/null
+++ b/R/speciesTree.R
@@ -0,0 +1,20 @@
+speciesTree <- function(x, FUN = min)
+### FUN = min => MAXTREE (Liu et al. 2010)
+### FUN = sum => shallowest divergence (Maddison & Knowles 2006)
+{
+    test.ultra <- which(!unlist(lapply(x, is.ultrametric)))
+    if (length(test.ultra))
+        stop(paste("the following trees were not ultrametric:\n",
+                   paste(test.ultra, collapse = " ")))
+
+    Ntree <- length(x)
+    D <- lapply(x, cophenetic.phylo)
+    nms <- rownames(D[[1]])
+    n <- length(nms)
+    M <- matrix(0, n*(n - 1)/2, Ntree)
+    for (i in 1:Ntree) M[, i] <- as.dist(D[[i]][nms, nms])
+    Y <- apply(M, 1, FUN)
+    attributes(Y) <- list(Size = n, Labels = nms, Diag = FALSE,
+                          Upper = FALSE, class = "dist")
+    as.phylo(stats::hclust(Y, "single"))
+}

diff --git a/man/speciesTree.Rd b/man/speciesTree.Rd
new file mode 100644 (file)
index 0000000..4ecbf21
--- /dev/null
+++ b/man/speciesTree.Rd
@@ -0,0 +1,57 @@
+\name{speciesTree}
+\alias{speciesTree}
+\title{Species Tree Estimation}
+\description{
+  This function calculates the species tree from a set of gene trees.
+}
+\usage{
+speciesTree(x, FUN = min)
+}
+\arguments{
+  \item{x}{a list of trees, e.g., an object of class
+    \code{"multiPhylo"}.}
+  \item{FUN}{a function used to compute the divergence times of each
+    pair of tips.}
+}
+\details{
+  For all trees in \code{x}, the divergence time of each pair of tips is
+  calculated: these are then `summarized' with \code{FUN} to build a new
+  distance matrix used to calculate the species tree with a
+  single-linkage hierarchical clustering. The default for \code{FUN}
+  computes the maximum tree (maxtree) of Liu et al. (2010). Using
+  \code{FUN = sum} gives the shallowest divergence tree of Maddison and
+  Knowles (2006).
+}
+\value{
+  an object of class \code{"phylo"}.
+}
+\references{
+  Liu, L., Yu, L. and Pearl, D. K. (2010) Maximum tree: a consistent
+  estimator of the species tree. \emph{Journal of Mathematical Biology},
+  \bold{60}, 95--106.
+
+  Maddison, W. P. and Knowles, L. L. (2006) Inferring phylogeny despite
+  incomplete lineage sorting. \emph{Systematic Biology}, \bold{55}, 21--30.
+}
+\author{Emmanuel Paradis}
+\examples{
+### example in Liu et al. (2010):
+tr1 <- read.tree(text = "(((B:0.05,C:0.05):0.01,D:0.06):0.04,A:0.1);")
+tr2 <- read.tree(text = "(((A:0.07,C:0.07):0.02,D:0.09):0.03,B:0.12);")
+TR <- c(tr1, tr2)
+TSmax <- speciesTree(TR) # MAXTREE
+TSsha <- speciesTree(TR, sum) # shallowest divergence
+
+layout(matrix(1:4, 1))
+## playing with 'x.lim' is not so complicated
+## but this will be improved someday
+plot(tr1, "c", d = "u", y.lim = c(-0.07, 0.1), font = 1)
+axisPhylo(4); title("Gene tree 1")
+plot(tr2, "c", d = "u", y.lim = c(-0.05, 0.12), font = 1)
+axisPhylo(4); title("Gene tree 2")
+plot(TSmax, "c", d = "u", y.lim = c(-0.1, 0.07), font = 1)
+axisPhylo(4); title("Species tree inferred\nby MAXTREE")
+plot(TSsha, "c", d = "u", y.lim = c(0, 0.17), font = 1)
+axisPhylo(4); title("Species tree inferred\nby Shallowest Divergence")
+}
+\keyword{models}