Package 'ape'

Description

ape provides functions for reading, writing, manipulating, analysing, and simulating phylogenetic trees and DNA sequences, computing DNA distances, translating into AA sequences, estimating trees with distance-based methods, and a range of methods for comparative analyses and analysis of diversification. Functionalities are also provided for programming new phylogenetic methods.

The complete list of functions can be displayed with library(help = ape).

More information on ape can be found at https://emmanuelparadis.github.io.

Author(s)

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 <[email protected]>

References

Paradis, E. (2012) Analysis of Phylogenetics and Evolution with R (Second Edition). New York: Springer.

Paradis, E., Claude, J. and Strimmer, K. (2004) APE: analyses of phylogenetics and evolution in R language. Bioinformatics, 20, 289–290.

Popescu, A.-A., Huber, K. T. and Paradis, E. (2012) ape 3.0: new tools for distance based phylogenetics and evolutionary analysis in R. Bioinformatics, 28, 1536–1537.

Paradis, E. and Schliep, K. (2019) ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics, 35, 526–528.

Description

These functions help to create and manipulate AA sequences.

Usage

## S3 method for class 'AAbin'
print(x, ...)

## S3 method for class 'AAbin'
x[i, j, drop = FALSE]

## S3 method for class 'AAbin'
c(..., recursive = FALSE)

## S3 method for class 'AAbin'
rbind(...)
## S3 method for class 'AAbin'
cbind(..., check.names = TRUE, fill.with.Xs = FALSE,
             quiet = FALSE)

## S3 method for class 'AAbin'
as.character(x, ...)

## S3 method for class 'AAbin'
labels(object, ...)

## S3 method for class 'AAbin'
image(x, what, col, bg = "white", xlab = "", ylab = "",
      show.labels = TRUE, cex.lab = 1, legend = TRUE, grid = FALSE,
      show.aa = FALSE, aa.cex = 1, aa.font = 1, aa.col = "black",
      scheme = "Ape_AA",...)

as.AAbin(x, ...)
## S3 method for class 'character'
as.AAbin(x, ...)

## S3 method for class 'list'
as.AAbin(x, ...)

## S3 method for class 'AAString'
as.AAbin(x, ...)

## S3 method for class 'AAStringSet'
as.AAbin(x, ...)

## S3 method for class 'AAMultipleAlignment'
as.AAbin(x, ...)

## S3 method for class 'AAbin'
as.list(x, ...)

## S3 method for class 'AAbin'
as.matrix(x, ...)

## S3 method for class 'AAbin'
as.phyDat(x, ...)

dist.aa(x, pairwise.deletion = FALSE, scaled = FALSE)
AAsubst(x)

Arguments

Details

These functions help to manipulate amino acid sequences of class "AAbin". These objects are stored in vectors, matrices, or lists which can be manipulated with the usual [ operator.

There is a conversion function to and from characters.

The function dist.aa computes the number of AA differences between each pair of sequences in a matrix; this can be scaled by the sequence length. See the function dist.ml in phangorn for evolutionary distances with AA sequences.

The function AAsubst returns the indices of the polymorphic sites (similar to seg.sites for DNA sequences; see examples below).

The two functions cbind.AAbin and rbind.AAbin work in the same way than the similar methods for the class "DNAbin": see cbind.DNAbin for more explanations about their respective behaviours.

Value

an object of class "AAbin", "character", "dist", or "numeric", depending on the function.

Author(s)

Emmanuel Paradis, Franz Krah

See Also

read.FASTA, trans, alview

Examples

data(woodmouse)
AA <- trans(woodmouse, 2)
seg.sites(woodmouse)
AAsubst(AA)

Description

ace estimates ancestral character states, and the associated uncertainty, for continuous and discrete characters. If marginal = TRUE, a marginal estimation procedure is used. With this method, the likelihood values at a given node are computed using only the information from the tips (and branches) descending from this node.

The present implementation of marginal reconstruction for discrete characters does not calculate the most likely state for each node, integrating over all the possible states, over all the other nodes in the tree, in proportion to their probability. For more details, see the Note below.

logLik, deviance, and AIC are generic functions used to extract the log-likelihood, the deviance, or the Akaike information criterion of a fitted object. If no such values are available, NULL is returned.

anova is another generic function which is used to compare nested models: the significance of the additional parameter(s) is tested with likelihood ratio tests. You must ensure that the models are effectively nested (if they are not, the results will be meaningless). It is better to list the models from the smallest to the largest.

Usage

ace(x, phy, type = "continuous", method = if (type == "continuous")
   "REML" else "ML", CI = TRUE,
    model = if (type == "continuous") "BM" else "ER",
    scaled = TRUE, kappa = 1, corStruct = NULL, ip = 0.1,
    use.expm = FALSE, use.eigen = TRUE, marginal = FALSE)
## S3 method for class 'ace'
print(x, digits = 4, ...)
## S3 method for class 'ace'
logLik(object, ...)
## S3 method for class 'ace'
deviance(object, ...)
## S3 method for class 'ace'
AIC(object, ..., k = 2)
## S3 method for class 'ace'
anova(object, ...)

Arguments

Details

If 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 (Felsenstein 1973, Schluter et al. 1997), least squares (method = "pic", Felsenstein 1985), or generalized least squares (method = "GLS", Martins and Hansen 1997, Cunningham et al. 1998). In the last case, the specification of phy and model are actually ignored: it is instead given through a correlation structure with the option corStruct.

In the setting method = "ML" and 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 method = "pic" or "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 while ML gives a downward bias. Therefore the former is recommanded.

For discrete characters (type = "discrete"), only maximum likelihood estimation is available (Pagel 1994) (see MPR for an alternative method). The model is specified through a numeric matrix with integer values taken as indices of the parameters. The numbers of rows and of columns of this matrix must be equal, and are taken to give the number of states of the character. For instance, matrix(c(0, 1, 1, 0), 2) will represent a model with two character states and equal rates of transition, matrix(c(0, 1, 2, 0), 2) a model with unequal rates, matrix(c(0, 1, 1, 1, 0, 1, 1, 1, 0), 3) a model with three states and equal rates of transition (the diagonal is always ignored). There are short-cuts to specify these models: "ER" is an equal-rates model (e.g., the first and third examples above), "ARD" is an all-rates-different model (the second example), and "SYM" is a symmetrical model (e.g., 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.

By default, the likelihood of the different ancestral states of discrete characters are computed with a joint estimation procedure using a procedure similar to the one described in Pupko et al. (2000). If marginal = TRUE, a marginal estimation procedure is used (this was the only choice until ape 3.1-1). With this method, the likelihood values at a given node are computed using only the information from the tips (and branches) descending from this node. With the joint estimation, all information is used for each node. The difference between these two methods is further explained in Felsenstein (2004, pp. 259-260) and in Yang (2006, pp. 121-126). The present implementation of the joint estimation uses a “two-pass” algorithm which is much faster than stochastic mapping while the estimates of both methods are very close.

With discrete characters it is necessary to compute the exponential of the rate matrix. The only possibility until ape 3.0-7 was the function matexpo in ape. If use.expm = TRUE and use.eigen = FALSE, the function expm, in the package of the same name, is used. matexpo is faster but quite inaccurate for large and/or asymmetric matrices. In case of doubt, use the latter. Since ape 3.0-10, it is possible to use an eigen decomposition avoiding the need to compute the matrix exponential; see details in Lebl (2013, sect. 3.8.3). This is much faster and is now the default.

Since version 5.2 of ape, ace can take state uncertainty for discrete characters into account: this should be coded with R's NA only. More details:

https://www.mail-archive.com/[email protected]/msg05286.html

Value

an object of class "ace" with the following elements:

Note

Liam Revell points out that for discrete characters the ancestral likelihood values returned with marginal = FALSE are actually the marginal estimates, while setting marginal = TRUE returns the conditional (scaled) likelihoods of the subtree:

http://blog.phytools.org/2015/05/about-how-acemarginaltrue-does-not.html

Author(s)

Emmanuel Paradis, Ben Bolker

References

Cunningham, C. W., Omland, K. E. and Oakley, T. H. (1998) Reconstructing ancestral character states: a critical reappraisal. Trends in Ecology & Evolution, 13, 361–366.

Felsenstein, J. (1973) Maximum likelihood estimation of evolutionary trees from continuous characters. American Journal of Human Genetics, 25, 471–492.

Felsenstein, J. (1985) Phylogenies and the comparative method. American Naturalist, 125, 1–15.

Felsenstein, J. (2004) Inferring Phylogenies. Sunderland: Sinauer Associates.

Lebl, J. (2013) Notes on Diffy Qs: Differential Equations for Engineers. https://www.jirka.org/diffyqs/.

Martins, E. P. and Hansen, T. F. (1997) Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data. American Naturalist, 149, 646–667.

Pagel, M. (1994) Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters. Proceedings of the Royal Society of London. Series B. Biological Sciences, 255, 37–45.

Pupko, T., Pe'er, I, Shamir, R., and Graur, D. (2000) A fast algorithm for joint reconstruction of ancestral amino acid sequences. Molecular Biology and Evolution, 17, 890–896.

Schluter, D., Price, T., Mooers, A. O. and Ludwig, D. (1997) Likelihood of ancestor states in adaptive radiation. Evolution, 51, 1699–1711.

Yang, Z. (2006) Computational Molecular Evolution. Oxford: Oxford University Press.

See Also

MPR, corBrownian, compar.ou, anova

Reconstruction of ancestral sequences can be done with the package phangorn (see function ?ancestral.pml).

Examples

### Some random data...
data(bird.orders)
x <- rnorm(23)
### Compare the three methods for continuous characters:
ace(x, bird.orders)
ace(x, bird.orders, method = "pic")
ace(x, bird.orders, method = "GLS",
    corStruct = corBrownian(1, bird.orders))
### For discrete characters:
x <- factor(c(rep(0, 5), rep(1, 18)))
ans <- ace(x, bird.orders, type = "d")
#### Showing the likelihoods on each node:
plot(bird.orders, type = "c", FALSE, label.offset = 1)
co <- c("blue", "yellow")
tiplabels(pch = 22, bg = co[as.numeric(x)], cex = 2, adj = 1)
nodelabels(thermo = ans$lik.anc, piecol = co, cex = 0.75)

Description

This function adds a horizontal bar giving the scale of the branch lengths to a plot of a phylogenetic tree on the current graphical device.

Usage

add.scale.bar(x, y, length = NULL, ask = FALSE,
              lwd = 1, lcol = "black", ...)

Arguments

Details

By default, the bar is placed in a corner of the graph depending on the direction of the tree. Otherwise both x and y must be specified (if only one is given it is ignored).

The further arguments (...) are used to format the text. They may be font, cex, col, and so on (see examples below, and the help page on text).

The function locator may be used to determine the x and y arguments.

Author(s)

Emmanuel Paradis

See Also

plot.phylo, axisPhylo, locator

Examples

tr <- rtree(10)
layout(matrix(1:2, 2, 1))
plot(tr)
add.scale.bar()
plot(tr)
add.scale.bar(cex = 0.7, font = 2, col = "red")
layout(1)

Description

Fills missing entries from incomplete distance matrix using the additive or the ultrametric procedure (see reference for details).

Usage

additive(X)
ultrametric(X)

Arguments

Value

a distance matrix.

Author(s)

Andrei Popescu

References

Makarenkov, V. and Lapointe, F.-J. (2004) A weighted least-squares approach for inferring phylogenies from incomplete distance matrices. Bioinformatics, 20, 2113–2121.

Description

This function helps to explore DNA alignments by zooming in. The user clicks twice defining the opposite corners of the portion which is extracted and drawned on a new window.

Usage

alex(x, ...)

Arguments

Details

This function works with a DNA alignment (freshly) plotted on an interactive graphical device (i.e., not a file) with image. After calling alex, the user clicks twice defining a rectangle in the alignment, then this portion of the alignment is extacted and plotted on a new window. The user can click as many times on the alignment. The process is stopped by a right-click. If the user clicks twice outside the alignment, a message “Try again!” is printed.

Each time alex is called, the alignment is plotted on a new window without closing or deleting those possibly already plotted.

In all cases, the device where x is plotted is the active window after the operation. It should not be closed during the whole process.

Value

NULL

Author(s)

Emmanuel Paradis

See Also

image.DNAbin, trex, alview

Examples

## Not run: 
data(woodmouse)
image(woodmouse)
alex(woodmouse)

## End(Not run)

Description

Comparison of DNA sequence sets, particularly when aligned.

Usage

## S3 method for class 'DNAbin'
all.equal(target, current, plot = FALSE, ...)

Arguments

Details

If the two sets of DNA sequences are exactly identical, this function returns TRUE. Otherwise, a detailed comparison is made only if the labels (i.e., rownames) of target and current are the same (possibly in different orders). In all other cases, a brief description of the differences is returned (sometimes with recommendations to make further comparisons).

This function can be used for testing in programs using isTRUE (see examples below).

Value

TRUE if the two sets are identical; a list with two elements (message and different.sites) if a detailed comparison is done; or a vector of mode character.

Author(s)

Emmanuel Paradis

See Also

image.DNAbin, clustal, checkAlignment, the generic function: all.equal

Examples

data(woodmouse)
woodm2 <- woodmouse
woodm2[1, c(1:5, 10:12, 30:40)] <- as.DNAbin("g")
res <- all.equal(woodmouse, woodm2, plot = TRUE)
str(res)

## if used for testing in R programs:
isTRUE(all.equal(woodmouse, woodmouse)) # TRUE
isTRUE(all.equal(woodmouse, woodm2)) # FALSE

all.equal(woodmouse, woodmouse[15:1, ])
all.equal(woodmouse, woodmouse[-1, ])
all.equal(woodmouse, woodmouse[, -1])

## Not run: 
## To run the followings you need internet and Clustal and MUSCLE
## correctly installed.
## Data from Johnson et al. (2006, Science)
refs <- paste("DQ082", 505:545, sep = "")
DNA <- read.GenBank(refs)
DNA.clustal <- clustal(DNA)
DNA.muscle <- muscle(DNA)
isTRUE(all.equal(DNA.clustal, DNA.muscle)) # FALSE
all.equal(DNA.clustal, DNA.muscle, TRUE)

## End(Not run)

Description

This function makes a global comparison of two phylogenetic trees.

Usage

## S3 method for class 'phylo'
all.equal(target, current, use.edge.length = TRUE,
                   use.tip.label = TRUE, index.return = FALSE,
                   tolerance = .Machine$double.eps ^ 0.5,
                   scale = NULL, ...)

Arguments

Details

This function is meant to be an adaptation of the generic function all.equal for the comparison of phylogenetic trees.

A single phylogenetic tree may have several representations in the Newick format and in the "phylo" class of objects used in ‘ape’. One aim of the present function is to be able to identify whether two objects of class "phylo" represent the same phylogeny.

Value

A logical value, or a two-column matrix.

Note

The algorithm used here does not work correctly for the comparison of topologies (i.e., ignoring tip labels) of unrooted trees. This also affects unique.multiPhylo which calls the present function. See:

https://www.mail-archive.com/[email protected]/msg01445.html.

Author(s)

Benoît Durand [email protected]

See Also

all.equal for the generic R function, comparePhylo

Examples

### maybe the simplest example of two representations
### for the same rooted tree...:
t1 <- read.tree(text = "(a:1,b:1);")
t2 <- read.tree(text = "(b:1,a:1);")
all.equal(t1, t2)
### ... compare with this:
identical(t1, t2)
### one just slightly more complicated...:
t3 <- read.tree(text = "((a:1,b:1):1,c:2);")
t4 <- read.tree(text = "(c:2,(a:1,b:1):1);")
all.equal(t3, t4) # == all.equal.phylo(t3, t4)
### ... here we force the comparison as lists:
all.equal.list(t3, t4)

Description

This function displays in the console or a file an alignment of DNA or AAsequences. The first sequence is printed on the first row and the bases of the other sequences are replaced by dots if they are identical with the first sequence.

Usage

alview(x, file = "", uppercase = TRUE, showpos = TRUE)

Arguments

Details

The first line of the output shows the position of the last column of the printed alignment.

Author(s)

Emmanuel Paradis

See Also

DNAbin, image.DNAbin, alex, clustal, checkAlignment, all.equal.DNAbin

Examples

data(woodmouse)
alview(woodmouse[, 1:50])
alview(woodmouse[, 1:50], uppercase = FALSE)
## display only some sites:
j <- c(10, 49, 125, 567) # just random
x <- woodmouse[, j]
alview(x, showpos = FALSE) # no site position displayed
alview(x, showpos = j)
## Not run: 
alview(woodmouse, file = "woodmouse.txt")

## End(Not run)

Description

These functions help to find files on the local disk.

Usage

Xplorefiles(from = "HOME", recursive = TRUE, ignore.case = TRUE)
editFileExtensions()
bydir(x)
Xplor(from = "HOME")

Arguments

Details

Xplorefiles looks for all files with a specified extension in their names. The default is to look for the following file types: CLUSTAL (.aln), FASTA (.fas, .fasta), FASTQ (.fq, .fastq), NEWICK (.nwk, .newick, .tre, .tree), NEXUS (.nex, .nexus), and PHYLIP (.phy). This list can be modified with editFileExtensions.

bydir sorts the list of files by directories.

Xplor combines the other operations and opens the results in a Web browser with clickable links to the directories and files.

Value

Xplorefiles returns a list. bydir prints the file listings on the console.

Author(s)

Emmanuel Paradis

Examples

## Not run: 
x <- Xplorefiles()
x # all data files on your disk
bydir(x) # sorted by directories
bydir(x["fasta"]) # only the FASTA files
Xplorefiles(getwd(), recursive = FALSE) # look only in current dir
Xplor()

## End(Not run)

Description

These functions transform a set of DNA sequences among various internal formats.

Usage

as.alignment(x)
as.DNAbin(x, ...)

## S3 method for class 'character'
as.DNAbin(x, ...)

## S3 method for class 'list'
as.DNAbin(x, ...)

## S3 method for class 'alignment'
as.DNAbin(x, ...)

## S3 method for class 'DNAString'
as.DNAbin(x, ...)

## S3 method for class 'DNAStringSet'
as.DNAbin(x, ...)

## S3 method for class 'PairwiseAlignmentsSingleSubject'
as.DNAbin(x, ...)

## S3 method for class 'DNAMultipleAlignment'
as.DNAbin(x, ...)

## S3 method for class 'DNAbin'
as.character(x, ...)

Arguments

Details

For as.alignment, the sequences given as argument should be stored as matrices or lists of single-character strings (the format used in ape before version 1.10). The returned object is in the format used in the package seqinr to store aligned sequences.

as.DNAbin is a generic function with methods so that it works with sequences stored into vectors, matrices, or lists. It can convert some S4 classes from the package Biostrings in BioConductor. For consistency within ape, this uses an S3-style syntax. To convert objects of class "DNAStringSetList", see the examples.

as.character is a generic function: the present method converts objects of class "DNAbin" into the format used before ape 1.10 (matrix of single characters, or list of vectors of single characters). This function must be used first to convert objects of class "DNAbin" into the class "alignment".

Value

an object of class "alignment" in the case of "as.alignment"; an object of class "DNAbin" in the case of "as.DNAbin"; a matrix of mode character or a list containing vectors of mode character in the case of "as.character".

Author(s)

Emmanuel Paradis

See Also

DNAbin, read.dna, read.GenBank, write.dna

Examples

data(woodmouse)
x <- as.character(woodmouse)
x[, 1:20]
str(as.alignment(x))
identical(as.DNAbin(x), woodmouse)
### conversion from BioConductor:
## Not run: 
if (require(Biostrings)) {
data(phiX174Phage)
X <- as.DNAbin(phiX174Phage)
## base frequencies:
base.freq(X) # from ape
alphabetFrequency(phiX174Phage) # from Biostrings
### for objects of class "DNAStringSetList"
X <- lapply(x, as.DNAbin) # a list of lists
### to put all sequences in a single list:
X <- unlist(X, recursive = FALSE)
class(X) <- "DNAbin"
}

## End(Not run)

Description

bitsplits returns the bipartitions (aka splits) for a single tree or a list of trees. If at least one tree is rooted, an error is returned.

countBipartitions returns the frequencies of the bipartitions from a reference tree (phy) observed in a list of trees (X), all unrooted.

as.bitsplits and as.prop.part are generic functions for converting between the "bitsplits" and "prop.part" classes.

Usage

bitsplits(x)
countBipartitions(phy, X)
as.bitsplits(x)
## S3 method for class 'prop.part'
as.bitsplits(x)
## S3 method for class 'bitsplits'
print(x, ...)
## S3 method for class 'bitsplits'
sort(x, decreasing = FALSE, ...)
as.prop.part(x, ...)
## S3 method for class 'bitsplits'
as.prop.part(x, include.trivial = FALSE, ...)

Arguments

Details

These functions count bipartitions as defined by internal branches, so they work only with unrooted trees. The structure of the class "bitsplits" is described in a separate document on ape's web site.

This data structure has a memory requirement proportional to $n^2$, so it can be inefficient with large trees (> 1000 tips), particularly if they are very different (i.e., with few shared splits). In any case, an error occurs if the product of the number of tips by the number of nodes is greater than $2^{31}-1$ (~2.1 billion). A warning message is given if the tree(s) has(ve) more than 46,341 tips. It may happen that the search for splits is interrupted if the data structure is full (with a warning message).

Value

bitsplits, as.bitsplits, and sort return an object of class "bitsplits".

countBipartitions returns a vector of integers.

as.prop.part returns an object of class "prop.part".

Author(s)

Emmanuel Paradis

See Also

prop.part, is.compatible

Examples

tr <- rtree(20)
pp <- prop.part(tr)
as.bitsplits(pp)
## works only with unrooted trees (ape 5.5):
countBipartitions(rtree(10, rooted = FALSE), rmtree(100, 10, rooted = FALSE))

Description

These functions convert objects between the classes "phylo" and "matching".

Usage

as.matching(x, ...)
## S3 method for class 'phylo'
as.matching(x, labels = TRUE, ...)
## S3 method for class 'matching'
as.phylo(x, ...)

Arguments

Details

A matching is a representation where each tip and each node are given a number, and sibling groups are grouped in a “matching pair” (see Diaconis and Holmes 1998, for details). This coding system can be used only for binary (fully dichotomous) trees.

Diaconis and Holmes (1998) gave some conventions to insure that a given tree has a unique representation as a matching. I have tried to follow them in the present functions.

Value

as.matching returns an object of class "matching" with the following component:

as.phylo.matching returns an object of class "phylo".

Note

Branch lengths are not supported in the present version.

Author(s)

Emmanuel Paradis

References

Diaconis, P. W. and Holmes, S. P. (1998) Matchings and phylogenetic trees. Proceedings of the National Academy of Sciences USA, 95, 14600–14602.

See Also

as.phylo

Examples

data(bird.orders)
m <- as.matching(bird.orders)
str(m)
m
tr <- as.phylo(m)
all.equal(tr, bird.orders, use.edge.length = FALSE)

Description

as.phylo is a generic function which converts an object into a tree of class "phylo". There are currently two methods for objects of class "hclust" and of class "phylog" (implemented in the package ade4). The default method is for any object inheriting the class "phylo" which is returned unchanged.

as.hclust.phylo is a method of the generic as.hclust which converts an object of class "phylo" into one of class "hclust". This can used to convert an object of class "phylo" into one of class "dendrogram" (see examples).

as.network and as.igraph convert trees of class "phylo" into these respective classes defined in the packages of the same names (where the generics are defined).

old2new.phylo and new2old.phylo are utility functions for converting between the old and new coding of the class "phylo".

Usage

as.phylo(x, ...)
## Default S3 method:
as.phylo(x, ...)
## S3 method for class 'hclust'
as.phylo(x, ...)
## S3 method for class 'phylog'
as.phylo(x, ...)
## S3 method for class 'phylo'
as.hclust(x, ...)
old2new.phylo(phy)
new2old.phylo(phy)
## S3 method for class 'phylo'
as.network(x, directed = is.rooted(x), ...)
## S3 method for class 'phylo'
as.igraph(x, directed = is.rooted(x), use.labels = TRUE, ...)

Arguments

Value

An object of class "hclust", "phylo", "network", or "igraph".

Note

In an object of class "hclust", the height gives the distance between the two sets that are being agglomerated. So these distances are divided by two when setting the branch lengths of a phylogenetic tree.

Author(s)

Emmanuel Paradis

See Also

hclust, as.hclust, dendrogram, as.phylo.formula

Examples

data(bird.orders)
hc <- as.hclust(bird.orders)
tr <- as.phylo(hc)
all.equal(bird.orders, tr) # TRUE

### shows the three plots for tree objects:
dend <- as.dendrogram(hc)
layout(matrix(c(1:3, 3), 2, 2))
plot(bird.orders, font = 1)
plot(hc)
par(mar = c(8, 0, 0, 0)) # leave space for the labels
plot(dend)

### how to get identical plots with
### plot.phylo and plot.dendrogram:
layout(matrix(1:2, 2, 1))
plot(bird.orders, font = 1, no.margin = TRUE, label.offset = 0.4)
par(mar = c(0, 0, 0, 8))
plot(dend, horiz = TRUE)
layout(1)

## Not run: 
### convert into networks:
if (require(network)) {
    x <- as.network(rtree(10))
    print(x)
    plot(x, vertex.cex = 1:4)
    plot(x, displaylabels = TRUE)
}
tr <- rtree(5)
if (require(igraph)) {
    print((x <- as.igraph(tr)))
    plot(x)
    print(as.igraph(tr, TRUE, FALSE))
    print(as.igraph(tr, FALSE, FALSE))
}

## End(Not run)

Description

The function as.phylo.formula (short form as.phylo) builds a phylogenetic tree (an object of class phylo) from a set of nested taxonomic variables.

Usage

## S3 method for class 'formula'
as.phylo(x, data = parent.frame(), collapse = TRUE, ...)

Arguments

Details

Taxonomic variables must be nested and passed in the correct order: the higher clade must be on the left of the formula, for instance ~Order/Family/Genus/Species. In most cases, the resulting tree will be unresolved and will contain polytomies.

The option collapse = FALSE has for effect to add single nodes in the tree when a given higher level has only one element in the level below (e.g., a monospecific genus); see the example below.

Value

an object of class "phylo".

Author(s)

Julien Dutheil [email protected], Eric Marcon and Klaus Schliep

See Also

as.phylo, read.tree for a description of "phylo" objects, multi2di

Examples

data(carnivora)
frm <- ~SuperFamily/Family/Genus/Species
tr <- as.phylo(frm, data = carnivora, collapse=FALSE)
tr$edge.length <- rep(1, nrow(tr$edge))
plot(tr, show.node.label=TRUE)
Nnode(tr)
## compare with:
Nnode(as.phylo(frm, data = carnivora, collapse = FALSE))

Description

This function adds a scaled axis on the side of a phylogeny plot.

Usage

axisPhylo(side = NULL, root.time = NULL, backward = TRUE, ...)

Arguments

Details

The further arguments (...) are used to format the axis. They may be font, cex, col, las, and so on (see the help pages on axis and par).

Author(s)

Emmanuel Paradis, Klaus Schliep

See Also

plot.phylo, add.scale.bar, axis, par

Examples

tr <- rtree(30)
ch <- rcoal(30)
plot(ch)
axisPhylo()
plot(tr, "c", FALSE, direction = "u")
axisPhylo(las = 1)

Description

This function computes the balance of a phylogenetic tree, that is for each node of the tree the numbers of descendants (i.e. tips) on each of its daughter-branch. The tree must be fully dichotomous.

Usage

balance(phy)

Arguments

Value

a numeric matrix with two columns and one row for each node of the tree. The columns give the numbers of descendants on each daughter-branches (the order of both columns being arbitrary). If the phylogeny phy has an element node.label, this is used as rownames for the returned matrix; otherwise the numbers (of mode character) of the matrix edge of phy are used as rownames.

Author(s)

Emmanuel Paradis

References

Aldous, D. J. (2001) Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today. Statistical Science, 16, 23–34.

Description

base.freq computes the frequencies (absolute or relative) of the four DNA bases (adenine, cytosine, guanine, and thymidine) from a sample of sequences.

GC.content computes the proportion of G+C (using the previous function). All missing or unknown sites are ignored.

Ftab computes the contingency table with the absolute frequencies of the DNA bases from a pair of sequences.

Usage

base.freq(x, freq = FALSE, all = FALSE)
GC.content(x)
Ftab(x, y = NULL)

Arguments

Details

The base frequencies are computed over all sequences in the sample.

For Ftab, if the argument y is given then both x and y are coerced as vectors and must be of equal length. If y is not given, x must be a matrix or a list and only the two first sequences are used.

Value

A numeric vector with names c("a", "c", "g", "t") (and possibly "r", "m", ..., a single numeric value, or a four by four matrix with similar dimnames.

Author(s)

Emmanuel Paradis

See Also

seg.sites, nuc.div (in pegas), DNAbin

Examples

data(woodmouse)
base.freq(woodmouse)
base.freq(woodmouse, TRUE)
base.freq(woodmouse, TRUE, TRUE)
GC.content(woodmouse)
Ftab(woodmouse)
Ftab(woodmouse[1, ], woodmouse[2, ]) # same than above
Ftab(woodmouse[14:15, ]) # between the last two

Description

This function fits by maximum likelihood a birth-death model to the combined phylogenetic and taxonomic data of a given clade. The phylogenetic data are given by a tree, and the taxonomic data by the number of species for the its tips.

Usage

bd.ext(phy, S, conditional = TRUE)

Arguments

Details

A re-parametrization of the birth-death model studied by Kendall (1948) so that the likelihood has to be maximized over d/b and b - d, where b is the birth rate, and d the death rate.

The standard-errors of the estimated parameters are computed using a normal approximation of the maximum likelihood estimates.

If the argument S has names, then they are matched to the tip labels of phy. The user must be careful here since the function requires that both series of names perfectly match, so this operation may fail if there is a typing or syntax error. If both series of names do not match, the values S are taken to be in the same order than the tip labels of phy, and a warning message is issued.

Note that the function does not check that the tree is effectively ultrametric, so if it is not, the returned result may not be meaningful.

If conditional = TRUE, the probabilities of the taxonomic data are calculated conditioned on no extinction (Rabosky et al. 2007). In previous versions of the present function (until ape 2.6-1), unconditional probabilities were used resulting in underestimated extinction rate. Though it does not make much sense to use conditional = FALSE, this option is provided to compare results from previous analyses: if the species richnesses are relatively low, both versions will give similar results (see examples).

Author(s)

Emmanuel Paradis

References

Paradis, E. (2003) Analysis of diversification: combining phylogenetic and taxonomic data. Proceedings of the Royal Society of London. Series B. Biological Sciences, 270, 2499–2505.

Rabosky, D. L., Donnellan, S. C., Talaba, A. L. and Lovette, I. J. (2007) Exceptional among-lineage variation in diversification rates during the radiation of Australia's most diverse vertebrate clade. Proceedings of the Royal Society of London. Series B. Biological Sciences, 274, 2915–2923.

See Also

birthdeath, branching.times, diversi.gof, diversi.time, ltt.plot, yule, yule.cov, bd.time

Examples

### An example from Paradis (2003) using the avian orders:
data(bird.orders)
### Number of species in each order from Sibley and Monroe (1990):
S <- c(10, 47, 69, 214, 161, 17, 355, 51, 56, 10, 39, 152,
       6, 143, 358, 103, 319, 23, 291, 313, 196, 1027, 5712)
bd.ext(bird.orders, S)
bd.ext(bird.orders, S, FALSE) # same than older versions

Description

This function fits a used-defined time-dependent birth-death model.

Usage

bd.time(phy, birth, death, BIRTH = NULL, DEATH = NULL,
        ip, lower, upper, fast = FALSE, boot = 0, trace = 0)

Arguments

Details

Details on how to specify the birth and death functions and their primitives can be found in the help page of yule.time.

The model is fitted by minimizing the least squares deviation between the observed and the predicted distributions of branching times. These computations rely heavily on numerical integrations. If fast = FALSE, integrations are done with R's integrate function. If fast = TRUE, a faster but less accurate function provided in ape is used. If fitting a complex model to a large phylogeny, a strategy might be to first use the latter option, and then to use the estimates as starting values with fast = FALSE.

Value

A list with the following components:

• par: a vector of estimates with names taken from the parameters in the specified functions.

• SS: the minimized sum of squares.

• convergence: output convergence criterion from nlminb.

• message: id.

• iterations: id.

• evaluations: id.

Author(s)

Emmanuel Paradis

References

Paradis, E. (2011) Time-dependent speciation and extinction from phylogenies: a least squares approach. Evolution, 65, 661–672.

See Also

ltt.plot, birthdeath, yule.time, LTT

Examples

set.seed(3)
tr <- rbdtree(0.1, 0.02)
bd.time(tr, 0, 0) # fits a simple BD model
bd.time(tr, 0, 0, ip = c(.1, .01)) # 'ip' is useful here
## the classic logistic:
birth.logis <- function(a, b) 1/(1 + exp(-a*t - b))
## Not run: 
bd.time(tr, birth.logis, 0, ip = c(0, -2, 0.01))
## slow to get:
## $par
##            a            b        death
## -0.003486961 -1.995983179  0.016496454
##
## $SS
## [1] 20.73023

## End(Not run)

Description

binaryPGLMM performs linear regression for binary phylogenetic data, estimating regression coefficients with approximate standard errors. It simultaneously estimates the strength of phylogenetic signal in the residuals and gives an approximate conditional likelihood ratio test for the hypothesis that there is no signal. Therefore, when applied without predictor (independent) variables, it gives a test for phylogenetic signal for binary data. The method uses a GLMM approach, alternating between penalized quasi-likelihood (PQL) to estimate the "mean components" and restricted maximum likelihood (REML) to estimate the "variance components" of the model.

binaryPGLMM.sim is a companion function that simulates binary phylogenetic data of the same structure analyzed by binaryPGLMM.

Usage

binaryPGLMM(formula, data = list(), phy, s2.init = 0.1,
            B.init = NULL, tol.pql = 10^-6, maxit.pql = 200,
            maxit.reml = 100)

binaryPGLMM.sim(formula, data = list(), phy, s2 = NULL, B = NULL, nrep = 1)

## S3 method for class 'binaryPGLMM'
print(x, digits = max(3, getOption("digits") - 3), ...)

Arguments

Details

The function estimates parameters for the model

$Pr(Y = 1) = q$

$q = inverse.logit(b0 + b1 * x1 + b2 * x2 + \dots + \epsilon)$

$\epsilon ~ Gaussian(0, s2 * V)$

where $V$ is a variance-covariance matrix derived from a phylogeny (typically under the assumption of Brownian motion evolution). Although mathematically there is no requirement for $V$ to be ultrametric, forcing $V$ into ultrametric form can aide in the interpretation of the model, because in regression for binary dependent variables, only the off-diagonal elements (i.e., covariances) of matrix $V$ are biologically meaningful (see Ives & Garland 2014).

The function converts a phylo tree object into a variance-covariance matrix, and further standardizes this matrix to have determinant = 1. This in effect standardizes the interpretation of the scalar s2. Although mathematically not required, it is a very good idea to standardize the predictor (independent) variables to have mean 0 and variance 1. This will make the function more robust and improve the interpretation of the regression coefficients. For categorical (factor) predictor variables, you will need to construct 0-1 dummy variables, and these should not be standardized (for obvious reasons).

The estimation method alternates between PQL to obtain estimates of the mean components of the model (this is the standard approach to estimating GLMs) and REML to obtain estimates of the variance components. This method gives relatively fast and robust estimation. Nonetheless, the estimates of the coefficients B will generally be upwards bias, as is typical of estimation for binary data. The standard errors of B are computed from the PQL results conditional on the estimate of s2 and therefore should tend to be too small. The function returns an approximate P-value for the hypothesis of no phylogenetic signal in the residuals (i.e., H0:s2 = 0) using an approximate likelihood ratio test based on the conditional REML likelihood (rather than the marginal likelihood). Simulations have shown that these P-values tend to be high (giving type II errors: failing to identify variances that in fact are statistically significantly different from zero).

It is a good idea to confirm statistical inferences using parametric bootstrapping, and the companion function binaryPGLMM.sim gives a simply tool for this. See Examples below.

Value

An object of class "binaryPGLMM".

Author(s)

Anthony R. Ives

References

Ives, A. R. and Helmus, M. R. (2011) Generalized linear mixed models for phylogenetic analyses of community structure. Ecological Monographs, 81, 511–525.

Ives, A. R. and Garland, T., Jr. (2014) Phylogenetic regression for binary dependent variables. Pages 231–261 in L. Z. Garamszegi, editor. Modern Phylogenetic Comparative Methods and Their Application in Evolutionary Biology. Springer-Verlag, Berlin Heidelberg.

See Also

package pez and its function communityPGLMM; package phylolm and its function phyloglm; package MCMCglmm

Examples

## Illustration of binaryPGLMM() with simulated data

# Generate random phylogeny

n <- 100
phy <- compute.brlen(rtree(n=n), method = "Grafen", power = 1)

# Generate random data and standardize to have mean 0 and variance 1
X1 <- rTraitCont(phy, model = "BM", sigma = 1)
X1 <- (X1 - mean(X1))/var(X1)

# Simulate binary Y
sim.dat <- data.frame(Y=array(0, dim=n), X1=X1, row.names=phy$tip.label)
sim.dat$Y <- binaryPGLMM.sim(Y ~ X1, phy=phy, data=sim.dat, s2=.5,
                             B=matrix(c(0,.25),nrow=2,ncol=1), nrep=1)$Y

# Fit model
binaryPGLMM(Y ~ X1, phy=phy, data=sim.dat)

## Not run: 
# Compare with phyloglm()
library(phylolm)
summary(phyloglm(Y ~ X1, phy=phy, data=sim.dat))

# Compare with glm() that does not account for phylogeny
summary(glm(Y ~ X1, data=sim.dat, family="binomial"))

# Compare with logistf() that does not account
# for phylogeny but is less biased than glm()
library(logistf)
logistf(Y ~ X1, data=sim.dat)

# Compare with MCMCglmm
library(MCMCglmm)

V <- vcv(phy)
V <- V/max(V)
detV <- exp(determinant(V)$modulus[1])
V <- V/detV^(1/n)

invV <- Matrix(solve(V),sparse=T)
sim.dat$species <- phy$tip.label
rownames(invV) <- sim.dat$species

nitt <- 43000
thin <- 10
burnin <- 3000

prior <- list(R=list(V=1, fix=1), G=list(G1=list(V=1, nu=1000, alpha.mu=0, alpha.V=1)))
summary(MCMCglmm(Y ~ X1, random=~species, ginvers=list(species=invV),
    data=sim.dat, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,
    family="categorical", prior=prior, verbose=FALSE))

## Examine bias in estimates of B1 and s2 from binaryPGLMM with
# simulated data. Note that this will take a while.

Reps = 1000

s2 <- 0.4
B1 <- 1

meanEsts <- data.frame(n = Inf, B1 = B1, s2 = s2, Pr.s2 = 1, propconverged = 1)

for (n in c(160, 80, 40, 20)) {

  meanEsts.n <- data.frame(B1 = 0, s2 = 0, Pr.s2 = 0, convergefailure = 0)
    for (rep in 1:Reps) {
      phy <- compute.brlen(rtree(n = n), method = "Grafen", power = 1)
      X <- rTraitCont(phy, model = "BM", sigma = 1)
      X <- (X - mean(X))/var(X)

      sim.dat <- data.frame(Y = array(0, dim = n), X = X, row.names = phy$tip.label)
      sim <- binaryPGLMM.sim(Y ~ 1 + X, phy = phy, data = sim.dat, s2 = s2,
                                       B = matrix(c(0,B1), nrow = 2, ncol = 1), nrep = 1)
      sim.dat$Y <- sim$Y

      z <- binaryPGLMM(Y ~ 1 + X, phy = phy, data = sim.dat)

      meanEsts.n[rep, ] <- c(z$B[2], z$s2, z$P.H0.s2, z$convergeflag == "converged")
  }
converged <- meanEsts.n[,4]
meanEsts <- rbind(meanEsts,
                  c(n, mean(meanEsts.n[converged==1,1]),
                            mean(meanEsts.n[converged==1,2]),
                            mean(meanEsts.n[converged==1, 3] < 0.05),
                            mean(converged)))
}
meanEsts

# Results output for B1 = 0.5, s2 = 0.4; n-Inf gives the values used to
# simulate the data
#    n       B1        s2      Pr.s2 propconverged
# 1 Inf 1.000000 0.4000000 1.00000000         1.000
# 2 160 1.012719 0.4479946 0.36153072         0.993
# 3  80 1.030876 0.5992027 0.24623116         0.995
# 4  40 1.110201 0.7425203 0.13373860         0.987
# 5  20 1.249886 0.8774708 0.05727377         0.873


## Examine type I errors for estimates of B0 and s2 from binaryPGLMM()
# with simulated data. Note that this will take a while.

Reps = 1000

s2 <- 0
B0 <- 0
B1 <- 0

H0.tests <- data.frame(n = Inf, B0 = B0, s2 = s2, Pr.B0 = .05,
                       Pr.s2 = .05, propconverged = 1)
for (n in c(160, 80, 40, 20)) {

  ests.n <- data.frame(B1 = 0, s2 = 0, Pr.B0 = 0, Pr.s2 = 0, convergefailure = 0)
  for (rep in 1:Reps) {
    phy <- compute.brlen(rtree(n = n), method = "Grafen", power = 1)
    X <- rTraitCont(phy, model = "BM", sigma = 1)
    X <- (X - mean(X))/var(X)

    sim.dat <- data.frame(Y = array(0, dim = n), X = X, row.names = phy$tip.label)
    sim <- binaryPGLMM.sim(Y ~ 1, phy = phy, data = sim.dat, s2 = s2,
                           B = matrix(B0, nrow = 1, ncol = 1), nrep = 1)
    sim.dat$Y <- sim$Y

    z <- binaryPGLMM(Y ~ 1, phy = phy, data = sim.dat)

    ests.n[rep, ] <- c(z$B[1], z$s2, z$B.pvalue, z$P.H0.s2, z$convergeflag == "converged")
  }

converged <- ests.n[,5]
H0.tests <- rbind(H0.tests,
                  c(n, mean(ests.n[converged==1,1]),
                    mean(ests.n[converged==1,2]),
                    mean(ests.n[converged==1, 3] < 0.05),
                    mean(ests.n[converged==1, 4] < 0.05),
                    mean(converged)))
}
H0.tests

# Results for type I errors for B0 = 0 and s2 = 0; n-Inf gives the values
# used to simulate the data. These results show that binaryPGLMM() tends to
# have lower-than-nominal p-values; fewer than 0.05 of the simulated
# data sets have H0:B0=0 and H0:s2=0 rejected at the alpha=0.05 level.
#     n            B0         s2      Pr.B0      Pr.s2 propconverged
# 1 Inf  0.0000000000 0.00000000 0.05000000 0.05000000         1.000
# 2 160 -0.0009350357 0.07273163 0.02802803 0.04804805         0.999
# 3  80 -0.0085831477 0.12205876 0.04004004 0.03403403         0.999
# 4  40  0.0019303847 0.25486307 0.02206620 0.03711133         0.997
# 5  20  0.0181394905 0.45949266 0.02811245 0.03313253         0.996

## End(Not run)

Description

This function binds together two phylogenetic trees to give a single object of class "phylo".

Usage

bind.tree(x, y, where = "root", position = 0, interactive = FALSE)
x + y

Arguments

Details

The argument x can be seen as the receptor tree, whereas y is the donor tree. The root of y is then grafted on a location of x specified by where and, possibly, position. If y has a root edge, this is added as in internal branch in the resulting tree.

x + y is a shortcut for:

    bind.tree(x, y, position = if (is.null(x$root.edge)) 0 else
    x$root.edge)
  

If only one of the trees has no branch length, the branch lengths of the other one are ignored with a warning.

If one (or both) of the trees has no branch length, it is possible to specify a value of 'position' to graft 'y' below the node of 'x' specified by 'where'. In this case, the exact value of 'position' is not important as long as it is greater than zero. The new node will be multichotomous if 'y' has no root edge. This can be solved by giving an arbitrary root edge to 'y' beforehand (e.g., y$root.edge <- 1): it will be deleted during the binding operation.

Value

an object of class "phylo".

Author(s)

Emmanuel Paradis

See Also

drop.tip, root

Examples

### binds the two clades of bird orders
treefile1 <- tempfile("tree", fileext = ".tre")
treefile2 <- tempfile("tree", fileext = ".tre")
cat("((Struthioniformes:21.8,Tinamiformes:21.8):4.1,",
    "((Craciformes:21.6,Galliformes:21.6):1.3,Anseriformes:22.9):3.0):2.1;",
    file = treefile1, sep = "\n")
cat("(Turniciformes:27.0,(Piciformes:26.3,((Galbuliformes:24.4,",
    "((Bucerotiformes:20.8,Upupiformes:20.8):2.6,",
    "(Trogoniformes:22.1,Coraciiformes:22.1):1.3):1.0):0.6,",
    "(Coliiformes:24.5,(Cuculiformes:23.7,(Psittaciformes:23.1,",
    "(((Apodiformes:21.3,Trochiliformes:21.3):0.6,",
    "(Musophagiformes:20.4,Strigiformes:20.4):1.5):0.6,",
    "((Columbiformes:20.8,(Gruiformes:20.1,Ciconiiformes:20.1):0.7):0.8,",
    "Passeriformes:21.6):0.9):0.6):0.6):0.8):0.5):1.3):0.7):1.0;",
    file = treefile2, sep = "\n")
tree.bird1 <- read.tree(treefile1)
tree.bird2 <- read.tree(treefile2)
unlink(c(treefile1, treefile2)) # clean-up
(birds <- tree.bird1 + tree.bird2)
layout(matrix(c(1, 2, 3, 3), 2, 2))
plot(tree.bird1)
plot(tree.bird2)
plot(birds)

### examples with random trees
x <- rtree(4, tip.label = LETTERS[1:4])
y <- rtree(4, tip.label = LETTERS[5:8])
x <- makeNodeLabel(x, prefix = "x_")
y <- makeNodeLabel(y, prefix = "y_")
x$root.edge <- y$root.edge <- .2

z <- bind.tree(x, y, po=.2)
plot(y, show.node.label = TRUE, font = 1, root.edge = TRUE)
title("y")
plot(x, show.node.label = TRUE, font = 1, root.edge = TRUE)
title("x")
plot(z, show.node.label = TRUE, font = 1, root.edge = TRUE)
title("z <- bind.tree(x, y, po=.2)")

## make sure the terminal branch length is long enough:
x$edge.length[x$edge[, 2] == 2] <- 0.2

z <- bind.tree(x, y, 2, .1)
plot(y, show.node.label = TRUE, font = 1, root.edge = TRUE)
title("y")
plot(x, show.node.label = TRUE, font = 1, root.edge = TRUE)
title("x")
plot(z, show.node.label = TRUE, font = 1, root.edge = TRUE)
title("z <- bind.tree(x, y, 2, .1)")

x <- rtree(50)
y <- rtree(50)
x$root.edge <- y$root.edge <- .2
z <- x + y
plot(y, show.tip.label = FALSE, root.edge = TRUE); axisPhylo()
title("y")
plot(x, show.tip.label = FALSE, root.edge = TRUE); axisPhylo()
title("x")
plot(z, show.tip.label = FALSE, root.edge = TRUE); axisPhylo()
title("z <- x + y")
layout(1)

Description

This function performs the BIONJ algorithm of Gascuel (1997).

Usage

bionj(X)

Arguments

Value

an object of class "phylo".

Author(s)

original C code by Hoa Sien Cuong and Olivier Gascuel; adapted and ported to R by Vincent Lefort [email protected]

References

Gascuel, O. (1997) BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data. Molecular Biology and Evolution, 14:, 685–695.

See Also

nj, fastme, mvr, bionjs, SDM, dist.dna

Examples

### From Saitou and Nei (1987, Table 1):
x <- c(7, 8, 11, 13, 16, 13, 17, 5, 8, 10, 13,
       10, 14, 5, 7, 10, 7, 11, 8, 11, 8, 12,
       5, 6, 10, 9, 13, 8)
M <- matrix(0, 8, 8)
M[lower.tri(M)] <- x
M <- t(M)
M[lower.tri(M)] <- x
dimnames(M) <- list(1:8, 1:8)
tr <- bionj(M)
plot(tr, "u")
### a less theoretical example
data(woodmouse)
trw <- bionj(dist.dna(woodmouse))
plot(trw)

Description

This data set describes the phylogenetic relationships of the families of birds as reported by Sibley and Ahlquist (1990). Sibley and Ahlquist inferred this phylogeny from an extensive number of DNA/DNA hybridization experiments. The “tapestry” reported by these two authors (more than 1000 species out of the ca. 9000 extant bird species) generated a lot of debates.

The present tree is based on the relationships among families. A few families were not included in the figures in Sibley and Ahlquist, and thus are not included here as well. The branch lengths were calculated from the values of $\Delta T_{50}H$ as found in Sibley and Ahlquist (1990, figs. 354, 355, 356, and 369).

Usage

data(bird.families)

Format

The data are stored as an object of class "phylo" which structure is described in the help page of the function read.tree.

Source

Sibley, C. G. and Ahlquist, J. E. (1990) Phylogeny and classification of birds: a study in molecular evolution. New Haven: Yale University Press.

See Also

read.tree, bird.orders

Examples

data(bird.families)
op <- par(cex = 0.3)
plot(bird.families)
par(op)

Description

This data set describes the phylogenetic relationships of the orders of birds as reported by Sibley and Ahlquist (1990). Sibley and Ahlquist inferred this phylogeny from an extensive number of DNA/DNA hybridization experiments. The “tapestry” reported by these two authors (more than 1000 species out of the ca. 9000 extant bird species) generated a lot of debates.

The present tree is based on the relationships among orders. The branch lengths were calculated from the values of $\Delta T_{50}H$ as found in Sibley and Ahlquist (1990, fig. 353).

Usage

data(bird.orders)

Format

The data are stored as an object of class "phylo" which structure is described in the help page of the function read.tree.

Source

Sibley, C. G. and Ahlquist, J. E. (1990) Phylogeny and classification of birds: a study in molecular evolution. New Haven: Yale University Press.

See Also

read.tree, bird.families

Examples

data(bird.orders)
plot(bird.orders)

Description

This function fits by maximum likelihood a birth-death model to the branching times computed from a phylogenetic tree using the method of Nee et al. (1994).

Usage

birthdeath(phy)
## S3 method for class 'birthdeath'
print(x, ...)

Arguments

Details

Nee et al. (1994) used a re-parametrization of the birth-death model studied by Kendall (1948) so that the likelihood has to be maximized over d/b and b - d, where b is the birth rate, and d the death rate. This is the approach used by the present function.

This function computes the standard-errors of the estimated parameters using a normal approximations of the maximum likelihood estimates: this is likely to be inaccurate because of asymmetries of the likelihood function (Nee et al. 1995). In addition, 95 intervals of both parameters are computed using profile likelihood: they are particularly useful if the estimate of d/b is at the boundary of the parameter space (i.e. 0, which is often the case).

Note that the function does not check that the tree is effectively ultrametric, so if it is not, the returned result may not be meaningful.

Value

An object of class "birthdeath" which is a list with the following components:

Author(s)

Emmanuel Paradis

References

Kendall, D. G. (1948) On the generalized “birth-and-death” process. Annals of Mathematical Statistics, 19, 1–15.

Nee, S., May, R. M. and Harvey, P. H. (1994) The reconstructed evolutionary process. Philosophical Transactions of the Royal Society of London. Series B. Biological Sciences, 344, 305–311.

Nee, S., Holmes, E. C., May, R. M. and Harvey, P. H. (1995) Estimating extinctions from molecular phylogenies. in Extinction Rates, eds. Lawton, J. H. and May, R. M., pp. 164–182, Oxford University Press.

See Also

branching.times, diversi.gof, diversi.time, ltt.plot, yule, bd.ext, yule.cov, bd.time

Description

These functions analyse bipartitions found in a series of trees.

prop.part counts the number of bipartitions found in a series of trees given as .... If a single tree is passed, the returned object is a list of vectors with the tips descending from each node (i.e., clade compositions indexed by node number).

prop.clades counts the number of times the bipartitions present in phy are present in a series of trees given as ... or in the list previously computed and given with part.

boot.phylo performs a bootstrap analysis.

Usage

boot.phylo(phy, x, FUN, B = 100, block = 1,
           trees = FALSE, quiet = FALSE,
           rooted = is.rooted(phy), jumble = TRUE,
            mc.cores = 1)
prop.part(..., check.labels = TRUE)
prop.clades(phy, ..., part = NULL, rooted = FALSE)
## S3 method for class 'prop.part'
print(x, ...)
## S3 method for class 'prop.part'
summary(object, ...)
## S3 method for class 'prop.part'
plot(x, barcol = "blue", leftmar = 4, col = "red", ...)

Arguments

Details

The argument FUN in boot.phylo must be the function used to estimate the tree from the original data matrix. Thus, if the tree was estimated with neighbor-joining (see nj), one maybe wants something like FUN = function(xx) nj(dist.dna(xx)).

block in boot.phylo specifies the number of columns to be resampled altogether. For instance, if one wants to resample at the codon-level, then block = 3 must be used.

Using check.labels = FALSE in prop.part decreases computing times. This requires that (i) all trees have the same tip labels, and (ii) these labels are ordered similarly in all trees (in other words, the element tip.label are identical in all trees).

The plot function represents a contingency table of the different partitions (on the x-axis) in the lower panel, and their observed numbers in the upper panel. Any further arguments (...) are used to change the aspects of the points in the lower panel: these may be pch, col, bg, cex, etc. This function works only if there is an attribute labels in the object.

The print method displays the partitions and their numbers. The summary method extracts the numbers only.

Value

prop.part returns an object of class "prop.part" which is a list with an attribute "number". The elements of this list are the observed clades, and the attribute their respective numbers. If the default check.labels = FALSE is used, an attribute "labels" is added, and the vectors of the returned object contains the indices of these labels instead of the labels themselves.

prop.clades and boot.phylo return a numeric vector which ith element is the number associated to the ith node of phy. If trees = TRUE, boot.phylo returns a list whose first element (named "BP") is like before, and the second element ("trees") is a list with the bootstraped trees.

summary returns a numeric vector.

Note

prop.clades calls internally prop.part with the option check.labels = TRUE, which may be very slow. If the trees passed as ... fulfills conditions (i) and (ii) above, then it might be faster to first call, e.g., pp <- prop.part(...), then use the option part: prop.clades(phy, part = pp).

Since ape 3.5, prop.clades should return sensible results for all values of rooted: if FALSE, the numbers of bipartitions (or splits); if TRUE, the number of clades (of hopefully rooted trees).

Author(s)

Emmanuel Paradis

References

Efron, B., Halloran, E. and Holmes, S. (1996) Bootstrap confidence levels for phylogenetic trees. Proceedings of the National Academy of Sciences USA, 93, 13429–13434.

Felsenstein, J. (1985) Confidence limits on phylogenies: an approach using the bootstrap. Evolution, 39, 783–791.

See Also

as.bitsplits, dist.topo, consensus, nodelabels

Examples

data(woodmouse)
f <- function(x) nj(dist.dna(x))
tr <- f(woodmouse)
### Are bootstrap values stable?
for (i in 1:5)
  print(boot.phylo(tr, woodmouse, f, quiet = TRUE))
### How many partitions in 100 random trees of 10 labels?...
TR <- rmtree(100, 10)
pp10 <- prop.part(TR)
length(pp10)
### ... and in 100 random trees of 20 labels?
TR <- rmtree(100, 20)
pp20 <- prop.part(TR)
length(pp20)
plot(pp10, pch = "x", col = 2)
plot(pp20, pch = "x", col = 2)

set.seed(2)
tr <- rtree(10) # rooted
## the following used to return a wrong result with ape <= 3.4:
prop.clades(tr, tr)
prop.clades(tr, tr, rooted = TRUE)
tr <- rtree(10, rooted = FALSE)
prop.clades(tr, tr) # correct

### an illustration of the use of prop.clades with bootstrap trees:

fun <- function(x) as.phylo(hclust(dist.dna(x), "average")) # upgma() in phangorn
tree <- fun(woodmouse)
## get 100 bootstrap trees:
bstrees <- boot.phylo(tree, woodmouse, fun, trees = TRUE)$trees
## get proportions of each clade:
clad <- prop.clades(tree, bstrees, rooted = TRUE)
## get proportions of each bipartition:
boot <- prop.clades(tree, bstrees)
layout(1)
par(mar = rep(2, 4))
plot(tree, main = "Bipartition vs. Clade Support Values")
drawSupportOnEdges(boot)
nodelabels(clad)
legend("bottomleft", legend = c("Bipartitions", "Clades"), pch = 22,
       pt.bg = c("green", "lightblue"), pt.cex = 2.5)

## Not run: 
## an example of double bootstrap:
nrep1 <- 100
nrep2 <- 100
p <- ncol(woodmouse)
DB <- 0

for (b in 1:nrep1) {
    X <- woodmouse[, sample(p, p, TRUE)]
    DB <- DB + boot.phylo(tr, X, f, nrep2, quiet = TRUE)
}
DB
## to compare with:
boot.phylo(tr, woodmouse, f, 1e4)

## End(Not run)

Description

This function computes the branching times of a phylogenetic tree, that is the distance from each node to the tips, under the assumption that the tree is ultrametric. Note that the function does not check that the tree is effectively ultrametric, so if it is not, the returned result may not be meaningful.

Usage

branching.times(phy)

Arguments

Value

a numeric vector with the branching times. If the phylogeny phy has an element node.label, this is used as names for the returned vector; otherwise the numbers (of mode character) of the matrix edge of phy are used as names.

Author(s)

Emmanuel Paradis

See Also

is.ultrametric

Description

These functions help to build lists of trees of class "multiPhylo".

Usage

## S3 method for class 'phylo'
c(..., recursive = TRUE)
## S3 method for class 'multiPhylo'
c(..., recursive = TRUE)
.compressTipLabel(x, ref = NULL)
.uncompressTipLabel(x)

Arguments

Details

These c methods check all the arguments, and return by default a list of single trees unless some objects are not trees or lists of trees, in which case recursive is switched to FALSE and a warning message is given. If recursive = FALSE, the objects are simply concatenated into a list. Before ape 4.0, recursive was always set to FALSE.

.compressTipLabel transforms an object of class "multiPhylo" by checking that all trees have the same tip labels and renumbering the tips in the edge matrix so that the tip numbers are also the same taking the first tree as the reference (duplicated labels are not allowed). The returned object has a unique vector of tip labels (attr(x, "TipLabel")).

.uncompressTipLabel does the reverse operation.

Value

An object of class "multiPhylo".

Author(s)

Emmanuel Paradis

See Also

summary.phylo, multiphylo

Examples

x <- c(rtree(4), rtree(2))
x
y <- c(rtree(4), rtree(4))
z <- c(x, y)
z
print(z, TRUE)
try(.compressTipLabel(x)) # error
a <- .compressTipLabel(y)
.uncompressTipLabel(a) # back to y
## eventually compare str(a) and str(y)

Description

Function CADM.global compute and test the coefficient of concordance among several distance matrices through a permutation test.

Function CADM.post carries out a posteriori permutation tests of the contributions of individual distance matrices to the overall concordance of the group.

Use in phylogenetic analysis: to identify congruence among distance matrices (D) representing different genes or different types of data. Congruent D matrices correspond to data tables that can be used together in a combined phylogenetic or other type of multivariate analysis.

Usage

CADM.global(Dmat, nmat, n, nperm=99, make.sym=TRUE, weights=NULL,
            silent=FALSE)
CADM.post  (Dmat, nmat, n, nperm=99, make.sym=TRUE, weights=NULL,
             mult="holm", mantel=FALSE, silent=FALSE)

Arguments

Details

Dmat must contain two or more distance matrices, listed one after the other, all of the same size, and corresponding to the same objects in the same order. Raw data tables can be transformed into distance matrices before comparison with other such distance matrices, or with data that have been obtained as distance matrices, e.g. serological or DNA hybridization data. The distances will be transformed to ranks before computation of the coefficient of concordance and other statistics.

CADM.global tests the global null hypothesis that all matrices are incongruent. If the global null is rejected, function CADM.post can be used to identify the concordant (H0 rejected) and discordant matrices (H0 not rejected) in the group. If a distance matrix has a negative value for the Mantel.mean statistic, that matrix clearly does not belong to the group. Remove that matrix (if there are more than one, remove first the matrix that has the most strongly negative value for Mantel.mean) and run the analysis again.

The corrections used for multiple testing are applied to the list of P-values (P) produced in the a posteriori tests; they take into account the number of tests (k) carried out simulatenously (number of matrices, parameter nmat).

The Holm correction is computed after ordering the P-values in a list with the smallest value to the left. Compute adjusted P-values as:

$P_{corr} = (k-i+1)*P$

where i is the position in the ordered list. Final step: from left to right, if an adjusted $P_{corr}$ in the ordered list is smaller than the one occurring at its left, make the smallest one equal to the largest one.

The Sidak correction is:

$P_{corr} = 1 - (1 - P)^k$

The Bonferonni correction is:

$P_{corr} = k*P$

Value

CADM.global produces a small table containing the W, Chi2, and Prob.perm statistics described in the following list. CADM.post produces a table stored in element A_posteriori_tests, containing Mantel.mean, Prob, and Corrected.prob statistics in rows; the columns correspond to the k distance matrices under study, labeled Dmat.1 to Dmat.k. If parameter mantel is TRUE, tables of Mantel statistics and P-values are computed among the matrices.

Author(s)

Pierre Legendre, Universite de Montreal

References

Campbell, V., Legendre, P. and Lapointe, F.-J. (2009) Assessing congruence among ultrametric distance matrices. Journal of Classification, 26, 103–117.

Campbell, V., Legendre, P. and Lapointe, F.-J. (2011) The performance of the Congruence Among Distance Matrices (CADM) test in phylogenetic analysis. BMC Evolutionary Biology, 11, 64.

Friedman, M. (1937) The use of ranks to avoid the assumption of normality implicit in the analysis of variance. Journal of the American Statistical Association, 32, 675–701.

Kendall, M. G. and Babington Smith, B. (1939) The problem of m rankings. Annals of Mathematical Statistics, 10, 275–287.

Lapointe, F.-J., Kirsch, J. A. W. and Hutcheon, J. M. (1999) Total evidence, consensus, and bat phylogeny: a distance-based approach. Molecular Phylogenetics and Evolution, 11, 55–66.

Legendre, P. (2010) Coefficient of concordance. Pp. 164-169 in: Encyclopedia of Research Design, Vol. 1. N. J. Salkind, ed. SAGE Publications, Inc., Los Angeles.

Legendre, P. and Lapointe, F.-J. (2004) Assessing congruence among distance matrices: single malt Scotch whiskies revisited. Australian and New Zealand Journal of Statistics, 46, 615–629.

Legendre, P. and Lapointe, F.-J. (2005) Congruence entre matrices de distance. P. 178-181 in: Makarenkov, V., G. Cucumel et F.-J. Lapointe [eds] Comptes rendus des 12emes Rencontres de la Societe Francophone de Classification, Montreal, 30 mai - 1er juin 2005.

Siegel, S. and Castellan, N. J., Jr. (1988) Nonparametric statistics for the behavioral sciences. 2nd edition. New York: McGraw-Hill.

Examples

# Examples 1 and 2: 5 genetic distance matrices computed from simulated DNA
# sequences representing 50 taxa having evolved along additive trees with
# identical evolutionary parameters (GTR+ Gamma + I). Distance matrices were
# computed from the DNA sequence matrices using a p distance corrected with the
# same parameters as those used to simulate the DNA sequences. See Campbell et
# al. (2009) for details.

# Example 1: five independent additive trees. Data provided by V. Campbell.

data(mat5Mrand)
res.global <- CADM.global(mat5Mrand, 5, 50)

# Example 2: three partly similar trees, two independent trees.
# Data provided by V. Campbell.

data(mat5M3ID)
res.global <- CADM.global(mat5M3ID, 5, 50)
res.post   <- CADM.post(mat5M3ID, 5, 50, mantel=TRUE)

# Example 3: three matrices respectively representing Serological
# (asymmetric), DNA hybridization (asymmetric) and Anatomical (symmetric)
# distances among 9 families. Data from Lapointe et al. (1999).

data(mat3)
res.global <- CADM.global(mat3, 3, 9, nperm=999)
res.post   <- CADM.post(mat3, 3, 9, nperm=999, mantel=TRUE)

# Example 4, showing how to bind two D matrices (cophenetic matrices
# in this example) into a file using rbind(), then run the global test.

a <- rtree(5)
b <- rtree(5)
A <- cophenetic(a)
B <- cophenetic(b)
x <- rownames(A)
B <- B[x, x]
M <- rbind(A, B)
CADM.global(M, 2, 5)

Description

Dataset adapted from Gittleman (1986), including 2 morphological variables (body and brain sizes), 8 life history traits variables and 4 taxonomic variables.

Usage

data(carnivora)

Format

A data frame with 112 observations on 17 variables.

Source

Gittleman, J. L. (1986) Carnivore life history patterns: allometric, phylogenetic and ecological associations. American Naturalist, 127: 744–771.

Examples

data(carnivora)
## Fig. 1 in Gittleman (1986):
plot(carnivora$BW ~ carnivora$FW, pch = (1:8)[carnivora$Family], log = "xy",
     xlab = "Female body weight (kg)", ylab = "Birth weigth (g)",
     ylim = c(1, 2000))
legend("bottomright", legend = levels(carnivora$Family), pch = 1:8)
plot(carnivora$BW ~ carnivora$FB, pch = (1:8)[carnivora$Family], log = "xy",
     xlab = "Female brain weight (g)", ylab = "Birth weigth (g)",
     ylim = c(1, 2000))
legend("bottomright", legend = levels(carnivora$Family), pch = 1:8)

Description

This function performs a series of diagnostics on a DNA alignement.

Usage

checkAlignment(x, check.gaps = TRUE, plot = TRUE, what = 1:4)

Arguments

Details

This function prints on the console a series of diagnostics on the set a aligned DNA sequences. If alignment gaps are present, their width distribution is analysed, as well as the width of contiguous base segments. The pattern of nucleotide diversity on each site is also analysed, and a relevant table is printed.

If plot = TRUE, four plots are done: an image of the alignement, the distribution of gap widths (if present), the Shannon index of nucleotide diversity along the sequence, and the number of observed bases along the sequence.

If the sequences contain many gaps, it might be better to set check.gaps = FALSE to skip the analysis of contiguous segments.

Value

NULL

Author(s)

Emmanuel Paradis

See Also

alview, image.DNAbin, all.equal.DNAbin

Examples

data(woodmouse)
checkAlignment(woodmouse)
layout(1)

Description

Checking and correcting character strings, particularly before writing a Newick tree.

Usage

checkLabel(x)

Arguments

Details

This function deletes the leading and trailing spaces (including tabulations, new lines, and left or right parentheses at the beginning or end of the strings), substitutes the spaces inside the strings by underscores, and substitutes commas, colons, semicolons, and parentheses inside the strings by dashes.

Value

a vector of mode character.

Author(s)

Emmanuel Paradis

See Also

makeLabel, makeNodeLabel, mixedFontLabel, stripLabel, updateLabel

Examples

checkLabel(" Homo sapiens\t(Primates; World)   ")

Description

This function takes as single argument an object (phy), checks its elements, and prints a diagnostic. All problems are printed with a label: FATAL (will likely cause an error or a crash) or MODERATE (may cause some problems).

This function is mainly intended for developers creating "phylo" objects from scratch.

Usage

checkValidPhylo(phy)

Arguments

Value

NULL.

Author(s)

Emmanuel Paradis

Examples

tr <- rtree(3)
checkValidPhylo(tr)
tr$edge[1] <- 0
checkValidPhylo(tr)

Description

This function calculates the number of cherries (see definition below) on a phylogenetic tree, and tests the null hypotheses whether this number agrees with those predicted from two null models of trees (the Yule model, and the uniform model).

Usage

cherry(phy)

Arguments

Details

A cherry is a pair of adjacent tips on a tree. The tree can be either rooted or unrooted, but the present function considers only rooted trees. The probability distribution function of the number of cherries on a tree depends on the speciation/extinction model that generated the tree.

McKenzie and Steel (2000) derived the probability distribution function of the number of cherries for two models: the Yule model and the uniform model. Broadly, in the Yule model, each extant species is equally likely to split into two daughter-species; in the uniform model, a branch is added to tree on any of the already existing branches with a uniform probability.

The probabilities are computed using recursive formulae; however, for both models, the probability density function converges to a normal law with increasing number of tips in the tree. The function uses these normal approximations for a number of tips greater than or equal to 20.

Value

A NULL value is returned, the results are simply printed.

Author(s)

Emmanuel Paradis

References

McKenzie, A. and Steel, M. (2000) Distributions of cherries for two models of trees. Mathematical Biosciences, 164, 81–92.

See Also

gammaStat

Description

This phylogeny of bats (Mammalia: Chiroptera) is a supertree (i.e. a composite phylogeny constructed from several sources; see source for details).

Usage

data(chiroptera)

Format

The data are stored in RData (binary) format.

Source

Jones, K. E., Purvis, A., MacLarnon, A., Bininda-Emonds, O. R. P. and Simmons, N. B. (2002) A phylogenetic supertree of the bats (Mammalia: Chiroptera). Biological Reviews of the Cambridge Philosophical Society, 77, 223–259.

See Also

read.nexus, zoom

Examples

data(chiroptera)
str(chiroptera)
op <- par(cex = 0.3)
plot(chiroptera, type = "c")
par(op)

Description

This function estimates the node ages of a tree using the mean path lengths method of Britton et al. (2002). The branch lengths of the input tree are interpreted as (mean) numbers of substitutions.

Usage

chronoMPL(phy, se = TRUE, test = TRUE)

Arguments

Details

The mean path lengths (MPL) method estimates the age of a node with the mean of the distances from this node to all tips descending from it. Under the assumption of a molecular clock, standard-errors of the estimates node ages can be computed (Britton et al. 2002).

The tests performed if test = TRUE is a comparison of the MPL of the two subtrees originating from a node; the null hypothesis is that the rate of substitution was the same in both subtrees (Britton et al. 2002). The test statistic follows, under the null hypothesis, a standard normal distribution. The returned P-value is the probability of observing a greater absolute value (i.e., a two-sided test). No correction for multiple testing is applied: this is left to the user.

Absolute dating can be done by multiplying the edge lengths found by calibrating one node age.

Value

an object of class "phylo" with branch lengths as estimated by the function. There are, by default, two attributes:

Note

The present version requires a dichotomous tree.

Author(s)

Emmanuel Paradis

References

Britton, T., Oxelman, B., Vinnersten, A. and Bremer, K. (2002) Phylogenetic dating with confidence intervals using mean path lengths. Molecular Phylogenetics and Evolution, 24, 58–65.

See Also

chronopl

Examples

tr <- rtree(10)
tr$edge.length <- 5*tr$edge.length
chr <- chronoMPL(tr)
layout(matrix(1:4, 2, 2, byrow = TRUE))
plot(tr)
title("The original tree")
plot(chr)
axisPhylo()
title("The dated MPL tree")
plot(chr)
nodelabels(round(attr(chr, "stderr"), 3))
title("The standard-errors")
plot(tr)
nodelabels(round(attr(chr, "Pval"), 3))
title("The tests")
layout(1)

Description

This function estimates the node ages of a tree using a semi-parametric method based on penalized likelihood (Sanderson 2002). The branch lengths of the input tree are interpreted as mean numbers of substitutions (i.e., per site).

Usage

chronopl(phy, lambda, age.min = 1, age.max = NULL,
         node = "root", S = 1, tol = 1e-8,
         CV = FALSE, eval.max = 500, iter.max = 500, ...)

Arguments

Details

The idea of this method is to use a trade-off between a parametric formulation where each branch has its own rate, and a nonparametric term where changes in rates are minimized between contiguous branches. A smoothing parameter (lambda) controls this trade-off. If lambda = 0, then the parametric component dominates and rates vary as much as possible among branches, whereas for increasing values of lambda, the variation are smoother to tend to a clock-like model (same rate for all branches).

lambda must be given. The known ages are given in age.min, and the correponding node numbers in node. These two arguments must obviously be of the same length. By default, an age of 1 is assumed for the root, and the ages of the other nodes are estimated.

If age.max = NULL (the default), it is assumed that age.min gives exactly known ages. Otherwise, age.max and age.min must be of the same length and give the intervals for each node. Some node may be known exactly while the others are known within some bounds: the values will be identical in both arguments for the former (e.g., age.min = c(10, 5), age.max = c(10, 6), node = c(15, 18) means that the age of node 15 is 10 units of time, and the age of node 18 is between 5 and 6).

If two nodes are linked (i.e., one is the ancestor of the other) and have the same values of age.min and age.max (say, 10 and 15) this will result in an error because the medians of these values are used as initial times (here 12.5) giving initial branch length(s) equal to zero. The easiest way to solve this is to change slightly the given values, for instance use age.max = 14.9 for the youngest node, or age.max = 15.1 for the oldest one (or similarly for age.min).

The input tree may have multichotomies. If some internal branches are of zero-length, they are collapsed (with a warning), and the returned tree will have less nodes than the input one. The presence of zero-lengthed terminal branches of results in an error since it makes little sense to have zero-rate branches.

The cross-validation used here is different from the one proposed by Sanderson (2002). Here, each tip is dropped successively and the analysis is repeated with the reduced tree: the estimated dates for the remaining nodes are compared with the estimates from the full data. For the $i$th tip the following is calculated:

$\sum_{j=1}^{n-2}{\frac{(t_j - t_j^{-i})^2}{t_j}}$

,

where $t_j$ is the estimated date for the $j$th node with the full phylogeny, $t_j^{-i}$ is the estimated date for the $j$th node after removing tip $i$ from the tree, and $n$ is the number of tips.

The present version uses the nlminb to optimise the penalized likelihood function: see its help page for details on parameters controlling the optimisation procedure.

Value

an object of class "phylo" with branch lengths as estimated by the function. There are three or four further attributes:

Note

The new function chronos replaces the present one which is no more maintained.

Author(s)

Emmanuel Paradis

References

Sanderson, M. J. (2002) Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach. Molecular Biology and Evolution, 19, 101–109.

See Also

chronos, chronoMPL

Description

chronos is the main function fitting a chronogram to a phylogenetic tree whose branch lengths are in number of substitution per sites.

makeChronosCalib is a tool to prepare data frames with the calibration points of the phylogenetic tree.

chronos.control creates a list of parameters to be passed to chronos.

Usage

chronos(phy, lambda = 1, model = "correlated", quiet = FALSE,
        calibration = makeChronosCalib(phy),
        control = chronos.control())
## S3 method for class 'chronos'
print(x, ...)
makeChronosCalib(phy, node = "root", age.min = 1,
   age.max = age.min, interactive = FALSE, soft.bounds = FALSE)
chronos.control(...)

Arguments

Details

chronos replaces chronopl but with a different interface and some extensions (see References).

The known dates (argument calibration) must be given in a data frame with the following column names: node, age.min, age.max, and soft.bounds (the last one is yet unused). For each row, these are, respectively: the number of the node in the “phylo” coding standard, the minimum age for this node, the maximum age, and a logical value specifying whether the bounds are soft. If age.min = age.max, this means that the age is exactly known. This data frame can be built with makeChronosCalib which returns by default a data frame with a single row giving age = 1 for the root. The data frame can be built interactively by clicking on the plotted tree.

The argument control allows one to change some parameters of the optimisation procedure. This must be a list with names. The available options with their default values are:

• tol = 1e-8: tolerance for the estimation of the substitution rates.

• iter.max = 1e4: the maximum number of iterations at each optimization step.

• eval.max = 1e4: the maximum number of function evaluations at each optimization step.

• nb.rate.cat = 10: the number of rate categories if model = "discrete" (set this parameter to 1 to fit a strict clock model).

• dual.iter.max = 20: the maximum number of alternative iterations between rates and dates.

• epsilon = 1e-6: the convergence diagnostic criterion.

Using model = "clock" is actually a short-cut to model = "discrete" and setting nb.rate.cat = 1 in the list passed to control.

The command chronos.control() returns a list with the default values of these parameters. They may be modified by passing them to this function, or directly in the list.

Value

chronos returns an object of class c("chronos", "phylo"). There is a print method for it. There are additional attributes which can be visualised with str or extracted with attr.

makeChronosCalib returns a data frame.

chronos.control returns a list.

Author(s)

Emmanuel Paradis, Santiago Claramunt, Guillaume Louvel

References

Kim, J. and Sanderson, M. J. (2008) Penalized likelihood phylogenetic inference: bridging the parsimony-likelihood gap. Systematic Biology, 57, 665–674.

Paradis, E. (2013) Molecular dating of phylogenies by likelihood methods: a comparison of models and a new information criterion. Molecular Phylogenetics and Evolution, 67, 436–444.

Sanderson, M. J. (2002) Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach. Molecular Biology and Evolution, 19, 101–109.

See Also

chronoMPL

Examples

library(ape)
tr <- rtree(10)
### the default is the correlated rate model:
chr <- chronos(tr)
### strict clock model:
ctrl <- chronos.control(nb.rate.cat = 1)
chr.clock <- chronos(tr, model = "discrete", control = ctrl)
### How different are the rates?
attr(chr, "rates")
attr(chr.clock, "rates")
## Not run: 
cal <- makeChronosCalib(tr, interactive = TRUE)
cal
### if you made mistakes, you can edit the data frame with:
### fix(cal)
chr <- chronos(tr, calibration = cal)

## End(Not run)

Description

These functions call their respective program from R to align a set of nucleotide sequences of class "DNAbin" or "AAbin". The application(s) must be installed seperately and it is highly recommended to do this so that the executables are in a directory located on the PATH of the system.

This version includes an experimental version of muscle5 which calls MUSCLE5 (see the link to the documentation in the References below); muscle still calls MUSCLE version 3. Note that the executable of MUSCLE5 is also named ‘muscle’ by the default compilation setting.

The functions efastats and letterconf require MUSCLE5.

Usage

clustal(x, y, guide.tree, pw.gapopen = 10, pw.gapext = 0.1,
        gapopen = 10, gapext = 0.2, exec = NULL, MoreArgs = "",
        quiet = TRUE, original.ordering = TRUE, file)
clustalomega(x, y, guide.tree, exec = NULL,MoreArgs = "",
              quiet = TRUE, original.ordering = TRUE, file)
muscle(x, y, guide.tree, exec, MoreArgs = "",
        quiet = TRUE, original.ordering = TRUE, file)
muscle5(x, exec = "muscle", MoreArgs = "", quiet = FALSE,
        file, super5 = FALSE, mc.cores = 1)
tcoffee(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE,
        original.ordering = TRUE)

efastats(X, exec = "muscle", quiet = FALSE)
letterconf(X, exec = "muscle")

Arguments

Details

It is highly recommended to install the executables properly so that they are in a directory located on the PATH (i.e., accessible from any other directory). Alternatively, the full path to the executable may be given (e.g., exec = "~/muscle/muscle"), or a (symbolic) link may be copied in the working directory. For Debian and its derivatives (e.g., Ubuntu), it is recommended to use the binaries distributed by Debian.

clustal tries to guess the name of the executable program depending on the operating system. Specifically, the followings are used: “clustalw” under Linux, “clustalw2” under MacOS, and “clustalw2.exe” under Windows. For clustalomega, “clustalo[.exe]” is the default on all systems (with no specific path).

When called without arguments (i.e., clustal(), ...), the function prints the options of the program which may be passed to MoreArgs.

Since ape 5.1, clustal, clustalomega, and muscle can align AA sequences as well as DNA sequences.

Value

an object of class "DNAbin" or "AAbin" with the aligned sequences.

efastats returns a data frame.

letterconf opens the default Web brower.

Author(s)

Emmanuel Paradis, Franz Krah

References

Chenna, R., Sugawara, H., Koike, T., Lopez, R., Gibson, T. J., Higgins, D. G. and Thompson, J. D. (2003) Multiple sequence alignment with the Clustal series of programs. Nucleic Acids Research 31, 3497–3500. http://www.clustal.org/

Edgar, R. C. (2004) MUSCLE: Multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Research, 32, 1792–1797. http://www.drive5.com/muscle/muscle_userguide3.8.html

Notredame, C., Higgins, D. and Heringa, J. (2000) T-Coffee: A novel method for multiple sequence alignments. Journal of Molecular Biology, 302, 205–217. https://tcoffee.org/

Sievers, F., Wilm, A., Dineen, D., Gibson, T. J., Karplus, K., Li, W., Lopez, R., McWilliam, H., Remmert, M., S\"oding, J., Thompson, J. D. and Higgins, D. G. (2011) Fast, scalable generation of high-quality protein multiple sequence alignments using Clustal Omega. Molecular Systems Biology, 7, 539. http://www.clustal.org/

https://drive5.com/muscle5/

See Also

image.DNAbin, del.gaps, all.equal.DNAbin, alex, alview, checkAlignment

Examples

## Not run: 
### display the options:
clustal()
clustalomega()
muscle()
tcoffee()

data(woodmouse)
### open gaps more easily:
clustal(woodmouse, pw.gapopen = 1, pw.gapext = 1)
### T-Coffee requires negative values (quite slow; muscle() is much faster):
tcoffee(woodmouse,  MoreArgs = "-gapopen=-10 -gapext=-2")

## End(Not run)

Description

This function extracts or generates information about coalescent intervals (number of lineages, interval lengths, interval count, total depth) from a phylogenetic tree or a list of internode distances. The input tree needs to be ultra-metric (i.e. clock-like).

Usage

coalescent.intervals(x)

Arguments

Value

An object of class "coalescentIntervals" with the following entries:

Author(s)

Korbinian Strimmer

See Also

branching.times, collapsed.intervals, read.tree.

Examples

data("hivtree.newick") # example tree in NH format
tree.hiv <- read.tree(text = hivtree.newick) # load tree
ci <- coalescent.intervals(tree.hiv) # from tree
ci
data("hivtree.table") # same tree, but in table format
ci <- coalescent.intervals(hivtree.table$size) # from vector of interval lengths
ci

Description

collapse.singles deletes the single nodes (i.e., with a single descendant) in a tree.

has.singles tests for the presence of single node(s) in a tree.

Usage

collapse.singles(tree, root.edge = FALSE)
has.singles(tree)

Arguments

Value

an object of class "phylo".

Author(s)

Emmanuel Paradis, Klaus Schliep

See Also

plot.phylo, read.tree

Examples

## a tree with 3 tips and 3 nodes:
e <- c(4L, 6L, 6L, 5L, 5L, 6L, 1L, 5L, 3L, 2L)
dim(e) <- c(5, 2)
tr <- structure(list(edge = e, tip.label = LETTERS[1:3], Nnode = 3L),
                class = "phylo")
tr
has.singles(tr)
## the following shows that node #4 (ie, the root) is a singleton
## and node #6 is the first bifurcating node
tr$edge
## A bifurcating tree has less nodes than it has tips:
## the following used to fail with ape 4.1 or lower:
plot(tr)
collapse.singles(tr) # only 2 nodes
## give branch lengths to use the 'root.edge' option:
tr$edge.length <- runif(5)
str(collapse.singles(tr, TRUE)) # has a 'root.edge'