From 17d005bd203296df262aa76d75d165e4953494c6 Mon Sep 17 00:00:00 2001
From: paradis <paradis@6e262413-ae40-0410-9e79-b911bd7a66b7>
Date: Wed, 20 Oct 2010 15:42:48 +0000
Subject: [PATCH] new speciesTree()

git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@133 6e262413-ae40-0410-9e79-b911bd7a66b7
---
 ChangeLog          |  7 ++++++
 DESCRIPTION        |  2 +-
 R/speciesTree.R    | 20 ++++++++++++++++
 man/speciesTree.Rd | 57 ++++++++++++++++++++++++++++++++++++++++++++++
 4 files changed, 85 insertions(+), 1 deletion(-)
 create mode 100644 R/speciesTree.R
 create mode 100644 man/speciesTree.Rd

diff --git a/ChangeLog b/ChangeLog
index b4a7e8a..16201fb 100644
--- 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 e165f27..27ce5f3 100644
--- 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
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
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}
-- 
2.39.2