Package: ape
Version: 3.0-8
-Date: 2013-01-31
+Date: 2013-03-30
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 <Emmanuel.Paradis@ird.fr>
NEW FEATURES
+ o The new function ewLasso whether tests an incomplete set of
+ distances uniquely determines the edge weights of a given
+ unrooted topology using the 'Lasso' method by Dress et
+ al. (2012, J. Math. Biol. 65:77).
+
o ace() gains a new option 'use.expm' to use expm() from the
package of the same name in place of matexpo().
OTHER CHANGES
o The files CDAM.global.R and CDAM.post.R have been renamed
- CADM.global.R and CADM.post.R
+ CADM.global.R and CADM.post.R.
+
+ o ace() has a new default for its option 'method': this is "REML"
+ for continuous characters and "ML" for discrete ones.
+
+ o ape now imports the package expm so this one must be installed.
--- /dev/null
+## ewLasso.R (2013-03-30)
+
+## Lasso Tree
+
+## Copyright 2013 Andrei-Alin Popescu
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+ewLasso <- function(X, phy)
+{
+ if (is.matrix(X)) X <- as.dist(X)
+ X[is.na(X)] <- -1
+ X[X < 0] <- -1
+ X[is.nan(X)] <- -1
+
+ N <- attr(X, "Size")
+ labels <- attr(X, "Labels")
+ if (is.null(labels)) labels <- as.character(1:N)
+ ans <- .C("ewLasso", as.double(X), as.integer(N),
+ phy$edge[, 1], phy$edge[, 2],
+ DUP = FALSE, NAOK = TRUE, PACKAGE = "ape")
+}
-## plotPhyloCoor.R (2012-02-14)
+## plotPhyloCoor.R (2013-03-30)
## Coordinates of a Tree Plot
-## Copyright 2008 Damien de Vienne
+## Copyright 2008 Damien de Vienne, 2013 Klaus Schliep
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
plotPhyloCoor <-
function (x, type = "phylogram", use.edge.length = TRUE, node.pos = NULL,
- direction = "rightwards", ...)
+ direction = "rightwards", tip.order = NULL, ...)
{
Ntip <- length(x$tip.label)
if (Ntip == 1)
phyloORclado <- type %in% c("phylogram", "cladogram")
horizontal <- direction %in% c("rightwards", "leftwards")
if (phyloORclado) {
- if (!is.null(attr(x, "order")))
- if (attr(x, "order") == "pruningwise")
- x <- reorder(x)
+ ## changed by KS:
yy <- numeric(Ntip + Nnode)
- TIPS <- x$edge[x$edge[, 2] <= Ntip, 2]
- yy[TIPS] <- 1:Ntip
-
+ if (!is.null(tip.order)) {
+ yy[tip.order] <- 1:length(tip.order)
+ } else {
+ x <- reorder(x)
+ TIPS <- x$edge[x$edge[, 2] <= Ntip, 2]
+ yy[TIPS] <- 1:Ntip
+ }
}
xe <- x$edge
## first reorder the tree in cladewise order to avoid cophyloplot() hanging:
- x <- reorder(reorder(x), order = "pruningwise")
+ ## x <- reorder(reorder(x), order = "pruningwise") ... maybe not needed anymore (EP)
+ x <- reorder(x, order = "postorder")
ereorder <- match(x$edge[, 2], xe[, 2])
if (phyloORclado) {
If \code{type = "continuous"}, the default model is Brownian motion
where characters evolve randomly following a random walk. This model
can be fitted by residual maximum likelihood (the default), maximum
- likelihood (Schluter et al. 1997), least squares (\code{method =
- "pic"}, Felsenstein 1985), or generalized least squares (\code{method
- = "GLS"}, Martins and Hansen 1997, Cunningham et al. 1998). In the
- last case, the specification of \code{phy} and \code{model} are
- actually ignored: it is instead given through a correlation structure
- with the option \code{corStruct}.
+ likelihood (Felsenstein 1973, Schluter et al. 1997), least squares
+ (\code{method = "pic"}, Felsenstein 1985), or generalized least
+ squares (\code{method = "GLS"}, Martins and Hansen 1997, Cunningham et
+ al. 1998). In the last case, the specification of \code{phy} and
+ \code{model} are actually ignored: it is instead given through a
+ correlation structure with the option \code{corStruct}.
In the setting \code{method = "ML"} and \code{model = "BM"} (this used
- to be the default until ape 3.0-7) the maximum likelihood estimation
- is done simultaneously on the ancestral values and the variance of the
- Brownian motion process; these estimates are then used to compute the
- confidence intervals in the standard way. The REML method first
- estimates the ancestral value at the root (aka, the phylogenetic
- mean), then the variance of the Brownian motion process is estimated
- by optimizing the residual log-likelihood. The ancestral values are
- finally inferred from the likelihood function giving these two
- parameters. If \code{method = "pic"} or \code{"GLS"}, the confidence
- intervals are computed using the expected variances under the model,
- so they depend only on the tree.
+ to be the default until \pkg{ape} 3.0-7) the maximum likelihood
+ estimation is done simultaneously on the ancestral values and the
+ variance of the Brownian motion process; these estimates are then used
+ to compute the confidence intervals in the standard way. The REML
+ method first estimates the ancestral value at the root (aka, the
+ phylogenetic mean), then the variance of the Brownian motion process
+ is estimated by optimizing the residual log-likelihood. The ancestral
+ values are finally inferred from the likelihood function giving these
+ two parameters. If \code{method = "pic"} or \code{"GLS"}, the
+ confidence intervals are computed using the expected variances under
+ the model, so they depend only on the tree.
It could be shown that, with a continous character, REML results in
unbiased estimates of the variance of the Brownian motion process
an equal-rates model (e.g., the first and third examples above),
\code{"ARD"} is an all-rates-different model (the second example), and
\code{"SYM"} is a symmetrical model (e.g., \code{matrix(c(0, 1, 2, 1,
- 0, 3, 2, 3, 0), 3)}). If a short-cut is used, the number of states
- is determined from the data.
+ 0, 3, 2, 3, 0), 3)}). If a short-cut is used, the number of states is
+ determined from the data.
With discrete characters it is necessary to compute the exponential of
the rate matrix. By default (and the only possible choice until
reappraisal. \emph{Trends in Ecology & Evolution}, \bold{13},
361--366.
+ Felsenstein, J. (1973) Maximum likelihood estimation
+ of evolutionary trees from continuous characters. \emph{American Journal of Human Genetics}, \bold{25}, 471--492.
+
Felsenstein, J. (1985) Phylogenies and the comparative
method. \emph{American Naturalist}, \bold{125}, 1--15.
\value{
A logical value, or a two-column matrix.
}
-\author{\enc{Benoît}{Benoit} \email{b.durand@alfort.AFSSA.FR}}
+\author{\enc{Benoît}{Benoit} Durand \email{b.durand@alfort.AFSSA.FR}}
\seealso{
\code{\link[base]{all.equal}} for the generic \R function
}
--- /dev/null
+\name{ewLasso}
+\alias{ewLasso}
+\title{
+ Incomplete distances and edge weights of unrooted topology
+}
+\description{
+ This function implements a method for checking whether an incomplete
+ set of distances satisfy certain conditions that might make it
+ uniquely determine the edge weights of a given topology, T. It prints
+ information about whether the graph with vertex set the set of leaves,
+ denoted by X, and edge set the set of non-missing distance pairs,
+ denoted by L, is connected or strongly non-bipartite. It then also
+ checks whether L is a triplet cover for T.
+}
+\usage{
+ewLasso(X, phy)
+}
+\arguments{
+ \item{X}{a distance matrix.}
+ \item{phy}{an unrooted tree of class \code{"phylo"}.}
+}
+\details{
+ Missing values must be represented by either \code{NA} or a negative value.
+
+ This implements a method for checking whether an incomplete set of
+ distances satisfies certain conditions that might make it uniquely
+ determine the edge weights of a given topology, T. It prints
+ information about whether the graph, G, with vertex set the set of
+ leaves, denoted by X, and edge set the set of non-missing distance
+ pairs, denoted by L, is connected or strongly non-bipartite. It also
+ checks whether L is a triplet cover for T. If G is not connected, then
+ T does not need to be the only topology satisfying the input
+ incomplete distances. If G is not strongly non-bipartite then the
+ edge-weights of the edges of T are not the unique ones for which the
+ input distance is satisfied. If L is a triplet cover, then the input
+ distance matrix uniquely determines the edge weights of T. See Dress
+ et al. (2012) for details.
+}
+\value{
+ NULL, the results are printed in the console.
+}
+\references{
+ Dress, A. W. M., Huber, K. T., and Steel, M. (2012) `Lassoing' a
+ phylogentic tree I: basic properties, shellings and covers.
+ \emph{Journal of Mathematical Biology}, \bold{65(1)}, 77--105.
+}
+\author{Andrei Popescu \email{niteloserpopescu@gmail.com}}
+\keyword{multivariate}
}
\description{
The two FastME functions (balanced and OLS) perform the
- Minimum Evolution algorithm of Desper and Gascuel (2002).
+ minimum evolution algorithm of Desper and Gascuel (2002).
}
\usage{
fastme.bal(X, nni = TRUE, spr = TRUE, tbr = TRUE)
\itemize{
\item{``star''}{a star (or comb) tree with a single internal node.}
- \item{``balanced''}{a fully balanced dichotmous rooted tree;
- \code{n} must be of power of 2 (2, 4, 8, \dots).}
+ \item{``balanced''}{a fully balanced dichotomous rooted tree;
+ \code{n} must be a power of 2 (2, 4, 8, \dots).}
\item{``left''}{a fully unbalanced rooted tree where the largest
clade is on the left-hand side when the tree is plotted upwards.}
\item{``right''}{same than above but in the other direction.}
}
\note{
Specifying a negative value for `nbcol' (meaning that the nucleotides
- are printed on a single line) gives the same result for the
+ are printed on a single line) gives the same output for the
interleaved and sequential formats.
The names of the sequences may be truncated with the function
\code{\link{makeLabel}}. In particular, Clustal is limited to 30
- characters, whereas PHYML seems limited to 99 characters.
-
- The FASTA format may be used as input for Clustal.
+ characters, and PHYML seems limited to 99 characters.
}
\value{
None (invisible `NULL').
This function writes trees in a file with the NEXUS format.
}
\details{
- If several trees are given, they must have all the same tip labels.
+ If several trees are given, they must all have the same tip labels.
If among the objects given some are not trees of class \code{"phylo"},
- they are simply skipped and not written to the file.
+ they are simply skipped and not written in the file.
}
\value{
None (invisible `NULL').
be equal to the number of tips of the tree + the number of nodes. The
order is the following: first the values for the tips in the same
order than for the labels, then the values for the nodes sequentially
- from the root to the most terminal nodes (i.e. in the order given by
+ from the root to the most terminal nodes (i.e., in the order given by
\code{phy$edge}).
}
--- /dev/null
+/* ewLasso.c 2013-03-30 */
+
+/* Copyright 2013 Andrei-Alin Popescu */
+
+/* This file is part of the R-package `ape'. */
+/* See the file ../COPYING for licensing issues. */
+
+#include "ape.h"
+
+int isTripletCover(int nmb, int n, int** s, int stat, int sSoFar[n], int* a)//number of sides, number of leaves, sides, side under consideration, set so far, lasso
+{
+ int ret=0;
+ if(stat==nmb)return 1;
+ int i=0;
+
+ for(i=1;i<=n;i++)
+ {
+ if(!s[stat][i])continue;//not in set
+ int sw=1, j;
+
+ for(j=1;j<=n;j++)//check if all distances to previous candidates are present
+ {
+ if(!sSoFar[j])continue;//not in set so far
+ if(!a[i*(n+1)+j]){//if not, then i is not a good candidate for this side
+ //Rprintf("failed to find distance between %i and %i, a[%i][%i]=%i \n",i,j,i,j,a[i*(n+1)+j]);
+ sw=0;
+ }
+ }
+
+ if(!sw)continue;//not all required distances are present
+
+ sSoFar[i]=1;//try choosing i as representative for the side
+ ret+=isTripletCover(nmb,n,s,stat+1,sSoFar,a)>0?1:0;//see if, with i chosen, we can find leaves in other sides to satisfy the triplet cover condition
+ sSoFar[i]=0;
+ }
+ return ret;
+}
+
+void ewLasso(double *D, int *N, int *e1, int *e2)
+{
+ double *S, Sdist, Ndist, *new_dist, A, B, smallest_S, x, y;
+ int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
+
+ n=*N;
+ int tCov=1;
+ int* a = (int*)R_alloc((n+1)*(n+1), sizeof(int));//adjacency matrix of G_{\cL} graph
+
+ ij=0;
+
+ for(i=1;i<=n;i++)
+ {
+ for(j=1;j<=n;j++)
+ {
+ if(D[give_index(i,j,n)]==-1)//if missing value then no edge between pair of taxa (i,j) in G
+ {
+ a[i*(n+1)+j]=a[j*(n+1)+i]=0;
+ }
+ else
+ {
+ a[i*(n+1)+j]=a[j*(n+1)+i]=1;// otherwise edge between pair of taxa (i,j) in G
+ }
+ }
+ }
+ //check for connectedness of G
+
+ int *q = (int*)R_alloc(2*n-1, sizeof(int));//BFS queue
+ int *v = (int*)R_alloc(2*n-1, sizeof(int));//visited?
+
+ int p=0,u=1;//p-head of queue, u- position after last loaded element
+
+ for(i=1;i<=n;i++)v[i]=-1;
+
+ int stNBipartite=1, con=1, comp=1;
+ int ini=1;
+ ij=0;
+ /*for(i=1;i<=n;i++)
+ {
+ for(j=1;j<=n;j++)
+ {
+ Rprintf("a[%i][%i]=%i ",i,j,a[i*(n+1)+j]);
+ }
+ Rprintf("\n");
+ }*/
+
+ while(comp)
+ {
+ q[p]=ini;
+ v[ini]=1;
+ comp=0;
+ int stNBipartiteLoc=0;//check if current connected component is bipartite
+ while(p<u)//BFS
+ {
+ int head=q[p];
+ //Rprintf("head: %i\n",head);
+ //Rprintf("unvisited neighbours: \n");
+ for(i=1;i<=n;i++)
+ {
+ if(i==head)continue;
+ if(a[i*(n+1)+head]==0)continue;
+ if(v[i]==v[head])//same color as vertex from which we visit-> not bipartite
+ {
+ stNBipartiteLoc=1;
+ }
+ if(v[i]!=-1)continue;
+ //Rprintf("vertex %i \n",i);
+ q[u++]=i;
+ v[i]=1-v[head];
+ }
+ p++;
+ }
+ stNBipartite*=stNBipartiteLoc;//anding strngly-non-bipartite over all connected components
+
+
+
+ //check if all vertices have been visited
+ for(int i=1;i<=n;i++)
+ {
+ if(v[i]==-1)
+ {
+ comp=1;
+ p=0;
+ u=1;
+ ini=i;
+ con=0;
+ break;
+ }
+ }
+ }
+
+ Rprintf("connected: %i\n",con);
+ Rprintf("strongly non-bipartite: %i\n",stNBipartite);
+
+ //finally check if \cL is triplet cover of T
+
+ //adjencency matrix of tree, 1 to n are leaves
+
+ int *at= R_alloc((2*n-1)*(2*n-1), sizeof(int));
+
+ for(i=1;i<=2*n-2;i++)
+ {
+ for(j=1;j<=2*n-2;j++)at[i*(2*n-1)+j]=0;
+ }
+
+ for(i=0;i<2*n-3;i++)
+ {
+ //Rprintf("e1[%i]=%i e2[%i]=%i \n",i,e1[i],i,e2[i]);
+ at[e1[i]*(2*n-1)+e2[i]]=at[e2[i]*(2*n-1)+e1[i]]=1;
+ }
+
+ /*for(i=1;i<2*n-1;i++)
+ {
+ for(j=1;j<2*n-1;j++)
+ {
+ Rprintf("at[%i][%i]=%i ",i,j,at[i*(2*n-1)+j]);
+ }
+ Rprintf("\n");
+ }*/
+
+ for(i=n+1;i<=2*n-2;i++)//for each interior vertex
+ {
+ for(j=1;j<2*n-1;j++)//reset queue and visited veectors
+ {
+ v[j]=-1;
+ q[j]=0;
+ }
+
+ v[i]=1;//'disconnect' graph at i
+ int *l=(int*)R_alloc(2*n-2, sizeof(int));//vertices adjacent to i
+
+ int nmb=0;//number of found adjacent vertices of i
+
+ for(j=1;j<=2*n-2;j++)//find adjacent vertices
+ {
+ if(at[i*(2*n-1)+j]==1)
+ {
+ l[nmb++]=j;
+ }
+ }
+
+ int** s=(int**)R_alloc(nmb,sizeof(int*));//set of leaves in each side, stored as presence/absence
+
+ for(j=0;j<nmb;j++)s[j]=(int*)R_alloc(n+1,sizeof(int));
+
+ for(j=0;j<nmb;j++)
+ {
+ for(int k=1;k<=n;k++)s[j][k]=0;
+ }
+
+ /*Rprintf("for %i \n",i);
+
+ for(j=0;j<nmb;j++)Rprintf("l[%i]= %i ",j,l[j]);
+
+ Rprintf("\n");*/
+
+ for(j=0;j<nmb;j++)//for each adjancet vertex, find its leafset
+ {
+ p=0;
+ u=1;
+ ini=l[j];
+ q[p]=ini;
+ v[ini]=1;
+ if(ini<=n)s[j][ini]=1;
+ comp=0;
+ int nbr=0;
+ while(p<u)//BFS
+ {
+ int head=q[p];
+ //Rprintf("head: %i\n",head);
+ //Rprintf("unvisited neighbours: \n");
+ for(nbr=1;nbr<=2*n-1;nbr++)
+ {
+ if(nbr==head)continue;
+ if(at[nbr*(2*n-1)+head]==0)continue;
+ if(v[nbr]!=-1)continue;
+ //Rprintf("vertex %i \n",nbr);
+ if(nbr<=n)s[j][nbr]=1;//if leaf, then set visited to true
+ q[u++]=nbr;
+ v[nbr]=1;
+ }
+ p++;
+ }
+ }
+
+ /*Rprintf("sides for %i\n",i);
+ for(j=0;j<nmb;j++)
+ {
+ int ii;
+ for(ii=1;ii<=n;ii++)
+ {
+ if(s[j][ii])Rprintf("%i ",ii);
+ }
+ Rprintf("\n");
+ }*/
+
+ int* sSoFar= (int*)R_alloc(n+1,sizeof(int));
+
+ for(j=1;j<=n;j++)sSoFar[j]=0;
+
+ tCov*=isTripletCover(nmb,n,s,0,sSoFar,a)>0?1:0;
+ }
+ Rprintf("is triplet cover? %i \n",tCov);
+}
projects / ape.git / commitdiff