Movatterモバイル変換

Version:	5.8-1
Date:	2024-12-10
Title:	Analyses of Phylogenetics and Evolution
Depends:	R (≥ 3.2.0)
Suggests:	gee, expm, igraph, phangorn, xml2
Imports:	nlme, lattice, graphics, methods, stats, utils, parallel, Rcpp(≥ 0.12.0), digest
LinkingTo:	Rcpp
ZipData:	no
Description:	Functions for reading, writing, plotting, and manipulating phylogenetic trees, analyses of comparative data in a phylogenetic framework, ancestral character analyses, analyses of diversification and macroevolution, computing distances from DNA sequences, reading and writing nucleotide sequences as well as importing from BioConductor, and several tools such as Mantel's test, generalized skyline plots, graphical exploration of phylogenetic data (alex, trex, kronoviz), estimation of absolute evolutionary rates and clock-like trees using mean path lengths and penalized likelihood, dating trees with non-contemporaneous sequences, translating DNA into AA sequences, and assessing sequence alignments. Phylogeny estimation can be done with the NJ, BIONJ, ME, MVR, SDM, and triangle methods, and several methods handling incomplete distance matrices (NJ*, BIONJ*, MVR*, and the corresponding triangle method). Some functions call external applications (PhyML, Clustal, T-Coffee, Muscle) whose results are returned into R.
License:	GPL-2 |GPL-3
URL:	https://github.com/emmanuelparadis/ape
BugReports:	https://github.com/emmanuelparadis/ape/issues
Encoding:	UTF-8
NeedsCompilation:	yes
Packaged:	2024-12-10 17:29:50 UTC; paradis
Author:	Emmanuel Paradis [aut, cre, cph], Simon Blomberg [aut, cph], Ben Bolker [aut, cph], Joseph Brown [aut, cph], Santiago Claramunt [aut, cph], Julien Claude [aut, cph], Hoa Sien Cuong [aut, cph], Richard Desper [aut, cph], Gilles Didier [aut, cph], Benoit Durand [aut, cph], Julien Dutheil [aut, cph], RJ Ewing [aut, cph], Olivier Gascuel [aut, cph], Thomas Guillerme [aut, cph], Christoph Heibl [aut, cph], Anthony Ives [aut, cph], Bradley Jones [aut, cph], Franz Krah [aut, cph], Daniel Lawson [aut, cph], Vincent Lefort [aut, cph], Pierre Legendre [aut, cph], Jim Lemon [aut, cph], Guillaume Louvel [aut, cph], Federico Marotta [aut, cph], Eric Marcon [aut, cph], Rosemary McCloskey [aut, cph], Johan Nylander [aut, cph], Rainer Opgen-Rhein [aut, cph], Andrei-Alin Popescu [aut, cph], Manuela Royer-Carenzi [aut, cph], Klaus Schliep [aut, cph], Korbinian Strimmer [aut, cph], Damien de Vienne [aut, cph]
Maintainer:	Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
Repository:	CRAN
Date/Publication:	2024-12-16 00:00:02 UTC

Analyses of Phylogenetics and Evolution

Description

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

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

More information onape can be found athttps://emmanuelparadis.github.io.

Author(s)

Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, RichardDesper, Benoit Durand, Julien Dutheil, Olivier Gascuel, ChristophHeibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon,Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Andrei-AlinPopescu, Klaus Schliep, Korbinian Strimmer, Damien de Vienne

Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>

References

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

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

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

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

Amino Acid Sequences

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 listswhich can be manipulated with the usual[ operator.

There is a conversion function to and from characters.

The functiondist.aa computes the number of AA differencesbetween each pair of sequences in a matrix; this can be scaled by thesequence length. See the functiondist.ml inphangorn for evolutionary distances with AA sequences.

The functionAAsubst returns the indices of the polymorphic sites(similar toseg.sites for DNA sequences; see examples below).

The two functionscbind.AAbin andrbind.AAbin work in thesame way than the similar methods for the class"DNAbin": seecbind.DNAbin for more explanations about their respectivebehaviours.

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)

Tree Estimation Based on an Improved Version of the NJ Algorithm

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 andported toR by Vincent Lefortvincent.lefort@lirmm.fr

References

Gascuel, O. (1997) BIONJ: an improved version of the NJ algorithmbased 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)] <- xM <- t(M)M[lower.tri(M)] <- xdimnames(M) <- list(1:8, 1:8)tr <- bionj(M)plot(tr, "u")### a less theoretical exampledata(woodmouse)trw <- bionj(dist.dna(woodmouse))plot(trw)

Congruence among distance matrices

Description

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

FunctionCADM.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, functionCADM.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 theMantel.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 forMantel.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, parameternmat).

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 adjustedP_{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 elementA_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 parametermantel 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 whiskiesrevisited.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)

Manipulate DNA Sequences in Bit-Level Format

Description

These functions help to manipulate DNA sequences coded in thebit-level coding scheme.

Usage

## S3 method for class 'DNAbin'print(x, printlen = 6, digits = 3, ...)## S3 method for class 'DNAbin'rbind(...)## S3 method for class 'DNAbin'cbind(..., check.names = TRUE, fill.with.gaps = FALSE,             quiet = FALSE)## S3 method for class 'DNAbin'x[i, j, drop = FALSE]## S3 method for class 'DNAbin'as.matrix(x, ...)## S3 method for class 'DNAbin'c(..., recursive = FALSE)## S3 method for class 'DNAbin'as.list(x, ...)## S3 method for class 'DNAbin'labels(object, ...)

Arguments

Details

These are all ‘methods’ of generic functions which are here applied toDNA sequences stored as objects of class"DNAbin". They areused in the same way than the standardR functions to manipulatevectors, matrices, and lists. Additionally, the operators[[and$ may be used to extract a vector from a list. Note thatthe default ofdrop is not the same than the generic operator:this is to avoid dropping rownames when selecting a single sequence.

These functions are provided to manipulate easily DNA sequences codedwith the bit-level coding scheme. The latter allows much fastercomparisons of sequences, as well as storing them in less memorycompared to the format used beforeape 1.10.

Forcbind, the default behaviour is to keep only individuals(as indicated by the rownames) for which there are no missing data. Iffill.with.gaps = TRUE, a ‘complete’ matrix is returned,enventually with insertion gaps as missing data. Ifcheck.names = TRUE (the default), the rownames of each matrix are checked, andthe rows are reordered if necessary (if some rownames are duplicated,an error is returned). Ifcheck.names = FALSE, the matricesmust all have the same number of rows, and are simply binded; therownames of the first matrix are used. See the examples.

as.matrix may be used to convert DNA sequences (of the samelength) stored in a list into a matrix while keeping the names and theclass.as.list does the reverse operation.

Value

an object of class"DNAbin" in the case ofrbind,cbind, and[.

Author(s)

Emmanuel Paradis

References

Paradis, E. (2007) A Bit-Level Coding Scheme for Nucleotides.https://emmanuelparadis.github.io/misc/BitLevelCodingScheme_20April2007.pdf

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

See Also

as.DNAbin,read.dna,read.GenBank,write.dna,image.DNAbin,AAbin

The corresponding generic functions are documented in the packagebase.

Examples

data(woodmouse)woodmouseprint(woodmouse, 15, 6)print(woodmouse[1:5, 1:300], 15, 6)### Just to show how distances could be influenced by sampling:dist.dna(woodmouse[1:2, ])dist.dna(woodmouse[1:3, ])### cbind and its options:x <- woodmouse[1:2, 1:5]y <- woodmouse[2:4, 6:10]as.character(cbind(x, y)) # gives warningas.character(cbind(x, y, fill.with.gaps = TRUE))## Not run: as.character(cbind(x, y, check.names = FALSE)) # gives an error## End(Not run)

Recode Blocks of Indels

Description

This function scans a set of aligned DNA sequences and returns amatrix with information of the localisations and lengths on alignmentgaps.

Usage

DNAbin2indel(x)

Arguments

Details

The output matrix has the same dimensions than the input one with,either a numeric value where an alignment gap starts giving the lengthof the gap, or zero. The rownames are kept.

Value

a numeric matrix.

Author(s)

Emmanuel Paradis

See Also

DNAbin,as.DNAbin,del.gaps,seg.sites,image.DNAbin,checkAlignment

Tree Estimation Based on the Minimum Evolution Algorithm

Description

The two FastME functions (balanced and OLS) perform theminimum evolution algorithm of Desper and Gascuel (2002).

Usage

  fastme.bal(X, nni = TRUE, spr = TRUE, tbr = FALSE)  fastme.ols(X, nni = TRUE)

Arguments

Details

The code to perform topology searches based on TBR (tree bisection andreconnection) did not run correctly and has been removed after therelease ofape 5.3. A warning is issued iftbr = TRUE.

Value

an object of class"phylo".

Author(s)

original C code by Richard Desper; adapted and ported to Rby Vincent Lefortvincent.lefort@lirmm.fr

References

Desper, R. and Gascuel, O. (2002) Fast and accurate phylogenyreconstruction algorithms based on the minimum-evolution principle.Journal of Computational Biology,9, 687–705.

See Also

nj,bionj,write.tree,read.tree,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)] <- xM <- t(M)M[lower.tri(M)] <- xdimnames(M) <- list(1:8, 1:8)tr <- fastme.bal(M)plot(tr, "u")### a less theoretical exampledata(woodmouse)trw <- fastme.bal(dist.dna(woodmouse))plot(trw)

Initialize a ‘corPhyl’ Structure Object

Description

Initialize acorPhyl correlation structure object.Does the same asInitialize.corStruct, but also checks the row names of data and builds an index.

Usage

## S3 method for class 'corPhyl'Initialize(object, data, ...)

Arguments

Value

An initialized object of same class asobject.

Author(s)

Julien Dutheildutheil@evolbio.mpg.de

See Also

corClasses,Initialize.corStruct.

Theoretical Lineage-Through Time Plots

Description

This function draws the lineage-through time (LTT) plots predictedunder a speciation-extinction model (aka birth-death model) withspecified values of speciation and extinction rates (which may varywith time).

A prediction interval is plotted by default which requires to define asample size (100 by default), and different curves can be combined.

Usage

LTT(birth = 0.1, death = 0, N = 100, Tmax = 50, PI = 95,    scaled = TRUE, eps = 0.1, add = FALSE, backward = TRUE,    ltt.style = list("black", 1, 1), pi.style = list("blue", 1, 2), ...)

Arguments

Details

For the moment, this works well whenbirth anddeath areconstant. Some improvements are under progress for time-dependentrates (but see below for an example).

Author(s)

Emmanuel Paradis

References

Hallinan, N. (2012) The generalized time variable reconstructedbirth–death process.Journal of Theoretical Biology,300, 265–276.

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

Paradis, E. (2015) Random phylogenies and the distribution ofbranching times.Journal of Theoretical Biology,387,39–45.

See Also

ltt.plot

Examples

### predicted LTT plot under a Yule model with lambda = 0.1### and 50 species after 50 units of time...LTT(N = 50)### ... and with a birth-death model with the same rate of### diversification (try with N = 500):LTT(0.2, 0.1, N = 50, PI = 0, add = TRUE, ltt.style = list("red", 2, 1))### predictions under different tree sizes:layout(matrix(1:4, 2, 2, byrow = TRUE))for (N in c(50, 100, 500, 1000)) {    LTT(0.2, 0.1, N = N)    title(paste("N =", N))}layout(1)## Not run: ### speciation rate decreasing with timebirth.logis <- function(t) 1/(1 + exp(0.02 * t + 4))LTT(birth.logis)LTT(birth.logis, 0.05)LTT(birth.logis, 0.1)## End(Not run)

Most Parsimonious Reconstruction

Description

This function does ancestral character reconstruction by parsimony asdescribed in Hanazawa et al. (1995) and modified by Narushima andHanazawa (1997).

Usage

MPR(x, phy, outgroup)

Arguments

Details

Hanazawa et al. (1995) and Narushima and Hanazawa (1997) used Farris's(1970) and Swofford and Maddison's (1987) framework to reconstructancestral states using parsimony. The character is assumed to takeinteger values. The algorithm finds the sets of values for each nodeas intervals with lower and upper values.

It is recommended to root the tree with the outgroup before theanalysis, so plotting the values withnodelabels issimple.

Value

a matrix of integers with two columns named “lower” and “upper”giving the lower and upper values of the reconstructed sets for eachnode.

Author(s)

Emmanuel Paradis

References

Farris, J. M. (1970) Methods for computing Wagner trees.Systematic Zoology,19, 83–92.

Hanazawa, M., Narushima, H. and Minaka, N. (1995) Generating mostparsimonious reconstructions on a tree: a generalization of theFarris–Swofford–Maddison method.Discrete AppliedMathematics,56, 245–265.

Narushima, H. and Hanazawa, M. (1997) A more efficient algorithm forMPR problems in phylogeny.Discrete Applied Mathematics,80, 231–238.

Swofford, D. L. and Maddison, W. P. (1987) Reconstructing ancestralcharacter states under Wagner parsimony.MathematicalBiosciences,87, 199–229.

See Also

ace,root,nodelabels

Examples

## the example in Narushima and Hanazawa (1997):tr <- read.tree(text = "(((i,j)c,(k,l)b)a,(h,g)e,f)d;")x <- c(1, 3, 0, 6, 5, 2, 4)names(x) <- letters[6:12](o <- MPR(x, tr, "f"))plot(tr)nodelabels(paste0("[", o[, 1], ",", o[, 2], "]"))tiplabels(x[tr$tip.label], adj = -2)## some random data:x <- rpois(30, 1)tr <- rtree(30, rooted = FALSE)MPR(x, tr, "t1")

Moran's I Autocorrelation Index

Description

This function computes Moran's I autocorrelation coefficient ofx giving a matrix of weights using the method described byGittleman and Kot (1990).

Usage

  Moran.I(x, weight, scaled = FALSE, na.rm = FALSE,          alternative = "two.sided")

Arguments

Details

The matrixweight is used as “neighbourhood” weights, andMoran's I coefficient is computed using the formula:

I = \frac{n}{S_0} \frac{\sum_{i=1}^n\sum_{j=1}^n w_{i,j}(y_i - \overline{y})(y_j - \overline{y})}{\sum_{i=1}^n {(y_i -\overline{y})}^2}

with

• y_i = observations

• w_{i,j} = distance weight

• n = number of observations

• S_0 =\sum_{i=1}^n\sum_{j=1}^n wij

The null hypothesis of no phylogenetic correlation is tested assumingnormality of I under this null hypothesis. If the observed valueof I is significantly greater than the expected value, then the valuesofx are positively autocorrelated, whereas if Iobserved <Iexpected, this will indicate negative autocorrelation.

Value

A list containing the elements:

Author(s)

Julien Dutheildutheil@evolbio.mpg.de andEmmanuel Paradis

References

Gittleman, J. L. and Kot, M. (1990) Adaptation: statistics and a nullmodel for estimating phylogenetic effects.Systematic Zoology,39, 227–241.

See Also

weight.taxo

Examples

tr <- rtree(30)x <- rnorm(30)## weights w[i,j] = 1/d[i,j]:w <- 1/cophenetic(tr)## set the diagonal w[i,i] = 0 (instead of Inf...):diag(w) <- 0Moran.I(x, w)Moran.I(x, w, alt = "l")Moran.I(x, w, alt = "g")Moran.I(x, w, scaled = TRUE) # usualy the same

Construction of Consensus Distance Matrix With SDM

Description

This function implements the SDM method of Criscuolo et al. (2006) fora set of n distance matrices.

Usage

SDM(...)

Arguments

Details

Reconstructs a consensus distance matrix from a set of input distancematrices on overlapping sets of taxa. Potentially missing values inthe supermatrix are represented byNA. An error is returned ifthe input distance matrices can not resolve to a consensus matrix.

Value

a 2-element list containing a distance matrix labelled by the union ofthe set of taxa of the input distance matrices, and a variance matrixassociated to the returned distance matrix.

Author(s)

Andrei Popescu

References

Criscuolo, A., Berry, V., Douzery, E. J. P. , and Gascuel, O. (2006)SDM: A fast distance-based approach for (super)tree building inphylogenomics.Systematic Biology,55, 740–755.

See Also

bionj,fastme,njs,mvrs,triangMtd

Ancestral Character Estimation

Description

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

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

logLik,deviance, andAIC are generic functionsused to extract the log-likelihood, the deviance, or the Akaikeinformation criterion of a fitted object. If no such values areavailable,NULL is returned.

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

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

Iftype = "continuous", the default model is Brownian motionwhere characters evolve randomly following a random walk. This modelcan be fitted by residual maximum likelihood (the default), maximumlikelihood (Felsenstein 1973, Schluter et al. 1997), least squares(method = "pic", Felsenstein 1985), or generalized leastsquares (method = "GLS", Martins and Hansen 1997, Cunningham etal. 1998). In the last case, the specification ofphy andmodel are actually ignored: it is instead given through acorrelation structure with the optioncorStruct.

In the settingmethod = "ML" andmodel = "BM" (this usedto be the default untilape 3.0-7) the maximum likelihoodestimation is done simultaneously on the ancestral values and thevariance of the Brownian motion process; these estimates are then usedto compute the confidence intervals in the standard way. The REMLmethod first estimates the ancestral value at the root (aka, thephylogenetic mean), then the variance of the Brownian motion processis estimated by optimizing the residual log-likelihood. The ancestralvalues are finally inferred from the likelihood function giving thesetwo parameters. Ifmethod = "pic" or"GLS", theconfidence intervals are computed using the expected variances underthe model, so they depend only on the tree.

It could be shown that, with a continous character, REML results inunbiased estimates of the variance of the Brownian motion processwhile ML gives a downward bias. Therefore the former is recommanded.

For discrete characters (type = "discrete"), only maximumlikelihood estimation is available (Pagel 1994) (seeMPRfor an alternative method). The model is specified through a numericmatrix with integer values taken as indices of the parameters. Thenumbers of rows and of columns of this matrix must be equal, and aretaken to give the number of states of the character. For instance,matrix(c(0, 1, 1, 0), 2) will represent a model with twocharacter 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 oftransition (the diagonal is always ignored). There are short-cuts tospecify these models:"ER" is an equal-rates model (e.g., thefirst and third examples above),"ARD" is anall-rates-different model (the second example), and"SYM" is asymmetrical 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 fromthe data.

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

With discrete characters it is necessary to compute the exponential ofthe rate matrix. The only possibility untilape 3.0-7 was thefunctionmatexpo inape. Ifuse.expm = TRUEanduse.eigen = FALSE, the functionexpm,in the package of the same name, is used.matexpo is faster butquite inaccurate for large and/or asymmetric matrices. In case ofdoubt, use the latter. Sinceape 3.0-10, it is possible to usean eigen decomposition avoiding the need to compute the matrixexponential; see details in Lebl (2013, sect. 3.8.3). This is muchfaster and is now the default.

Since version 5.2 ofape,ace can take state uncertaintyfor discrete characters into account: this should be coded withR'sNA only. More details:

https://www.mail-archive.com/r-sig-phylo@r-project.org/msg05286.html

Value

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

Note

Liam Revell points out that for discrete characters the ancestrallikelihood values returned withmarginal = FALSE are actuallythe marginal estimates, while settingmarginal = TRUE returnsthe 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 criticalreappraisal.Trends in Ecology & Evolution,13,361–366.

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

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

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

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

Martins, E. P. and Hansen, T. F. (1997) Phylogenies and thecomparative method: a general approach to incorporating phylogeneticinformation into the analysis of interspecific data.AmericanNaturalist,149, 646–667.

Pagel, M. (1994) Detecting correlated evolution on phylogenies: ageneral method for the comparative analysis of discretecharacters.Proceedings of the Royal Society of London. SeriesB. Biological Sciences,255, 37–45.

Pupko, T., Pe'er, I, Shamir, R., and Graur, D. (2000) A fast algorithmfor 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 packagephangorn (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)

Add a Scale Bar to a Phylogeny Plot

Description

This function adds a horizontal bar giving the scale of the branchlengths to a plot of a phylogenetic tree on the current graphicaldevice.

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 onthe direction of the tree. Otherwise bothx andy mustbe specified (if only one is given it is ignored).

The further arguments (...) are used to format the text. Theymay befont,cex,col, and so on (see examplesbelow, and the help page ontext).

The functionlocator may be used todetermine thex andy 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)

Incomplete Distance Matrix Filling

Description

Fills missing entries from incomplete distance matrix using theadditive 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-squaresapproach for inferring phylogenies from incomplete distancematrices.Bioinformatics,20, 2113–2121.

Alignment Explorer With Multiple Devices

Description

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

Usage

alex(x, ...)

Arguments

Details

This function works with a DNA alignment (freshly) plotted on aninteractive graphical device (i.e., not a file) withimage.After callingalex, the user clicks twice defining a rectanglein the alignment, then this portion of the alignment is extacted andplotted on anew window. The user can click as many times onthe alignment. The process is stopped by a right-click. If the userclicks twice outside the alignment, a message “Try again!” isprinted.

Each timealex is called, the alignment is plotted on a newwindow without closing or deleting those possibly already plotted.

In all cases, the device wherex is plotted is the activewindow after the operation. It shouldnot be closed during thewhole 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)

Compare DNA Sets

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 functionreturnsTRUE. Otherwise, a detailed comparison is made only ifthe labels (i.e., rownames) oftarget andcurrent are thesame (possibly in different orders). In all other cases, a briefdescription of the differences is returned (sometimes withrecommendations to make further comparisons).

This function can be used for testing in programs usingisTRUE (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 avector of mode character.

Author(s)

Emmanuel Paradis

See Also

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

Examples

data(woodmouse)woodm2 <- woodmousewoodm2[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)) # TRUEisTRUE(all.equal(woodmouse, woodm2)) # FALSEall.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)) # FALSEall.equal(DNA.clustal, DNA.muscle, TRUE)## End(Not run)

Global Comparison of two Phylogenies

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 functionall.equal for the comparison of phylogenetic trees.

A single phylogenetic tree may have several representations in the Newickformat and in the"phylo" class of objects used in ‘ape’. Oneaim of the present function is to be able to identify whether twoobjects 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 oftopologies (i.e., ignoring tip labels) of unrooted trees. This alsoaffectsunique.multiPhylo which calls the present function. See:

https://www.mail-archive.com/r-sig-phylo@r-project.org/msg01445.html.

Author(s)

Benoît Durandb.durand@alfort.AFSSA.FR

See Also

all.equal for the genericR 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)

Print DNA or AA Sequence Alignement

Description

This function displays in the console or a file an alignment of DNA orAAsequences. The first sequence is printed on the first row and thebases of the other sequences are replaced by dots if they areidentical 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 randomx <- woodmouse[, j]alview(x, showpos = FALSE) # no site position displayedalview(x, showpos = j)## Not run: alview(woodmouse, file = "woodmouse.txt")## End(Not run)

Internal Ape Functions

Description

Internalape functions which are undocumented but still exportedbecause called by other packages. Use with care!

Tools to Explore Files

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 intheir 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 witheditFileExtensions.

bydir sorts the list of files by directories.

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

Value

Xplorefiles returns a list.bydir prints the filelistings on the console.

Author(s)

Emmanuel Paradis

Examples

## Not run: x <- Xplorefiles()x # all data files on your diskbydir(x) # sorted by directoriesbydir(x["fasta"]) # only the FASTA filesXplorefiles(getwd(), recursive = FALSE) # look only in current dirXplor()## End(Not run)

Conversion Among DNA Sequence Internal Formats

Description

These functions transform a set of DNA sequences among variousinternal 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

Foras.alignment, the sequences given as argument should bestored as matrices or lists of single-character strings (the formatused inape before version 1.10). The returned object is in theformat used in the packageseqinr to store aligned sequences.

as.DNAbin is a generic function with methods so that it workswith sequences stored into vectors, matrices, or lists. It can convertsome S4 classes from the packageBiostrings in BioConductor. Forconsistency withinape, this uses an S3-style syntax. To convertobjects of class"DNAStringSetList", see the examples.

as.character is a generic function: the present methodconverts objects of class"DNAbin" into the format usedbeforeape 1.10 (matrix of single characters, or list of vectorsof single characters). This function must be used first to convertobjects 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 caseof"as.DNAbin"; a matrix of mode character or a list containingvectors 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 apealphabetFrequency(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)

Split Frequencies and Conversion Among Split Classes

Description

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

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

as.bitsplits andas.prop.part are generic functions forconverting 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, sothey work only with unrooted trees. The structure of the class"bitsplits" is described in a separate document on ape's website.

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

Value

bitsplits,as.bitsplits, andsort return an objectof 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))

Conversion Between Phylo and Matching Objects

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 givena number, and sibling groups are grouped in a “matching pair” (seeDiaconis and Holmes 1998, for details). This coding system can be usedonly for binary (fully dichotomous) trees.

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

Value

as.matching returns an object of class"matching" withthe 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 phylogenetictrees.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)mtr <- as.phylo(m)all.equal(tr, bird.orders, use.edge.length = FALSE)

Conversion Among Tree and Network Objects

Description

as.phylo is a generic function which converts an object into atree of class"phylo". There are currently two methods forobjects of class"hclust" and of class"phylog"(implemented in the packageade4). The default method is for anyobject inheriting the class"phylo" which is returned unchanged.

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

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

old2new.phylo andnew2old.phylo are utility functionsfor 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", theheight gives thedistance between the two sets that are being agglomerated. So thesedistances are divided by two when setting the branch lengths of aphylogenetic 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 labelsplot(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)

Conversion from Taxonomy Variables to Phylogenetic Trees

Description

The functionas.phylo.formula (short formas.phylo)builds a phylogenetic tree (an object of classphylo) froma 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 treewill be unresolved and will contain polytomies.

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

Value

an object of class"phylo".

Author(s)

Julien Dutheildutheil@evolbio.mpg.de, Eric Marcon andKlaus Schliep

See Also

as.phylo,read.tree for adescription of"phylo" objects,multi2di

Examples

data(carnivora)frm <- ~SuperFamily/Family/Genus/Speciestr <- 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))

Axis on Side of Phylogeny

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. Theymay befont,cex,col,las, and so on (seethe help pages onaxis andpar).

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)

Balance of a Dichotomous Phylogenetic Tree

Description

This function computes the balance of a phylogenetic tree, that is foreach node of the tree the numbers of descendants (i.e. tips) on eachof 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 thetree. The columns give the numbers of descendants on eachdaughter-branches (the order of both columns being arbitrary). If thephylogenyphy has an elementnode.label, this is used asrownames for the returned matrix; otherwise the numbers (of modecharacter) of the matrixedge ofphy are used as rownames.

Author(s)

Emmanuel Paradis

References

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

Base frequencies from DNA Sequences

Description

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

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

Ftab computes the contingency table with the absolutefrequencies 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 thesample.

ForFtab, if the argumenty is given then bothxandy are coerced as vectors and must be of equal length. Ify is not given,x must be a matrix or a list and onlythe two first sequences are used.

Value

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

Author(s)

Emmanuel Paradis

See Also

seg.sites,nuc.div (inpegas),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 aboveFtab(woodmouse[14:15, ]) # between the last two

Extended Version of the Birth-Death Models to Estimate Speciationand Extinction Rates

Description

This function fits by maximum likelihood a birth-death model to thecombined phylogenetic and taxonomic data of a given clade. Thephylogenetic data are given by a tree, and the taxonomic data by thenumber 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 overd/b andb - d, whereb is the birth rate, andd the deathrate.

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

If the argumentS has names, then they are matched to the tiplabels ofphy. The user must be careful here since the functionrequires that both series of names perfectly match, so this operationmay fail if there is a typing or syntax error. If both series of namesdo not match, the valuesS are taken to be in the same orderthan the tip labels ofphy, and a warning message is issued.

Note that the function does not check that the tree is effectivelyultrametric, so if it is not, the returned result may not bemeaningful.

Ifconditional = TRUE, the probabilities of the taxonomic dataare calculated conditioned on no extinction (Rabosky et al. 2007). Inprevious versions of the present function (until ape 2.6-1),unconditional probabilities were used resulting in underestimatedextinction rate. Though it does not make much sense to useconditional = FALSE, this option is provided to compare resultsfrom 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 phylogeneticand taxonomic data.Proceedings of the Royal Society ofLondon. 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 diversificationrates during the radiation of Australia's most diverse vertebrateclade.Proceedings of the Royal Society of London. SeriesB. 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

Time-Dependent Birth-Death Models

Description

This function fits a used-defined time-dependent birth-deathmodel.

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 theirprimitives can be found in the help page ofyule.time.

The model is fitted by minimizing the least squares deviation betweenthe observed and the predicted distributions of branching times. Thesecomputations rely heavily on numerical integrations. Iffast = FALSE, integrations are done with R'sintegratefunction. Iffast = TRUE, a faster but less accurate functionprovided inape is used. If fitting a complex model to a largephylogeny, a strategy might be to first use the latter option, andthen to use the estimates as starting values withfast = FALSE.

Value

A list with the following components:

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

• SS: the minimized sum of squares.

• convergence: output convergence criterion fromnlminb.

• message: id.

• iterations: id.

• evaluations: id.

Author(s)

Emmanuel Paradis

References

Paradis, E. (2011) Time-dependent speciation and extinction fromphylogenies: 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 modelbd.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)

Phylogenetic Generalized Linear Mixed Model for Binary Data

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)

whereV is a variance-covariance matrix derived from a phylogeny (typically under the assumption of Brownian motion evolution). Although mathematically there is no requirement forV to be ultrametric, forcingV 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 matrixV 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 modelsfor phylogenetic analyses of community structure.EcologicalMonographs,81, 511–525.

Ives, A. R. and Garland, T., Jr. (2014) Phylogenetic regression forbinary dependent variables. Pages 231–261in L. Z. Garamszegi,editor.Modern Phylogenetic Comparative Methods and TheirApplication in Evolutionary Biology. Springer-Verlag, BerlinHeidelberg.

See Also

packagepez and its functioncommunityPGLMM;packagephylolm and its functionphyloglm;packageMCMCglmm

Examples

## Illustration of binaryPGLMM() with simulated data# Generate random phylogenyn <- 100phy <- compute.brlen(rtree(n=n), method = "Grafen", power = 1)# Generate random data and standardize to have mean 0 and variance 1X1 <- rTraitCont(phy, model = "BM", sigma = 1)X1 <- (X1 - mean(X1))/var(X1)# Simulate binary Ysim.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 modelbinaryPGLMM(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 phylogenysummary(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 MCMCglmmlibrary(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.labelrownames(invV) <- sim.dat$speciesnitt <- 43000thin <- 10burnin <- 3000prior <- 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 = 1000s2 <- 0.4B1 <- 1meanEsts <- 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 = 1000s2 <- 0B0 <- 0B1 <- 0H0.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)

Binds Trees

Description

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

Usage

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

Arguments

Details

The argumentx can be seen as the receptor tree, whereasy is the donor tree. The root ofy is then grafted on alocation ofx specified bywhere and, possibly,position. Ify has a root edge, this is added as ininternal 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 ofthe other one are ignored with a warning.

If one (or both) of the trees has no branch length, it is possible tospecify a value of 'position' to graft 'y' below the node of 'x'specified by 'where'. In this case, the exact value of 'position' isnot important as long as it is greater than zero. The new node will bemultichotomous if 'y' has no root edge. This can be solved by givingan 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 orderstreefile1 <- 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 treesx <- 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 <- .2z <- 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.2z <- 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 <- .2z <- x + yplot(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)

Phylogeny of the Families of Birds From Sibley and Ahlquist

Description

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

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

Usage

data(bird.families)

Format

The data are stored as an object of class"phylo" whichstructure is described in the help page of the functionread.tree.

Source

Sibley, C. G. and Ahlquist, J. E. (1990) Phylogeny and classificationof 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)

Phylogeny of the Orders of Birds From Sibley and Ahlquist

Description

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

The present tree is based on the relationships among orders. Thebranch 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" whichstructure is described in the help page of the functionread.tree.

Source

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

See Also

read.tree,bird.families

Examples

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

Estimation of Speciation and Extinction Rates With Birth-Death Models

Description

This function fits by maximum likelihood a birth-death model to thebranching times computed from a phylogenetic tree using the method ofNee 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 modelstudied by Kendall (1948) so that the likelihood has to be maximizedoverd/b andb - d, whereb is the birth rate,andd the death rate. This is the approach used by the presentfunction.

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

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

Value

An object of class"birthdeath" which is a list with thefollowing 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 reconstructedevolutionary process.Philosophical Transactions of the RoyalSociety of London. Series B. Biological Sciences,344, 305–311.

Nee, S., Holmes, E. C., May, R. M. and Harvey, P. H. (1995) Estimatingextinctions from molecular phylogenies. inExtinction 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

Tree Bipartition and Bootstrapping Phylogenies

Description

These functions analyse bipartitions found in a series of trees.

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

prop.clades counts the number of times the bipartitions presentinphy are present in a series of trees given as... orin the list previously computed and given withpart.

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 argumentFUN inboot.phylo must be the function usedto estimate the tree from the original data matrix. Thus, if the treewas estimated with neighbor-joining (seenj), one maybe wantssomething likeFUN = function(xx) nj(dist.dna(xx)).

block inboot.phylo specifies the number of columns tobe resampled altogether. For instance, if one wants to resample at thecodon-level, thenblock = 3 must be used.

Usingcheck.labels = FALSE inprop.part decreasescomputing times. This requires that (i) all trees have the same tiplabels,and (ii) these labels are ordered similarly in alltrees (in other words, the elementtip.label are identical inall trees).

The plot function represents a contingency table of the differentpartitions (on thex-axis) in the lower panel, and their observednumbers in the upper panel. Any further arguments (...) are used tochange the aspects of the points in the lower panel: these may bepch,col,bg,cex, etc. This functionworks only if there is an attributelabels in the object.

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

Value

prop.part returns an object of class"prop.part" whichis a list with an attribute"number". The elements of this listare the observed clades, and the attribute their respectivenumbers. If the defaultcheck.labels = FALSE is used, anattribute"labels" is added, and the vectors of the returnedobject contains the indices of these labels instead of the labelsthemselves.

prop.clades andboot.phylo return a numeric vectorwhichith element is the number associated to theithnode ofphy. Iftrees = TRUE,boot.phylo returnsa list whose first element (named"BP") is like before, and thesecond element ("trees") is a list with the bootstrapedtrees.

summary returns a numeric vector.

Note

prop.clades calls internallyprop.part with the optioncheck.labels = TRUE, which may be very slow. If the treespassed as... fulfills conditions (i) and (ii) above, then itmight be faster to first call, e.g.,pp <- prop.part(...), thenuse the optionpart:prop.clades(phy, part = pp).

Sinceape 3.5,prop.clades should return sensible resultsfor all values ofrooted: ifFALSE, the numbers ofbipartitions (or splits); ifTRUE, the number of clades (ofhopefully rooted trees).

Author(s)

Emmanuel Paradis

References

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

Felsenstein, J. (1985) Confidence limits on phylogenies: an approachusing 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 phangorntree <- 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 <- 100nrep2 <- 100p <- ncol(woodmouse)DB <- 0for (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)

Branching Times of a Phylogenetic Tree

Description

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

Usage

branching.times(phy)

Arguments

Value

a numeric vector with the branching times. If the phylogenyphyhas an elementnode.label, this is used as names for thereturned vector; otherwise the numbers (of mode character) of thematrixedge ofphy are used as names.

Author(s)

Emmanuel Paradis

See Also

is.ultrametric

Building Lists of Trees

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

Thesec methods check all the arguments, and return by defaulta list of single trees unless some objects are not trees or lists oftrees, in which caserecursive is switched to FALSE and awarning message is given. Ifrecursive = FALSE, the objects aresimply concatenated into a list. Beforeape 4.0,recursivewas always set to FALSE.

.compressTipLabel transforms an object of class"multiPhylo" by checking that all trees have the same tiplabels and renumbering the tips in theedge matrix so that thetip numbers are also the same taking the first tree as the reference(duplicated labels are not allowed). The returned object has a uniquevector 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))xy <- c(rtree(4), rtree(4))z <- c(x, y)zprint(z, TRUE)try(.compressTipLabel(x)) # errora <- .compressTipLabel(y).uncompressTipLabel(a) # back to y## eventually compare str(a) and str(y)

Carnivora body sizes and life history traits

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)

Check DNA Alignments

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 theset a aligned DNA sequences. If alignment gaps are present, theirwidth distribution is analysed, as well as the width of contiguousbase segments. The pattern of nucleotide diversity on each site isalso analysed, and a relevant table is printed.

Ifplot = TRUE, four plots are done: an image of thealignement, the distribution of gap widths (if present), the Shannonindex of nucleotide diversity along the sequence, and the number ofobserved bases along the sequence.

If the sequences contain many gaps, it might be better to setcheck.gaps = FALSE to skip the analysis of contiguoussegments.

Value

NULL

Author(s)

Emmanuel Paradis

See Also

alview,image.DNAbin,all.equal.DNAbin

Examples

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

Checking Labels

Description

Checking and correcting character strings, particularly before writinga Newick tree.

Usage

checkLabel(x)

Arguments

Details

This function deletes the leading and trailing spaces (includingtabulations, new lines, and left or right parentheses at the beginningor end of the strings), substitutes the spaces inside the strings byunderscores, and substitutes commas, colons, semicolons, andparentheses 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)   ")

Check the Structure of a "phylo" Object

Description

This function takes as single argument an object (phy), checks itselements, and prints a diagnostic. All problems are printed with alabel: FATAL (will likely cause an error or a crash) or MODERATE (maycause 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] <- 0checkValidPhylo(tr)

Number of Cherries and Null Models of Trees

Description

This function calculates the number of cherries (see definition below)on a phylogenetic tree, and tests the null hypotheses whether thisnumber agrees with those predicted from two null models of trees (theYule 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 eitherrooted or unrooted, but the present function considers only rootedtrees. The probability distribution function of the number of cherrieson a tree depends on the speciation/extinction model that generatedthe tree.

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

The probabilities are computed using recursive formulae; however, forboth models, the probability density function converges to a normallaw with increasing number of tips in the tree. The function usesthese normal approximations for a number of tips greater than or equalto 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 twomodels of trees.Mathematical Biosciences,164, 81–92.

See Also

gammaStat

Bat Phylogeny

Description

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

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. andSimmons, N. B. (2002) A phylogenetic supertree of the bats (Mammalia:Chiroptera).Biological Reviews of the Cambridge PhilosophicalSociety,77, 223–259.

See Also

read.nexus,zoom

Examples

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

Molecular Dating With Mean Path Lengths

Description

This function estimates the node ages of a tree using the mean pathlengths method of Britton et al. (2002). The branch lengths of theinput 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 withthe mean of the distances from this node to all tips descending fromit. Under the assumption of a molecular clock, standard-errors of theestimates node ages can be computed (Britton et al. 2002).

The tests performed iftest = TRUE is a comparison of the MPLof the two subtrees originating from a node; the null hypothesis isthat the rate of substitution was the same in both subtrees (Brittonet al. 2002). The test statistic follows, under the null hypothesis, astandard normal distribution. The returnedP-value is theprobability of observing a greater absolute value (i.e., a two-sidedtest). No correction for multiple testing is applied: this is left tothe user.

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

Value

an object of class"phylo" with branch lengths as estimated bythe 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 pathlengths.Molecular Phylogenetics and Evolution,24,58–65.

See Also

chronopl

Examples

tr <- rtree(10)tr$edge.length <- 5*tr$edge.lengthchr <- 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)

Molecular Dating With Penalized Likelihood

Description

This function estimates the node ages of a tree using asemi-parametric method based on penalized likelihood (Sanderson2002). The branch lengths of the input tree are interpreted as meannumbers 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 parametricformulation where each branch has its own rate, and a nonparametricterm where changes in rates are minimized between contiguousbranches. A smoothing parameter (lambda) controls this trade-off. Iflambda = 0, then the parametric component dominates and rates vary asmuch as possible among branches, whereas for increasing values oflambda, the variation are smoother to tend to a clock-like model (samerate for all branches).

lambda must be given. The known ages are given inage.min, and the correponding node numbers innode.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 nodesare estimated.

Ifage.max = NULL (the default), it is assumed thatage.min gives exactly known ages. Otherwise,age.max andage.min must be of the same length and give the intervals foreach node. Some node may be known exactly while the others areknown within some bounds: the values will be identical in botharguments 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 10units 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) andhave the same values ofage.min andage.max (say, 10 and15) this will result in an error because the medians of these valuesare 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 thegiven values, for instance useage.max = 14.9 for the youngestnode, orage.max = 15.1 for the oldest one (or similarly forage.min).

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

The cross-validation used here is different from the one proposed bySanderson (2002). Here, each tip is dropped successively and theanalysis is repeated with the reduced tree: the estimated dates forthe remaining nodes are compared with the estimates from the fulldata. For theith tip the following is calculated:

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

,

wheret_j is the estimated date for thejth nodewith the full phylogeny,t_j^{-i} is the estimated datefor thejth node after removing tipi from the tree,andn is the number of tips.

The present version uses thenlminb to optimisethe penalized likelihood function: see its help page for details onparameters controlling the optimisation procedure.

Value

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

Note

The new functionchronos replaces the present one whichis no more maintained.

Author(s)

Emmanuel Paradis

References

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

See Also

chronos,chronoMPL

Molecular Dating by Penalised Likelihood and Maximum Likelihood

Description

chronos is the main function fitting a chronogram to aphylogenetic tree whose branch lengths are in number of substitutionper sites.

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

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

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 replaceschronopl but with a different interfaceand some extensions (see References).

The known dates (argumentcalibration) must be given in a dataframe with the following column names: node, age.min, age.max, andsoft.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 valuespecifying whether the bounds are soft. If age.min = age.max, thismeans that the age is exactly known. This data frame can be built withmakeChronosCalib which returns by default a data frame with asingle row giving age = 1 for the root. The data frame can be builtinteractively by clicking on the plotted tree.

The argumentcontrol allows one to change some parameters ofthe optimisation procedure. This must be a list with names. Theavailable options with their default values are:

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

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

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

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

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

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

Usingmodel = "clock" is actually a short-cut tomodel = "discrete" and settingnb.rate.cat = 1 in the list passed tocontrol.

The commandchronos.control() returns a list with the defaultvalues of these parameters. They may be modified by passing them tothis function, or directly in the list.

Value

chronos returns an object of classc("chronos", "phylo"). There is a print method for it. There are additionalattributes which can be visualised withstr or extracted withattr.

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 phylogeneticinference: bridging the parsimony-likelihood gap.SystematicBiology,57, 665–674.

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

Sanderson, M. J. (2002) Estimating absolute rates of molecularevolution and divergence times: a penalized likelihoodapproach.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)

Multiple Sequence Alignment with External Applications

Description

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

This version includes an experimental version ofmuscle5 whichcalls MUSCLE5 (see the link to the documentation in the Referencesbelow);muscle still calls MUSCLE version 3. Note that theexecutable of MUSCLE5 is also named ‘muscle’ by the defaultcompilation setting.

The functionsefastats andletterconf 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 thatthey are in a directory located on the PATH (i.e., accessible from anyother directory). Alternatively, the full path to the executablemay be given (e.g.,exec = "~/muscle/muscle"), or a (symbolic)link may be copied in the working directory. For Debian and itsderivatives (e.g., Ubuntu), it is recommended to use the binariesdistributed by Debian.

clustal tries to guess the name of the executable programdepending on the operating system. Specifically, the followings areused: “clustalw” under Linux, “clustalw2” under MacOS, and“clustalw2.exe” under Windows. Forclustalomega,“clustalo[.exe]” is the default on all systems (with no specificpath).

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

Sinceape 5.1,clustal,clustalomega, andmuscle can align AA sequences as well as DNA sequences.

Value

an object of class"DNAbin" or"AAbin" with the alignedsequences.

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 alignmentwith the Clustal series of programs.Nucleic Acids Research31, 3497–3500.http://www.clustal.org/

Edgar, R. C. (2004) MUSCLE: Multiple sequence alignment with highaccuracy 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 novelmethod for multiple sequence alignments.Journal of MolecularBiology,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 ofhigh-quality protein multiple sequence alignments using ClustalOmega.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)

Coalescent Intervals

Description

This function extracts or generates information about coalescent intervals(number of lineages, interval lengths, interval count, total depth) froma phylogenetic tree or a list of internode distances. The input treeneeds 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 formattree.hiv <- read.tree(text = hivtree.newick) # load treeci <- coalescent.intervals(tree.hiv) # from treecidata("hivtree.table") # same tree, but in table formatci <- coalescent.intervals(hivtree.table$size) # from vector of interval lengthsci

Collapse Single Nodes

Description

collapse.singles deletes the single nodes (i.e., with a singledescendant) 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")trhas.singles(tr)## the following shows that node #4 (ie, the root) is a singleton## and node #6 is the first bifurcating nodetr$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'

Collapsed Coalescent Intervals

Description

This function takes a"coalescentIntervals" objects and collapses neighbouringcoalescent intervals into a single combined interval so that every collapsed interval islarger thanepsilon. Collapsed coalescent intervals are used, e.g., to obtain thegeneralized skyline plot (skyline). Forepsilon = 0 no intervalis collapsed.

Usage

collapsed.intervals(ci, epsilon=0)

Arguments

Details

Proceeding from the tips to the root of the tree each smallinterval is pooled with the neighboring interval closer to the root. If theneighboring interval is also small, then pooling continues until the compositeinterval is larger thanepsilon. Note that this approach prevents theoccurrence of zero-length intervals at the present.For more details see Strimmer and Pybus (2001).

Value

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

Author(s)

Korbinian Strimmer

References

Strimmer, K. and Pybus, O. G. (2001) Exploring the demographic historyof DNA sequences using the generalized skyline plot.MolecularBiology and Evolution,18, 2298–2305.

See Also

coalescent.intervals,skyline.

Examples

data("hivtree.table") # example tree# colescent intervals from vector of interval lengthsci <- coalescent.intervals(hivtree.table$size)ci# collapsed intervalscl1 <- collapsed.intervals(ci,0)cl2 <- collapsed.intervals(ci,0.0119)cl1cl2

Cheverud's Comparative Method

Description

This function computes the phylogenetic variance component and theresidual deviation for continous characters, taking into account thephylogenetic relationships among species, following the comparativemethod described in Cheverud et al. (1985). The correction proposed byRholf (2001) is used.

Usage

compar.cheverud(y, W, tolerance = 1e-06, gold.tol = 1e-04)

Arguments

Details

Model:

y = \rho W y + e

wheree is the error term, assumed to be normally distributed.\rho is estimated by the maximum likelihood procedure givenin Rohlf (2001), using a golden section search algorithm. The code ofthis function is indeed adapted from a MatLab code given in appendixin Rohlf's article, to correct a mistake in Cheverud's original paper.

Value

A list with the following components:

Author(s)

Julien Dutheildutheil@evolbio.mpg.de

References

Cheverud, J. M., Dow, M. M. and Leutenegger, W. (1985) The quantitativeassessment of phylogenetic constraints in comparative analyses: sexualdimorphism in body weight among primates.Evolution,39, 1335–1351.

Rohlf, F. J. (2001) Comparative methods for the analysis of continuousvariables: geometric interpretations.Evolution,55,2143–2160.

Harvey, P. H. and Pagel, M. D. (1991)The Comparative Method inEvolutionary Biology. Oxford University Press.

See Also

compar.lynch

Examples

### Example from Harvey and Pagel's book:y<-c(10,8,3,4)W <- matrix(c(1,1/6,1/6,1/6,1/6,1,1/2,1/2,1/6,1/2,1,1,1/6,1/2,1,1), 4)compar.cheverud(y,W)### Example from Rohlf's 2001 article:W<- matrix(c(  0,1,1,2,0,0,0,0,  1,0,1,2,0,0,0,0,  1,1,0,2,0,0,0,0,  2,2,2,0,0,0,0,0,  0,0,0,0,0,1,1,2,  0,0,0,0,1,0,1,2,  0,0,0,0,1,1,0,2,  0,0,0,0,2,2,2,0),8)W <- 1/WW[W == Inf] <- 0y<-c(-0.12,0.36,-0.1,0.04,-0.15,0.29,-0.11,-0.06)compar.cheverud(y,W)

Comparative Analysis with GEEs

Description

compar.gee performs the comparative analysis using generalizedestimating equations as described by Paradis and Claude (2002).

drop1 tests single effects of a fitted model output fromcompar.gee.

predict returns the predicted (fitted) values of the model.

Usage

compar.gee(formula, data = NULL, family = "gaussian", phy, corStruct,          scale.fix = FALSE, scale.value = 1)## S3 method for class 'compar.gee'drop1(object, scope, quiet = FALSE, ...)## S3 method for class 'compar.gee'predict(object, newdata = NULL, type = c("link", "response"), ...)

Arguments

Details

If a data frame is specified for the argumentdata, then itsrownames are matched to the tip labels ofphy. The user must becareful here since the function requires that both series of namesperfectly match, so this operation may fail if there is a typing orsyntax error. If both series of names do not match, the values in thedata frame are taken to be in the same order than the tip labels ofphy, and a warning message is issued.

Ifdata = NULL, then it is assumed that the variables are inthe same order than the tip labels ofphy.

Value

compar.gee returns an object of class"compar.gee" withthe following components:

drop1 returns an object of class"anova".

predict returns a vector or a data frame ifnewdata is used.

Note

The calculation of the phylogenetic degrees of freedom is likely to beapproximative for non-Brownian correlation structures (this will berefined soon).

The calculation of the quasilikelihood information criterion (QIC)needs to be tested.

Author(s)

Emmanuel Paradis

References

Pan, W. (2001) Akaike's information criterion in generalizedestimating equations.Biometrics,57, 120–125.

Paradis, E. and Claude J. (2002) Analysis of comparative data usinggeneralized estimating equations.Journal of theoreticalBiology,218, 175–185.

See Also

read.tree,pic,compar.lynch,drop1

Examples

### The example in Phylip 3.5c (originally from Lynch 1991)### (the same analysis than in help(pic)...)tr <- "((((Homo:0.21,Pongo:0.21):0.28,Macaca:0.49):0.13,Ateles:0.62):0.38,Galago:1.00);"tree.primates <- read.tree(text = tr)X <- c(4.09434, 3.61092, 2.37024, 2.02815, -1.46968)Y <- c(4.74493, 3.33220, 3.36730, 2.89037, 2.30259)### Both regressions... the results are quite close to those obtained### with pic().compar.gee(X ~ Y, phy = tree.primates)compar.gee(Y ~ X, phy = tree.primates)### Now do the GEE regressions through the origin: the results are quite### different!compar.gee(X ~ Y - 1, phy = tree.primates)compar.gee(Y ~ X - 1, phy = tree.primates)

Lynch's Comparative Method

Description

This function computes the heritable additive value and the residualdeviation for continous characters, taking into account thephylogenetic relationships among species, following the comparativemethod described in Lynch (1991).

Usage

compar.lynch(x, G, eps = 1e-4)

Arguments

Details

The parameter estimates are computed following the EM(expectation-maximization) algorithm. This algorithm usually leads toconvergence but may lead to local optima of the likelihoodfunction. It is recommended to run several times the function in orderto detect these potential local optima. The ‘optimal’ value foreps depends actually on the range of the data and may bechanged by the user in order to check the stability of the parameterestimates. Convergence occurs when the differences between twosuccessive iterations of the EM algorithm leads to differences betweenboth residual and additive values less than or equal toeps.

Value

A list with the following components:

Note

The present function does not perform the estimation of ancestralphentoypes as proposed by Lynch (1991). This will be implemented ina future version.

Author(s)

Julien Claudejulien.claude@umontpellier.fr

References

Lynch, M. (1991) Methods for the analysis of comparative data inevolutionary biology.Evolution,45, 1065–1080.

See Also

pic,compar.gee

Examples

### The example in Lynch (1991)x <- "((((Homo:0.21,Pongo:0.21):0.28,Macaca:0.49):0.13,Ateles:0.62):0.38,Galago:1.00);"tree.primates <- read.tree(text = x)X <- c(4.09434, 3.61092, 2.37024, 2.02815, -1.46968)Y <- c(4.74493, 3.33220, 3.36730, 2.89037, 2.30259)compar.lynch(cbind(X, Y),             G = vcv.phylo(tree.primates, cor = TRUE))

Ornstein–Uhlenbeck Model for Continuous Characters

Description

This function fits an Ornstein–Uhlenbeck model giving a phylogenetictree, and a continuous character. The user specifies the node(s) wherethe optimum changes. The parameters are estimated by maximumlikelihood; their standard-errors are computed assuming normality ofthese estimates.

Usage

compar.ou(x, phy, node = NULL, alpha = NULL)

Arguments

Details

The Ornstein–Uhlenbeck (OU) process can be seen as a generalizationof the Brownian motion process. In the latter, characters are assumedto evolve randomly under a random walk, that is change is equallylikely in any direction. In the OU model, change is more likelytowards the direction of an optimum (denoted\theta) witha strength controlled by a parameter denoted\alpha.

The present function fits a model where the optimum parameter\theta, is allowed to vary throughout the tree. This isspecified with the argumentnode:\theta changesafter each node whose number is given there. Note that the optimumchangesonly for the lineages which are descendants of thisnode.

Hansen (1997) recommends to not estimate\alpha togetherwith the other parameters. The present function allows this by givinga numeric value to the argumentalpha. By default, thisparameter is estimated, but this seems to yield very largestandard-errors, thus validating Hansen's recommendation. In practice,a “poor man estimation” of\alpha can be done byrepeating the function call with different values ofalpha, andselecting the one that minimizes the deviance (see Hansen 1997 for anexample).

Ifx has names, its values are matched to the tip labels ofphy, otherwise its values are taken to be in the same orderthan the tip labels ofphy.

The user must be careful here since the function requires that bothseries of names perfectly match, so this operation may fail if thereis a typing or syntax error. If both series of names do not match, thevalues in thex are taken to be in the same order than the tiplabels ofphy, and a warning message is issued.

Value

an object of class"compar.ou" which is list with the followingcomponents:

Note

The inversion of the variance-covariance matrix in the likelihoodfunction appeared as somehow problematic. The present implementationuses a Cholevski decomposition with the functionchol2inv instead of the usual functionsolve.

Author(s)

Emmanuel Paradis

References

Hansen, T. F. (1997) Stabilizing selection and the comparativeanalysis of adaptation.Evolution,51, 1341–1351.

See Also

ace,compar.lynch,corBrownian,corMartins,pic

Examples

data(bird.orders)### This is likely to give you estimates close to 0, 1, and 0### for alpha, sigma^2, and theta, respectively:compar.ou(x <- rnorm(23), bird.orders)### Much better with a fixed alpha:compar.ou(x, bird.orders, alpha = 0.1)### Let us 'mimick' the effect of different optima### for the two clades of birds...x <- c(rnorm(5, 0), rnorm(18, 5))### ... the model with two optima:compar.ou(x, bird.orders, node = 25, alpha = .1)### ... and the model with a single optimum:compar.ou(x, bird.orders, node = NULL, alpha = .1)### => Compare both models with the difference in deviances##     which follows a chi^2 with df = 1.

Compare Two "phylo" Objects

Description

This function compares two phylogenetic trees, rooted or unrooted, andreturns a detailed report of this comparison.

Usage

comparePhylo(x, y, plot = FALSE, force.rooted = FALSE,             use.edge.length = FALSE, commons = TRUE,             location = "bottomleft", ...)## S3 method for class 'comparePhylo'print(x, ...)

Arguments

Details

In all cases, the numbers of tips and of nodes and the tip labels arecompared.

If both trees are rooted, or ifforce.rooted = TRUE, the cladecompositions of each tree are compared. If both trees are alsoultrametric, their branching times are compared.

If both trees are unrooted and have the same number of nodes, thebipartitions (aka splits) are compared.

Ifplot = TRUE, the edge lengths are not used by defaultbecause in some situations with unrooted trees, some splits might notbe visible if the corresponding internal edge length is very short. Touse edge lengths, setuse.edge.length = TRUE.

Value

an object of class"comparePhylo" which is a list with messagesfrom the comparison and, optionally, tables comparing branching times.

Author(s)

Emmanuel Paradis, Klaus Schliep

See Also

all.equal.phylo

Examples

## two unrooted trees but force comparison as rooted:a <- read.tree(text = "(a,b,(c,d));")b <- read.tree(text = "(a,c,(b,d));")comparePhylo(a, b, plot = TRUE, force.rooted = TRUE)## two random unrooted trees:c <- rtree(5, rooted = FALSE)d <- rtree(5, rooted = FALSE)comparePhylo(c, d, plot = TRUE)

Branch Lengths Computation

Description

This function computes branch lengths of a tree using differentmethods.

Usage

compute.brlen(phy, method = "Grafen", power = 1, ...)

Arguments

Details

Grafen's (1989) computation of branch lengths: each node is given a‘height’, namely the number of leaves of the subtree minus one, 0 forleaves. Each height is scaled so that root height is 1, and thenraised at power 'rho' (> 0). Branch lengths are then computed as thedifference between height of lower node and height of upper node.

If one or several numeric values are provided asmethod, theyare recycled if necessary. If a function is given instead, furtherarguments are given in place of... (they must be named, seeexamples).

Zero-length branches are not treated as multichotomies, and thus mayneed to be collapsed (seedi2multi).

Value

An object of classphylo with branch lengths.

Author(s)

Julien Dutheildutheil@evolbio.mpg.de andEmmanuel Paradis

References

Grafen, A. (1989) The phylogenetic regression.PhilosophicalTransactions of the Royal society of London. Series B. BiologicalSciences,326, 119–157.

See Also

read.tree for a description ofphylo objects,di2multi,multi2di

Examples

data(bird.orders)plot(compute.brlen(bird.orders, 1))plot(compute.brlen(bird.orders, runif, min = 0, max = 5))layout(matrix(1:4, 2, 2))plot(compute.brlen(bird.orders, power=1), main=expression(rho==1))plot(compute.brlen(bird.orders, power=3), main=expression(rho==3))plot(compute.brlen(bird.orders, power=0.5), main=expression(rho==0.5))plot(compute.brlen(bird.orders, power=0.1), main=expression(rho==0.1))layout(1)

Compute and Set Branching Times

Description

This function computes the branch lengths of a tree giving itsbranching times (aka node ages or heights).

Usage

compute.brtime(phy, method = "coalescent", force.positive = NULL)

Arguments

Details

By default, a set of random branching times is generated from a simplecoalescent, and the optionforce.positive is set toTRUEso that no branch length is negative.

If a numeric vector is passed tomethod, it is taken as thebranching times of the nodes with respect to their numbers (i.e., thefirst element ofmethod is the branching time of the nodenumberedn + 1 [= the root], the second element of the nodenumberedn + 2, and so on), soforce.positive is set toFALSE. This may result in negative branch lengths. To avoidthis, one should useforce.positive = TRUE in which case thebranching times are eventually reordered.

Value

An object of class"phylo" with branch lengths and ultrametric.

Author(s)

Emmanuel Paradis

See Also

compute.brlen,branching.times

Examples

tr <- rtree(10)layout(matrix(1:4, 2))plot(compute.brtime(tr)); axisPhylo()plot(compute.brtime(tr, force.positive = FALSE)); axisPhylo()plot(compute.brtime(tr, 1:9)); axisPhylo() # a bit nonsenseplot(compute.brtime(tr, 1:9, TRUE)); axisPhylo()layout(1)

Concensus Trees

Description

Given a series of trees, this function returns the consensus tree. Bydefault, the strict-consensus tree is computed. To get themajority-rule consensus tree, usep = 0.5. Any value between0.5 and 1 can be used.

Usage

consensus(..., p = 1, check.labels = TRUE, rooted = FALSE)

Arguments

Details

Usingcheck.labels = FALSE results inconsiderable decrease in computing times. This requires that alltrees have the same tip labels,and these labels areordered similarly in all trees (in other words, the elementtip.label are identical in all trees).

Untilape 5.6-2, the trees passed to this function wereimplicitly treated as rooted, even when the optionrooted = FALSE was used. This is now fixed (see PR65 on GitHub) so that, bydefault, the trees are explicitly treated as unrooted (even ifis.rooted returnsTRUE). Thus, it couldbe that results now differ from previous analyses (settingrooted = TRUE might help to replicate previous results).

Value

an object of class"phylo".

Author(s)

Emmanuel Paradis

References

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

See Also

prop.part,dist.topo

Pairwise Distances from a Phylogenetic Tree

Description

cophenetic.phylo computes the pairwise distances between thepairs of tips from a phylogenetic tree using its branch lengths.

dist.nodes does the same but between all nodes, internal andterminal, of the tree.

Usage

## S3 method for class 'phylo'cophenetic(x)dist.nodes(x, fail.if.no.length = FALSE)

Arguments

Value

a numeric matrix with colnames and rownames set to the names of thetips (as given by the elementtip.label of the argumentphy), or, in the case ofdist.nodes, the numbers of thetips and the nodes (as given by the elementedge).

Author(s)

Emmanuel Paradis

See Also

read.tree to read tree files in Newick format,cophenetic for the generic function

Plots two phylogenetic trees face to face with links between the tips.

Description

This function plots two trees face to face with the links if specified. It is possible to rotate the branches of each tree around the nodes by clicking.

Usage

cophyloplot(x, y, assoc = NULL, use.edge.length = FALSE, space = 0,       length.line = 1, gap = 2, type = "phylogram", rotate = FALSE,       col = par("fg"), lwd = par("lwd"), lty = par("lty"),       show.tip.label = TRUE, font = 3, ...)

Arguments

Details

The aim of this function is to plot simultaneously two phylogenetic trees with associated taxa. The two trees do not necessarily have the same number of tips and more than one tip in one phylogeny can be associated with a tip in the other.

The association matrix used to draw the links has to be a matrix with two columns containing the names of the tips. One line in the matrix represents one link on the plot. The first column of the matrix has to contain tip labels of the first tree (phy1) and the second column of the matrix, tip labels of the second tree (phy2). There is no limit (low or high) for the number of lines in the matrix. A matrix with two colums and one line will give a plot with one link.

Argumentsgap,length.line andspace have to be changed to get a nice plot of the two phylogenies. Note that the function takes into account the length of the character strings corresponding to the names at the tips, so that the lines do not overwrite those names.

Therotate argument can be used to transform both phylogenies in order to get the more readable plot (typically by decreasing the number of crossing lines). This can be done by clicking on the nodes. The escape button or right click take back to the console.

Author(s)

Damien de Viennedamien.de-vienne@u-psud.fr

See Also

plot.phylo,rotate,rotateConstr

Examples

#two random treestree1 <- rtree(40)tree2 <- rtree(20)#creation of the association matrix:association <- cbind(tree2$tip.label, tree2$tip.label)cophyloplot(tree1, tree2, assoc = association,            length.line = 4, space = 28, gap = 3)#plot with rotations## Not run: cophyloplot(tree1, tree2, assoc=association, length.line=4, space=28, gap=3, rotate=TRUE)## End(Not run)

Blomberg et al.'s Correlation Structure

Description

The “ACDC” (accelerated/decelerated) model assumes that continuoustraits evolve under a Brownian motion model which rates accelerates(ifg < 1) or decelerates (ifg > 1) throughtime. Ifg = 1, then the model reduces to a Brownian motionmodel.

Usage

corBlomberg(value, phy, form = ~1, fixed = FALSE)## S3 method for class 'corBlomberg'corMatrix(object, covariate = getCovariate(object),                   corr = TRUE, ...)## S3 method for class 'corBlomberg'coef(object, unconstrained = TRUE, ...)

Arguments

Value

an object of class"corBlomberg", the coefficients from anobject of this class, or the correlation matrix of an initializedobject of this class. In most situations, onlycorBlomberg willbe called by the user.

Author(s)

Emmanuel Paradis

References

Blomberg, S. P., Garland, Jr, T., and Ives, A. R. (2003) Testing forphylogenetic signal in comparative data: behavioral traits are morelabile.Evolution,57, 717–745.

Brownian Correlation Structure

Description

Expected covariance under a Brownian model (Felsenstein 1985,Martinsand Hansen 1997)

V_{ij} = \gamma \times t_a

wheret_a is the distance on the phylogeny between the rootand the most recent common ancestor of taxai andjand\gamma is a constant.

Usage

corBrownian(value=1, phy, form=~1)## S3 method for class 'corBrownian'coef(object, unconstrained = TRUE, ...)## S3 method for class 'corBrownian'corMatrix(object, covariate = getCovariate(object), corr = TRUE, ...)

Arguments

Value

An object of classcorBrownian, or the coefficient from anobject of this class (actually sendsnumeric(0)), or thecorrelation matrix of an initialized object of this class.

Author(s)

Julien Dutheildutheil@evolbio.mpg.de

References

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

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

See Also

corClasses

Phylogenetic Correlation Structures

Description

Classes of phylogenetic correlation structures ("corPhyl")available inape.

• corBrownian: Brownian motion model (Felsenstein 1985)

• corMartins: The covariance matrix defined in Martins and Hansen(1997)

• corGrafen: The covariance matrix defined in Grafen (1989)

• corPagel: The covariance matrix defined in Freckelton et al. (2002)

• corBlomberg: The covariance matrix defined in Blomberg et al. (2003)

See the help page of each class for references and detaileddescription.

Author(s)

Julien Dutheildutheil@evolbio.mpg.de, EmmanuelParadis

See Also

corClasses andgls in thenlme librarie,corBrownian,corMartins,corGrafen,corPagel,corBlomberg,vcv,vcv2phylo

Examples

library(nlme)txt <- "((((Homo:0.21,Pongo:0.21):0.28,Macaca:0.49):0.13,Ateles:0.62):0.38,Galago:1.00);"tree.primates <- read.tree(text = txt)X <- c(4.09434, 3.61092, 2.37024, 2.02815, -1.46968)Y <- c(4.74493, 3.33220, 3.36730, 2.89037, 2.30259)Species <- c("Homo", "Pongo", "Macaca", "Ateles", "Galago")dat <- data.frame(Species = Species, X = X, Y = Y)m1 <- gls(Y ~ X, dat, correlation=corBrownian(1, tree.primates, form = ~Species))summary(m1)m2 <- gls(Y ~ X, dat, correlation=corMartins(1, tree.primates, form = ~Species))summary(m2)corMatrix(m2$modelStruct$corStruct)m3 <- gls(Y ~ X, dat, correlation=corGrafen(1, tree.primates, form = ~Species))summary(m3)corMatrix(m3$modelStruct$corStruct)

Grafen's (1989) Correlation Structure

Description

Grafen's (1989) covariance structure. Branch lengths are computed usingGrafen's method (seecompute.brlen). The covariancematrice is then the traditional variance-covariance matrix for aphylogeny.

Usage

corGrafen(value, phy, form=~1, fixed = FALSE)## S3 method for class 'corGrafen'coef(object, unconstrained = TRUE, ...)## S3 method for class 'corGrafen'corMatrix(object,                  covariate = getCovariate(object), corr = TRUE, ...)

Arguments

Value

An object of classcorGrafen or the rho coefficient from anobject of this class or the correlation matrix of an initializedobject of this class.

Author(s)

Julien Dutheildutheil@evolbio.mpg.de

References

Grafen, A. (1989) The phylogenetic regression.PhilosophicalTransactions of the Royal society of London. Series B. BiologicalSciences,326, 119–157.

See Also

corClasses,compute.brlen,vcv.phylo.

Martins's (1997) Correlation Structure

Description

Martins and Hansen's (1997) covariance structure:

V_{ij} = \gamma \times e^{-\alpha t_{ij}}

wheret_{ij} is the phylogenetic distance between taxai andj and\gamma is a constant.

Usage

corMartins(value, phy, form = ~1, fixed = FALSE)## S3 method for class 'corMartins'coef(object, unconstrained = TRUE, ...)## S3 method for class 'corMartins'corMatrix(object,covariate = getCovariate(object), corr = TRUE, ...)

Arguments

Value

An object of classcorMartins or the alpha coefficient from an object of this classor the correlation matrix of an initialized object of this class.

Author(s)

Julien Dutheildutheil@evolbio.mpg.de

References

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

See Also

corClasses.

Pagel's “lambda” Correlation Structure

Description

The correlation structure from the present model is derived from theBrownian motion model by multiplying the off-diagonal elements (i.e.,the covariances) by\lambda. The variances are thus thesame than for a Brownian motion model.

Usage

corPagel(value, phy, form = ~1, fixed = FALSE)## S3 method for class 'corPagel'corMatrix(object, covariate = getCovariate(object),                   corr = TRUE, ...)## S3 method for class 'corPagel'coef(object, unconstrained = TRUE, ...)

Arguments

Value

an object of class"corPagel", the coefficients from an objectof this class, or the correlation matrix of an initialized object ofthis class. In most situations, onlycorPagel will be calledby the user.

Author(s)

Emmanuel Paradis

References

Freckleton, R. P., Harvey, P. H. and M. Pagel, M. (2002) Phylogeneticanalysis and comparative data: a test and review of evidence.American Naturalist,160, 712–726.

Pagel, M. (1999) Inferring the historical patterns of biologicalevolution.Nature,401,877–884.

Correlations among Multiple Traits with Phylogenetic Signal

Description

This function calculates Pearson correlation coefficients for multiple continuous traits that may have phylogenetic signal, allowing users to specify measurement error as the standard error of trait values at the tips of the phylogenetic tree. Phylogenetic signal for each trait is estimated from the data assuming that trait evolution is given by a Ornstein-Uhlenbeck process. Thus, the function allows the estimation of phylogenetic signal in multiple traits while incorporating correlations among traits. It is also possible to include independent variables (covariates) for each trait to remove possible confounding effects. corphylo() returns the correlation matrix for trait values, estimates of phylogenetic signal for each trait, and regression coefficients for independent variables affecting each trait.

Usage

corphylo(X, U = list(), SeM = NULL, phy = NULL, REML = TRUE,method = c("Nelder-Mead", "SANN"), constrain.d = FALSE, reltol = 10^-6,maxit.NM = 1000, maxit.SA = 1000, temp.SA = 1, tmax.SA = 1, verbose = FALSE)## S3 method for class 'corphylo'print(x, digits = max(3, getOption("digits") - 3), ...)

Arguments

Details

For the case of two variables, the function estimates parameters for the model of the form, for example,

X[1] = B[1,0] + B[1,1] * u[1,1] + \epsilon[1]

X[2] = B[2,0] + B[2,1] * u[2,1] + \epsilon[2]

\epsilon ~ Gaussian(0, V)

whereB[1,0],B[1,1],B[2,0], andB[2,1] are regression coefficients, andV is a variance-covariance matrix containing the correlation coefficient r, parameters of the OU processd1 andd2, and diagonal matricesM1 andM2 of measurement standard errors forX[1] andX[2]. The matrixV is2n x 2n, withn x n blocks given by

V[1,1] = C[1,1](d1) + M1

V[1,2] = C[1,2](d1,d2)

V[2,1] = C[2,1](d1,d2)

V[2,2] = C[2,2](d2) + M2

whereC[i,j](d1,d2) are derived from phy under the assumption of joint OU evolutionary processes for each trait (see Zheng et al. 2009). This formulation extends in the obvious way to more than two traits.

Value

An object of class "corphylo".

Author(s)

Anthony R. Ives

References

Zheng, L., A. R. Ives, T. Garland, B. R. Larget, Y. Yu, and K. F. Cao. 2009. New multivariate tests for phylogenetic signal and trait correlations applied to ecophysiological phenotypes of nineManglietia species.Functional Ecology23:1059–1069.

Examples

## Simple example using data without correlations or phylogenetic## signal. This illustrates the structure of the input data.phy <- rcoal(10, tip.label = 1:10)X <- matrix(rnorm(20), nrow = 10, ncol = 2)rownames(X) <- phy$tip.labelU <- list(NULL, matrix(rnorm(10, mean = 10, sd = 4), nrow = 10, ncol = 1))rownames(U[[2]]) <- phy$tip.labelSeM <- matrix(c(0.2, 0.4), nrow = 10, ncol = 2)rownames(SeM) <- phy$tip.labelcorphylo(X = X, SeM = SeM, U = U, phy = phy, method = "Nelder-Mead")## Not run: ## Simulation example for the correlation between two variables. The## example compares the estimates of the correlation coefficients from## corphylo when measurement error is incorporated into the analyses with## three other cases: (i) when measurement error is excluded, (ii) when## phylogenetic signal is ignored (assuming a "star" phylogeny), and (iii)## neither measurement error nor phylogenetic signal are included.## In the simulations, variable 2 is associated with a single## independent variable. This requires setting up a list U that has 2## elements: element U[[1]] is NULL and element U[[2]] is a n x 1 vector## containing simulated values of the independent variable.# Set up parameter values for simulating datan <- 50phy <- rcoal(n, tip.label = 1:n)R <- matrix(c(1, 0.7, 0.7, 1), nrow = 2, ncol = 2)d <- c(0.3, .95)B2 <- 1Se <- c(0.2, 1)SeM <- matrix(Se, nrow = n, ncol = 2, byrow = T)rownames(SeM) <- phy$tip.label# Set up needed matrices for the simulationsp <- length(d)star <- stree(n)star$edge.length <- array(1, dim = c(n, 1))star$tip.label <- phy$tip.labelVphy <- vcv(phy)Vphy <- Vphy/max(Vphy)Vphy <- Vphy/exp(determinant(Vphy)$modulus[1]/n)tau <- matrix(1, nrow = n, ncol = 1) C <- matrix(0, nrow = p * n, ncol = p * n)for (i in 1:p) for (j in 1:p) {Cd <- (d[i]^tau * (d[j]^t(tau)) * (1 - (d[i] * d[j])^Vphy))/(1 - d[i] * d[j])C[(n * (i - 1) + 1):(i * n), (n * (j - 1) + 1):(j * n)] <- R[i, j] * Cd}MM <- matrix(SeM^2, ncol = 1)V <- C + diag(as.numeric(MM))## Perform a Cholesky decomposition of Vphy. This is used to generate## phylogenetic signal: a vector of independent normal random variables,## when multiplied by the transpose of the Cholesky deposition of Vphy will## have covariance matrix equal to Vphy.iD <- t(chol(V))# Perform Nrep simulations and collect the resultsNrep <- 100cor.list <- matrix(0, nrow = Nrep, ncol = 1)cor.noM.list <- matrix(0, nrow = Nrep, ncol = 1)cor.noP.list <- matrix(0, nrow = Nrep, ncol = 1)cor.noMP.list <- matrix(0, nrow = Nrep, ncol = 1)d.list <- matrix(0, nrow = Nrep, ncol = 2)d.noM.list <- matrix(0, nrow = Nrep, ncol = 2)B.list <- matrix(0, nrow = Nrep, ncol = 3)B.noM.list <- matrix(0, nrow = Nrep, ncol = 3)B.noP.list <- matrix(0, nrow = Nrep, ncol = 3)for (rep in 1:Nrep) {XX <- iD X <- matrix(XX, nrow = n, ncol = 2)rownames(X) <- phy$tip.labelU <- list(NULL, matrix(rnorm(n, mean = 2, sd = 10), nrow = n, ncol = 1))rownames(U[[2]]) <- phy$tip.labelcolnames(U[[2]]) <- "V1"X[,2] <- X[,2] + B2[1] * U[[2]][,1] - B2[1] * mean(U[[2]][,1])z <- corphylo(X = X, SeM = SeM, U = U, phy = phy, method = "Nelder-Mead")z.noM <- corphylo(X = X, U = U, phy = phy, method = "Nelder-Mead")z.noP <- corphylo(X = X, SeM = SeM, U = U, phy = star, method = "Nelder-Mead")cor.list[rep] <- z$cor.matrix[1, 2]cor.noM.list[rep] <- z.noM$cor.matrix[1, 2]cor.noP.list[rep] <- z.noP$cor.matrix[1, 2]cor.noMP.list[rep] <- cor(cbind(lm(X[,1] ~ 1)$residuals, lm(X[,2] ~ U[[2]])$residuals))[1,2]d.list[rep, ] <- z$dd.noM.list[rep, ] <- z.noM$dB.list[rep, ] <- z$BB.noM.list[rep, ] <- z.noM$BB.noP.list[rep, ] <- z.noP$Bshow(c(rep, z$convcode, z$cor.matrix[1, 2], z$d))}correlation <- rbind(R[1, 2], mean(cor.list), mean(cor.noM.list),                     mean(cor.noP.list), mean(cor.noMP.list))rownames(correlation) <- c("True", "With SeM and Phy", "Without SeM",                           "Without Phy", "Without Phy or SeM")correlationsignal.d <- rbind(d, colMeans(d.list), colMeans(d.noM.list))rownames(signal.d) <- c("True", "With SeM and Phy", "Without SeM")signal.dest.B <- rbind(c(0, 0, B2), colMeans(B.list), colMeans(B.noM.list),               colMeans(B.noP.list))rownames(est.B) <- c("True", "With SeM and Phy", "Without SeM", "Without Phy")colnames(est.B) <- rownames(z$B)est.B# Example simulation output# correlation                        # [,1]# True               0.7000000# With SeM and Phy   0.7055958# Without SeM        0.3125253# Without Phy        0.4054043# Without Phy or SeM 0.3476589# signal.d                     # [,1]      [,2]# True             0.300000 0.9500000# With SeM and Phy 0.301513 0.9276663# Without SeM      0.241319 0.4872675# est.B                        # B1.0      B2.0     B2.V1# True              0.00000000 0.0000000 1.0000000# With SeM and Phy -0.01285834 0.2807215 0.9963163# Without SeM       0.01406953 0.3059110 0.9977796# Without Phy       0.02139281 0.3165731 0.9942140## End(Not run)

Phylogenetic Correlogram

Description

This function computes a correlogram from taxonomic levels.

Usage

  correlogram.formula(formula, data = NULL, use = "all.obs")

Arguments

Details

See the vignette in R:vignette("MoranI").

Value

An object of classcorrelogram which is a data frame with threecolumns:

or an object of classcorrelogramList containing a list ofobjects of classcorrelogram if several variables are given asresponse informula.

Author(s)

Julien Dutheildutheil@evolbio.mpg.de andEmmanuel Paradis

See Also

plot.correlogram,Moran.I

Examples

data(carnivora)### Using the formula interface:co <- correlogram.formula(SW ~ Order/SuperFamily/Family/Genus,      data=carnivora)coplot(co)### Several correlograms on the same plot:cos <- correlogram.formula(SW + FW ~ Order/SuperFamily/Family/Genus,      data=carnivora)cosplot(cos)

NEXUS Data Example

Description

Example of Protein data in NEXUS format (Maddison et al., 1997).Data is written in interleaved format using a single DATA block.Original data from Rokas et al (2002).

Usage

data(cynipids)

Format

ASCII text in NEXUS format

References

Maddison, D. R., Swofford, D. L. and Maddison, W. P. (1997) NEXUS: anextensible file format for systematic information.SystematicBiology,46, 590–621.

Rokas, A., Nylander, J. A. A., Ronquist, F. and Stone, G. N. (2002) Amaximum likelihood analysis of eight phylogenetic markers in Gallwasps(Hymenoptera: Cynipidae): implications for insect phylogeneticstudies.Molecular Phylogenetics and Evolution,22,206–219.

Probability Density Under Birth–Death Models

Description

These functions compute the probability density under somebirth–death models, that is the probability of obtainingxspecies after a timet giving how speciation and extinctionprobabilities vary through time (these may be constant, or even equalto zero for extinction).

Usage

dyule(x, lambda = 0.1, t = 1, log = FALSE)dbd(x, lambda, mu, t, conditional = FALSE, log = FALSE)dbdTime(x, birth, death, t, conditional = FALSE,        BIRTH = NULL, DEATH = NULL, fast = FALSE)

Arguments

Details

These three functions compute the probabilities to observexspecies starting from a single one after timet (assumed to becontinuous). The first function is a short-cut for the second one withmu = 0 and with default values for the two other arguments.dbdTime is for time-varyinglambda andmuspecified asR functions.

dyule is vectorized simultaneously on its three argumentsx,lambda, andt, according toR's rules ofrecycling arguments.dbd is vectorized simultaneouslyxandt (to make likelihood calculations easy), anddbdTime is vectorized only onx; the other arguments areeventually shortened with a warning if necessary.

The returned value is, logically, zero for values ofx out ofrange, i.e., negative or zero fordyule or ifconditional = TRUE. However, it is not checked if the values ofx arepositive non-integers and the probabilities are computed and returned.

The details on the form of the argumentsbirth,death,BIRTH,DEATH, andfast can be found in the linksbelow.

Value

a numeric vector.

Note

If you use these functions to calculate a likelihood function, it isstrongly recommended to compute the log-likelihood with, for instancein the case of a Yule process,sum(dyule( , log = TRUE)) (seeexamples).

Author(s)

Emmanuel Paradis

References

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

See Also

bd.time,yule.time

Examples

x <- 0:10plot(x, dyule(x), type = "h", main = "Density of the Yule process")text(7, 0.85, expression(list(lambda == 0.1, t == 1)))y <- dbd(x, 0.1, 0.05, 10)z <- dbd(x, 0.1, 0.05, 10, conditional = TRUE)d <- rbind(y, z)colnames(d) <- xbarplot(d, beside = TRUE, ylab = "Density", xlab = "Number of species",        legend = c("unconditional", "conditional on\nno extinction"),        args.legend = list(bty = "n"))title("Density of the birth-death process")text(17, 0.4, expression(list(lambda == 0.1, mu == 0.05, t == 10)))## Not run: ### generate 1000 values from a Yule process with lambda = 0.05x <- replicate(1e3, Ntip(rlineage(0.05, 0)))### the correct way to calculate the log-likelihood...:sum(dyule(x, 0.05, 50, log = TRUE))### ... and the wrong way:log(prod(dyule(x, 0.05, 50)))### a third, less preferred, way:sum(log(dyule(x, 0.05, 50)))## End(Not run)

Definition of Vectors for Plotting or Annotating

Description

This function can be used to define vectors to annotate a set of taxonnames, labels, etc. It should facilitate the (re)definition of coloursor similar attributes for plotting trees or other graphics.

Usage

def(x, ..., default = NULL, regexp = FALSE)

Arguments

Details

The idea of this function is to make the definition of colours, etc.,simpler than what is done usually. A typical use is:

def(tr$tip.label, Homo_sapiens = "blue")

which will return a vector of character strings all "black" except onematching the tip label "Homo_sapiens" which will be "blue". Another usecould be:

def(tr$tip.label, Homo_sapiens = 2)

which will return a vector a numerical values all 1 except for"Homo_sapiens" which will be 2. Several definitions can be done, e.g.:

def(tr$tip.label, Homo_sapiens = "blue", Pan_paniscus = "red")

The default value is determined with respect to the mode of the valuesgiven with the... (either "black" or 1).

Ifregexp = TRUE is used, then the names of the statements must bequoted, e.g.:

def(tr$tip.label, "^Pan_" = "red", regexp = TRUE)

will return "red" for all labels starting with "Pan_".

Value

a vector of the same length thanx.

Author(s)

Emmanuel Paradis

Examples

data(bird.orders)a <- def(bird.orders$tip.label, Galliformes = 2)str(a) # numericplot(bird.orders, font = a)co <- def(bird.orders$tip.label, Passeriformes = "red", Trogoniformes = "blue")str(co) # characterplot(bird.orders, tip.color = co)### use of a regexp (so we need to quote it) to colour all orders### with names starting with "C" (and change the default):co2 <- def(bird.orders$tip.label, "^C" = "gold", default = "grey", regexp = TRUE)plot(bird.orders, tip.color = co2)

Vertex Degrees in Trees and Networks

Description

degree is a generic function to calculate the degree of allnodes in a tree or in a network.

Usage

degree(x, ...)## S3 method for class 'phylo'degree(x, details = FALSE, ...)## S3 method for class 'evonet'degree(x, details = FALSE, ...)

Arguments

Details

The degree of a node (or vertex) in a network is defined by the numberof branches (or edges) that connect to this node. In a phylogenetictree, the tips (or terminal nodes) are of degree one, and the(internal) nodes are of degree two or more.

There are currently two methods for the classes"phylo" and"evonet". The default of these functions is to return a summarytable with the degrees observed in the tree or network in the firstcolumn, and the number of nodes in the second column. Ifdetails = TRUE, a vector giving the degree of each node (as numbered in theedge matrix) is returned.

The validity of the object is not checked, sodegree can beused to check problems with badly conformed trees.

Value

a data frame ifdetails = FALSE, or a vector of integersotherwise.

Author(s)

Emmanuel Paradis

See Also

checkValidPhylo

Examples

data(bird.orders)degree(bird.orders)degree(bird.orders, details = TRUE)data(bird.families)degree(bird.families)degree(rtree(10)) # 10, 1, 8degree(rtree(10, rooted = FALSE)) # 10, 0, 8degree(stree(10)) # 10 + 1 node of degree 10

Delete Alignment Gaps in DNA or AA Sequences

Description

These functions remove gaps ("-") in a sample of DNA sequences.

Usage

del.gaps(x)del.colgapsonly(x, threshold = 1, freq.only = FALSE)del.rowgapsonly(x, threshold = 1, freq.only = FALSE)

Arguments

Details

del.gaps remove all gaps, so the returned sequences may nothave all the same lengths and are therefore returned in a list.

del.colgapsonly removes the columns with a proportion at leastthreshold of gaps. Thus by default, only the columns with gapsonly are removed (useful when a small matrix is extracted from a largealignment).del.rowgapsonly does the same for the rows.

The class of the input sequences is respected and kept unchanged,unless it contains neither"DNAbin" nor"AAbin" in whichcase the object is first converted into the class"DNAbin".

Value

del.gaps returns a vector (if there is only one input sequence)or a list of sequences;del.colgapsonly anddel.rowgapsonly return a matrix of sequences or a numericvector (with names for the second function) iffreq.only = TRUE.

Author(s)

Emmanuel Paradis

See Also

base.freq,seg.sites,image.DNAbin,checkAlignment

Delta Plots

Description

This function makes a\delta plot following Holland etal. (2002).

Usage

delta.plot(X, k = 20, plot = TRUE, which = 1:2)

Arguments

Details

See Holland et al. (2002) for details and interpretation.

The computing time of this function is proportional to the fourthpower of the number of observations (O(n^4)), so calculationsmay be very long with only a slight increase in sample size.

Value

This function returns invisibly a named list with two components:

• counts: the counts for the histogram of\delta_q values

• delta.bar: the mean\delta value for eachobservation

Author(s)

Emmanuel Paradis

References

Holland, B. R., Huber, K. T., Dress, A. and Moulton, V. (2002) Deltaplots: a tool for analyzing phylogenetic distance data.Molecular Biology and Evolution,12, 2051–2059.

See Also

dist.dna

Examples

data(woodmouse)d <- dist.dna(woodmouse)delta.plot(d)layout(1)delta.plot(d, 40, which = 1)

Pairwise Distances from DNA Sequences

Description

This function computes a matrix of pairwise distances from DNAsequences using a model of DNA evolution. Eleven substitution models(and the raw distance) are currently available.

Usage

dist.dna(x, model = "K80", variance = FALSE,         gamma = FALSE, pairwise.deletion = FALSE,         base.freq = NULL, as.matrix = FALSE)

Arguments

Details

The molecular evolutionary models available through the optionmodel have been extensively described in the literature. Abrief description is given below; more details can be found in thereferences.

• raw,N: This is simply the proportion or the number ofsites that differ between each pair of sequences. This may be usefulto draw “saturation plots”. The optionsvariance andgamma have no effect, butpairwise.deletion may have.

• TS,TV: These are the numbers of transitions andtransversions, respectively.

• JC69: This model was developed by Jukes and Cantor (1969). Itassumes that all substitutions (i.e. a change of a base by anotherone) have the same probability. This probability is the same for allsites along the DNA sequence. This last assumption can be relaxed byassuming that the substition rate varies among site following agamma distribution which parameter must be given by the user. Bydefault, no gamma correction is applied. Another assumption is thatthe base frequencies are balanced and thus equal to 0.25.

• K80: The distance derived by Kimura (1980), sometimes referredto as “Kimura's 2-parameters distance”, has the same underlyingassumptions than the Jukes–Cantor distance except that two kinds ofsubstitutions are considered: transitions (A <-> G, C <-> T), andtransversions (A <-> C, A <-> T, C <-> G, G <-> T). They are assumedto have different probabilities. A transition is the substitution ofa purine (C, T) by another one, or the substitution of a pyrimidine(A, G) by another one. A transversion is the substitution of apurine by a pyrimidine, or vice-versa. Both transition andtransversion rates are the same for all sites along the DNAsequence. Jin and Nei (1990) modified the Kimura model to allow forvariation among sites following a gamma distribution. Like for theJukes–Cantor model, the gamma parameter must be given by theuser. By default, no gamma correction is applied.

• F81: Felsenstein (1981) generalized the Jukes–Cantor modelby relaxing the assumption of equal base frequencies. The formulaeused in this function were taken from McGuire et al. (1999).

• K81: Kimura (1981) generalized his model (Kimura 1980) byassuming different rates for two kinds of transversions: A <-> C andG <-> T on one side, and A <-> T and C <-> G on the other. This iswhat Kimura called his “three substitution types model” (3ST), andis sometimes referred to as “Kimura's 3-parameters distance”.

• F84: This model generalizes K80 by relaxing the assumptionof equal base frequencies. It was first introduced by Felsenstein in1984 in Phylip, and is fully described by Felsenstein and Churchill(1996). The formulae used in this function were taken from McGuireet al. (1999).

• BH87: Barry and Hartigan (1987) developed a distance basedon the observed proportions of changes among the four bases. Thisdistance is not symmetric.

• T92: Tamura (1992) generalized the Kimura model by relaxingthe assumption of equal base frequencies. This is done by takinginto account the bias in G+C content in the sequences. Thesubstitution rates are assumed to be the same for all sites alongthe DNA sequence.

• TN93: Tamura and Nei (1993) developed a model which assumesdistinct rates for both kinds of transition (A <-> G versus C <->T), and transversions. The base frequencies are not assumed to beequal and are estimated from the data. A gamma correction of theinter-site variation in substitution rates is possible.

• GG95: Galtier and Gouy (1995) introduced a model where theG+C content may change through time. Different rates are assumed fortransitons and transversions.

• logdet: The Log-Det distance, developed by Lockhart etal. (1994), is related to BH87. However, this distance issymmetric. Formulae from Gu and Li (1996) are used.dist.logdet inphangorn uses a differentimplementation that gives substantially different distances forlow-diverging sequences.

• paralin: Lake (1994) developed the paralinear distance whichcan be viewed as another variant of the Barry–Hartigan distance.

• indel: this counts the number of sites where there is aninsertion/deletion gap in one sequence and not in the other.

• indelblock: same than before but contiguous gaps arecounted as a single unit. Note that the distance between-A- andA-- is 3 because there are three different blocks of gaps, whereasthe “indel” distance will be 2.

Value

an object of classdist (by default), or a numericmatrix ifas.matrix = TRUE. Ifmodel = "BH87", a numericmatrix is returned because the Barry–Hartigan distance is notsymmetric.

Ifvariance = TRUE an attribute called"variance" isgiven to the returned object.

Note

If the sequences are very different, most evolutionary distances areundefined and a non-finite value (Inf or NaN) is returned. You may dodist.dna(, model = "raw") to check whether some values arehigher than 0.75.

Author(s)

Emmanuel Paradis

References

Barry, D. and Hartigan, J. A. (1987) Asynchronous distance betweenhomologous DNA sequences.Biometrics,43, 261–276.

Felsenstein, J. (1981) Evolutionary trees from DNA sequences: amaximum likelihood approach.Journal of Molecular Evolution,17, 368–376.

Felsenstein, J. and Churchill, G. A. (1996) A Hidden Markov modelapproach to variation among sites in rate of evolution.Molecular Biology and Evolution,13, 93–104.

Galtier, N. and Gouy, M. (1995) Inferring phylogenies from DNAsequences of unequal base compositions.Proceedings of theNational Academy of Sciences USA,92, 11317–11321.

Gu, X. and Li, W.-H. (1996) Bias-corrected paralinear and LogDetdistances and tests of molecular clocks and phylogenies undernonstationary nucleotide frequencies.Molecular Biology andEvolution,13, 1375–1383.

Jukes, T. H. and Cantor, C. R. (1969) Evolution of proteinmolecules. inMammalian Protein Metabolism, ed. Munro, H. N.,pp. 21–132, New York: Academic Press.

Kimura, M. (1980) A simple method for estimating evolutionary rates ofbase substitutions through comparative studies of nucleotidesequences.Journal of Molecular Evolution,16, 111–120.

Kimura, M. (1981) Estimation of evolutionary distances betweenhomologous nucleotide sequences.Proceedings of the NationalAcademy of Sciences USA,78, 454–458.

Jin, L. and Nei, M. (1990) Limitations of the evolutionary parsimonymethod of phylogenetic analysis.Molecular Biology andEvolution,7, 82–102.

Lake, J. A. (1994) Reconstructing evolutionary trees from DNA andprotein sequences: paralinear distances.Proceedings of theNational Academy of Sciences USA,91, 1455–1459.

Lockhart, P. J., Steel, M. A., Hendy, M. D. and Penny, D. (1994)Recovering evolutionary trees under a more realistic model of sequenceevolution.Molecular Biology and Evolution,11,605–602.

McGuire, G., Prentice, M. J. and Wright, F. (1999). Improved errorbounds for genetic distances from DNA sequences.Biometrics,55, 1064–1070.

Tamura, K. (1992) Estimation of the number of nucleotide substitutionswhen there are strong transition-transversion and G + C-contentbiases.Molecular Biology and Evolution,9, 678–687.

Tamura, K. and Nei, M. (1993) Estimation of the number of nucleotidesubstitutions in the control region of mitochondrial DNA in humans andchimpanzees.Molecular Biology and Evolution,10, 512–526.

See Also

read.GenBank,read.dna,write.dna,DNAbin,dist.gene,cophenetic.phylo,dist

Pairwise Distances from Genetic Data

Description

This function computes a matrix of distances between pairs ofindividuals from a matrix or a data frame of genetic data.

Usage

dist.gene(x, method = "pairwise", pairwise.deletion = FALSE,          variance = FALSE)

Arguments

Details

This function is meant to be very general and accepts different kindsof data (alleles, haplotypes, SNP, DNA sequences, ...). The rows ofthe data matrix represent the individuals, and the columns the loci.

In the case of the pairwise method, the distanced between twoindividuals is the number of loci for which they differ, and theassociated variance isd(L - d)/L, whereL is the numberof loci.

In the case of the percentage method, this distance is divided byL,and the associated variance isd(1 - d)/L.

For more elaborate distances with DNA sequences, see the functiondist.dna.

Value

an object of classdist. Ifvariance = TRUE anattribute called"variance" is given to the returned object.

Note

Missing data (NA) are coded and treated in R's usual way.

Author(s)

Emmanuel Paradis

See Also

dist.dna,cophenetic.phylo,dist

Topological Distances Between Two Trees

Description

This function computes the topological distance between twophylogenetic trees or among trees in a list (ify = NULL usingdifferent methods.

Usage

dist.topo(x, y = NULL, method = "PH85", mc.cores = 1)

Arguments

Details

Two methods are available: the one by Penny and Hendy (1985,originally from Robinson and Foulds 1981), and the branch length scoreby Kuhner and Felsenstein (1994). The trees are always considered asunrooted.

The topological distance is defined as twice the number of internalbranches defining different bipartitions of the tips (Robinson andFoulds 1981; Penny and Hendy 1985). Rzhetsky and Nei (1992) proposed amodification of the original formula to take multifurcations intoaccount.

The branch length score may be seen as similar to the previousdistance but taking branch lengths into account. Kuhner andFelsenstein (1994) proposed to calculate the square root of the sum ofthe squared differences of the (internal) branch lengths definingsimilar bipartitions (or splits) in both trees.

Value

a single numeric value if bothx andy are used, anobject of class"dist" otherwise.

Note

The geodesic distance of Billera et al. (2001) has been disabled: seethe packagedistory on CRAN.

Author(s)

Emmanuel Paradis

References

Billera, L. J., Holmes, S. P. and Vogtmann, K. (2001) Geometry of thespace of phylogenetic trees.Advances in Applied Mathematics,27, 733–767.

Kuhner, M. K. and Felsenstein, J. (1994) Simulation comparison ofphylogeny algorithms under equal and unequal evolutionary rates.Molecular Biology and Evolution,11, 459–468.

Nei, M. and Kumar, S. (2000)Molecular Evolution andPhylogenetics. Oxford: Oxford University Press.

Penny, D. and Hendy, M. D. (1985) The use of tree comparisonmetrics.Systemetic Zoology,34, 75–82.

Robinson, D. F. and Foulds, L. R. (1981) Comparison of phylogenetictrees.Mathematical Biosciences,53, 131–147.

Rzhetsky, A. and Nei, M. (1992) A simple method for estimating andtesting minimum-evolution trees.Molecular Biology andEvolution,9, 945–967.

See Also

cophenetic.phylo,prop.part

Examples

ta <- rtree(30, rooted = FALSE)tb <- rtree(30, rooted = FALSE)dist.topo(ta, ta) # 0dist.topo(ta, tb) # unlikely to be 0## rmtopology() simulated unrooted trees by default:TR <- rmtopology(100, 10)## these trees have 7 internal branches, so the maximum distance## between two of them is 14:DTR <- dist.topo(TR)table(DTR)

Tests of Constant Diversification Rates

Description

This function computes two tests of the distribution of branchingtimes using the Cramér–von Mises and Anderson–Darlinggoodness-of-fit tests. By default, it is assumed that thediversification rate is constant, and an exponential distribution isassumed for the branching times. In this case, the expecteddistribution under this model is computed with a rate estimated fromthe data. Alternatively, the user may specify an expected cumulativedensity function (z): in this case,x andz mustbe of the same length. See the examples for how to compute the latterfrom a sample of expected branching times.

Usage

diversi.gof(x, null = "exponential", z = NULL)

Arguments

Details

The Cramér–von Mises and Anderson–Darling testscompare the empirical density function (EDF) of the observations to anexpected cumulative density function. By contrast to theKolmogorov–Smirnov test where the greatest difference between thesetwo functions is used, in both tests all differences are taken intoaccount.

The distributions of both test statistics depend on the nullhypothesis, and on whether or not some parameters were estimated fromthe data. However, these distributions are not known precisely andcritical values were determined by Stephens (1974) usingsimulations. These critical values were used for the present function.

Value

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

Author(s)

Emmanuel Paradis

References

Paradis, E. (1998) Testing for constant diversification rates usingmolecular phylogenies: a general approach based on statistical testsfor goodness of fit.Molecular Biology and Evolution,15, 476–479.

Stephens, M. A. (1974) EDF statistics for goodness of fit and somecomparisons.Journal of the American Statistical Association,69, 730–737.

See Also

branching.times,diversi.timeltt.plot,birthdeath,yule,yule.cov

Examples

data(bird.families)x <- branching.times(bird.families)### suppose we have a sample of expected branching times `y';### for simplicity, take them from a uniform distribution:y <- runif(500, 0, max(x) + 1) # + 1 to avoid A2 = Inf### now compute the expected cumulative distribution:x <- sort(x)N <- length(x)ecdf <- numeric(N)for (i in 1:N) ecdf[i] <- sum(y <= x[i])/500### finally do the test:diversi.gof(x, "user", z = ecdf)

Analysis of Diversification with Survival Models

Description

This functions fits survival models to a set of branching times, someof them may be known approximately (censored). Three models arefitted, Model A assuming constant diversification, Model B assumingthat diversification follows a Weibull law, and Model C assuming thatdiversification changes with a breakpoint at time ‘Tc’. The models arefitted by maximum likelihood.

Usage

diversi.time(x, census = NULL, censoring.codes = c(1, 0), Tc = NULL)

Arguments

Details

The principle of the method is to consider each branching time as anevent: if the branching time is accurately known, then it is a failureevent; if it is approximately knwon then it is a censoring event. Ananalogy is thus made between the failure (or hazard) rate estimated bythe survival models and the diversification rate of the lineage. Timeis here considered from present to past.

Model B assumes a monotonically changing diversification rate. Theparameter that controls the change of this rate is called beta. Ifbeta is greater than one, then the diversification rate decreasesthrough time; if it is lesser than one, the the rate increases throughtime. If beta is equal to one, then Model B reduces to Model A.

Value

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

Author(s)

Emmanuel Paradis

References

Paradis, E. (1997) Assessing temporal variations in diversificationrates from phylogenies: estimation and hypothesistesting.Proceedings of the Royal Society of London. SeriesB. Biological Sciences,264, 1141–1147.

See Also

branching.times,diversi.gofltt.plot,birthdeath,bd.ext,yule,yule.cov

Diversity Contrast Test

Description

This function performs the diversity contrast test comparing pairs ofsister-clades.

Usage

diversity.contrast.test(x, method = "ratiolog",        alternative = "two.sided", nrep = 0, ...)

Arguments

Details

Ifmethod = "ratiolog", the test described in Barraclough etal. (1996) is performed. Ifmethod = "proportion", the versionin Barraclough et al. (1995) is used. Ifmethod = "difference",the signed difference is used (Sargent 2004). Ifmethod = "logratio",then this is Wiegmann et al.'s (1993) version. Thesefour tests are essentially different versions of the same test (Vamosiand Vamosi 2005, Vamosi 2007). See Paradis (2012) for a comparison oftheir statistical performance with other tests.

Ifnrep = 0, a Wilcoxon test is done on the species diversitycontrasts with the null hypothesis is that they are distributed aroundzero. Ifnrep > 0, a randomization procedure is done where thesigns of the diversity contrasts are randomly chosen. This is used tocreate a distribution of the test statistic which is compared with theobserved value (the sum of the diversity contrasts).

Value

a single numeric value with theP-value.

Author(s)

Emmanuel Paradis

References

Barraclough, T. G., Harvey, P. H. and Nee, S. (1995) Sexualselection and taxonomic diversity in passerine birds.Proceedings of the Royal Society of London. Series B. BiologicalSciences,259, 211–215.

Barraclough, T. G., Harvey, P. H., and Nee, S. (1996) Rate ofrbcL gene sequence evolution and species diversification inflowering plants (angiosperms).Proceedings of the Royal Societyof London. Series B. Biological Sciences,263, 589–591.

Paradis, E. (2012) Shift in diversification in sister-cladecomparisons: a more powerful test.Evolution,66,288–295.

Sargent, R. D. (2004) Floral symmetry affects speciation rates inangiosperms.Proceedings of the Royal Society of London. SeriesB. Biological Sciences,271, 603–608.

Vamosi, S. M. (2007) Endless tests: guidelines for analysing non-nestedsister-group comparisons. An addendum.Evolutionary EcologyResearch,9, 717.

Vamosi, S. M. and Vamosi, J. C. (2005) Endless tests: guidelines foranalysing non-nested sister-group comparisons.EvolutionaryEcology Research,7, 567–579.

Wiegmann, B., Mitter, C. and Farrell, B. 1993. Diversification ofcarnivorous parasitic insects: extraordinary radiation or specializeddead end?American Naturalist,142, 737–754.

See Also

slowinskiguyer.test,mcconwaysims.testrichness.yule.test

Examples

### data from Vamosi & Vamosi (2005):fleshy <- c(1, 1, 1, 1, 1, 3, 3, 5, 9, 16, 33, 40, 50, 100, 216, 393, 850, 947,1700)dry <- c(2, 64, 300, 89, 67, 4, 34, 10, 150, 35, 2, 60, 81, 1, 3, 1, 11, 1, 18)x <- cbind(fleshy, dry)diversity.contrast.test(x)diversity.contrast.test(x, alt = "g")diversity.contrast.test(x, alt = "g", nrep = 1e4)slowinskiguyer.test(x)mcconwaysims.test(x)

dN/dS Ratio

Description

This function computes the pairwise ratios dN/dS for a set of alignedDNA sequences using Li's (1993) method.

Usage

dnds(x, code = 1, codonstart = 1, quiet = FALSE,     details = FALSE, return.categories = FALSE)

Arguments

Details

Sinceape 5.6, the degeneracy of each codon is calculateddirectly from the genetic code using the functiontrans. A consequence is that ambiguous bases are ignored(seesolveAmbiguousBases).

Ifdetails = TRUE, a table is printed for each pair ofsequences giving the numbers of transitions and transversions for eachcategory of degeneracy (nondegenerate, twofold, and fourfold). This ishelpful when non-meaningful values are returned (e.g., NaN, Inf,negative values).

Value

an object of class"dist", or a numeric matrix ifreturn.categories = TRUE.

Author(s)

Emmanuel Paradis

References

Li, W.-H. (1993) Unbiased estimation of the rates of synonymous andnonsynonymous substitution.Journal of Molecular Evolution,36, 96–99.

See Also

AAbin,trans,alview,solveAmbiguousBases

Examples

data(woodmouse)res <- dnds(woodmouse, quiet = TRUE) # NOT correctres2 <- dnds(woodmouse, code = 2, quiet = TRUE) # using the correct codeidentical(res, res2) # FALSE...cor(res, res2) # ... but very close## There a few N's in the woodmouse data, but this does not affect## greatly the results:res3 <- dnds(solveAmbiguousBases(woodmouse), code = 2, quiet = TRUE)cor(res, res3)## a simple example showing the usefulness of 'details = TRUE'X <- as.DNAbin(matrix(c("C", "A", "G", "G", "T", "T"), 2, 3))alview(X)dnds(X, quiet = TRUE) # NaNdnds(X, details = TRUE) # only a TV at a nondegenerate site

Remove Tips in a Phylogenetic Tree

Description

drop.tip removes the terminal branches of a phylogenetic tree,possibly removing the corresponding internal branches.keep.tipdoes the opposite operation (i.e., returns the induced tree).

extract.clade does the inverse operation: it keeps all the tipsfrom a given node, and deletes all the other tips.

Usage

drop.tip(phy, tip, ...)## S3 method for class 'phylo'drop.tip(phy, tip, trim.internal = TRUE, subtree = FALSE,         root.edge = 0, rooted = is.rooted(phy), collapse.singles = TRUE,         interactive = FALSE, ...)## S3 method for class 'multiPhylo'drop.tip(phy, tip, ...)keep.tip(phy, tip, ...)## S3 method for class 'phylo'keep.tip(phy, tip, ...)## S3 method for class 'multiPhylo'keep.tip(phy, tip, ...)extract.clade(phy, node, root.edge = 0, collapse.singles = TRUE,              interactive = FALSE)

Arguments

Details

The argumenttip can be either character or numeric. In thefirst case, it gives the labels of the tips to be deleted; in thesecond case the numbers of these labels in the vectorphy$tip.label are given.

This also applies tonode, but if this argument is characterand the tree has no node label, this results in an error. If more thanone value is given withnode (i.e., a vector of length two ormore), only the first one is used with a warning.

Iftrim.internal = FALSE, the new tips are given"NA" aslabels, unless there are node labels in the tree in which case theyare used.

Ifsubtree = TRUE, the returned tree has one or severalterminal branches named with node labels if available. Otherwise it isindicated how many tips have been removed (with a label"[x_tips]").This is done for as many monophyletic groups that have been deleted.

Note thatsubtree = TRUE impliestrim.internal = TRUE.

To undestand how the optionroot.edge works, see the examplesbelow. Ifrooted = FALSE and the tree has a root edge, thelatter is removed in the output.

Value

an object of class"phylo".

Author(s)

Emmanuel Paradis, Klaus Schliep, Joseph Brown

See Also

bind.tree,root

Examples

data(bird.families)tip <- c("Eopsaltriidae", "Acanthisittidae", "Pittidae", "Eurylaimidae","Philepittidae", "Tyrannidae", "Thamnophilidae", "Furnariidae","Formicariidae", "Conopophagidae", "Rhinocryptidae", "Climacteridae","Menuridae", "Ptilonorhynchidae", "Maluridae", "Meliphagidae","Pardalotidae", "Petroicidae", "Irenidae", "Orthonychidae","Pomatostomidae", "Laniidae", "Vireonidae", "Corvidae","Callaeatidae", "Picathartidae", "Bombycillidae", "Cinclidae","Muscicapidae", "Sturnidae", "Sittidae", "Certhiidae","Paridae", "Aegithalidae", "Hirundinidae", "Regulidae","Pycnonotidae", "Hypocoliidae", "Cisticolidae", "Zosteropidae","Sylviidae", "Alaudidae", "Nectariniidae", "Melanocharitidae","Paramythiidae","Passeridae", "Fringillidae")plot(drop.tip(bird.families, tip))plot(drop.tip(bird.families, tip, trim.internal = FALSE))data(bird.orders)plot(drop.tip(bird.orders, 6:23, subtree = TRUE))plot(drop.tip(bird.orders, c(1:5, 20:23), subtree = TRUE))plot(drop.tip(bird.orders, c(1:20, 23), subtree = TRUE))plot(drop.tip(bird.orders, c(1:20, 23), subtree = TRUE, rooted = FALSE))### Examples of the use of `root.edge'tr <- read.tree(text = "(A:1,(B:1,(C:1,(D:1,E:1):1):1):1):1;")drop.tip(tr, c("A", "B"), root.edge = 0) # = (C:1,(D:1,E:1):1);drop.tip(tr, c("A", "B"), root.edge = 1) # = (C:1,(D:1,E:1):1):1;drop.tip(tr, c("A", "B"), root.edge = 2) # = (C:1,(D:1,E:1):1):2;drop.tip(tr, c("A", "B"), root.edge = 3) # = (C:1,(D:1,E:1):1):3;

Draw Additional Edges on a Plotted Tree

Description

edges draws edges on a plotted tree.fancyarrowsenhancesarrows with triangle and harpoonheads; it can be called fromedges.

Usage

edges(nodes0, nodes1, arrows = 0, type = "classical", ...)fancyarrows(x0, y0, x1, y1, length = 0.25, angle = 30, code = 2,            col = par("fg"), lty = par("lty"), lwd = par("lwd"),            type = "triangle", ...)

Arguments

Details

The first function is helpful when drawing reticulations on a phylogeny,especially if computed from the edge matrix.

fancyarrows does not work with log-transformed scale(s).

Author(s)

Emmanuel Paradis

See Also

plot.phylo,nodelabels

Examples

set.seed(2)tr <- rcoal(6)plot(tr, "c")edges(10, 9, col = "red", lty = 2)edges(10:11, 8, col = c("blue", "green")) # recycling of 'nodes1'edges(1, 2, lwd = 2, type = "h", arrows = 3, col = "green")nodelabels()

Evolutionary Networks

Description

evonet builds a network from a tree of class"phylo". There areprint,plot, andreorder methods as well as a few conversion functions.

Usage

evonet(phy, from, to = NULL)## S3 method for class 'evonet'print(x, ...)## S3 method for class 'evonet'plot(x, col = "blue", lty = 1, lwd = 1, alpha = 0.5,              arrows = 0, arrow.type = "classical", ...)## S3 method for class 'evonet'Nedge(phy)## S3 method for class 'evonet'reorder(x, order = "cladewise", index.only = FALSE, ...)## S3 method for class 'evonet'as.phylo(x, ...)## S3 method for class 'evonet'as.networx(x, weight = NA, ...)## S3 method for class 'evonet'as.network(x, directed = TRUE, ...)## S3 method for class 'evonet'as.igraph(x, directed = TRUE, use.labels = TRUE, ...)as.evonet(x, ...)## S3 method for class 'phylo'as.evonet(x, ...)read.evonet(file = "", text = NULL, comment.char = "", ...)write.evonet(x, file = "", ...)

Arguments

Details

evonet is a constructor function that checks the arguments.

The classes"networx","network", and"igraph"are defined in the packagesphangorn,network, andigraph, respectively.

read.evonet reads networks from files in extended newick format(Cardona et al. 2008).

Value

an object of classc("evonet", "phylo") which is made of anobject of class"phylo" plus an elementreticulation coding additional edges among nodes and uses thesame coding rules than theedge matrix.

The conversion functions return an object of the appropriate class.

Author(s)

Emmanuel Paradis, Klaus Schliep

References

Cardona, G., Rossell, F., and Valiente, G. (2008) Extended Newick: itis time for a standard representation of phylogeneticnetworks.BMC Bioinformatics,9, 532.

See Also

as.networx in packagephangorn

Examples

tr <- rcoal(5)(x <- evonet(tr, 6:7, 8:9))plot(x)## simple example of extended Newick format:(enet <- read.evonet(text = "((a:2,(b:1)#H1:1):1,(#H1,c:1):2);"))plot(enet, arrows=1)## from Fig. 2 in Cardona et al. 2008:z <- read.evonet(text ="((1,((2,(3,(4)Y#H1)g)e,(((Y#H1, 5)h,6)f)X#H2)c)a,((X#H2,7)d,8)b)r;")zplot(z)## Not run: if (require(igraph)) {    plot(as.igraph(z))}## End(Not run)

Incomplete distances and edge weights of unrooted topology

Description

This function implements a method for checking whether an incompleteset of distances satisfy certain conditions that might make ituniquely determine the edge weights of a given topology, T. It printsinformation 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 alsochecks whether L is a triplet cover for T.

Usage

ewLasso(X, phy)

Arguments

Details

Missing values must be represented by eitherNA or a negative value.

This implements a method for checking whether an incomplete set ofdistances satisfies certain conditions that might make it uniquelydetermine the edge weights of a given topology, T. It printsinformation about whether the graph, G, with vertex set the set ofleaves, denoted by X, and edge set the set of non-missing distancepairs, denoted by L, is connected or strongly non-bipartite. It alsochecks whether L is a triplet cover for T. If G is not connected, thenT does not need to be the only topology satisfying the inputincomplete distances. If G is not strongly non-bipartite then theedge-weights of the edges of T are not the unique ones for which theinput distance is satisfied. If L is a triplet cover, then the inputdistance matrix uniquely determines the edge weights of T. See Dresset al. (2012) for details.

Value

NULL, the results are printed in the console.

Author(s)

Andrei Popescu

References

Dress, A. W. M., Huber, K. T., and Steel, M. (2012) ‘Lassoing’ aphylogentic tree I: basic properties, shellings and covers.Journal of Mathematical Biology,65(1), 77–105.

Gamma-Statistic of Pybus and Harvey

Description

This function computes the gamma-statistic which summarizes theinformation contained in the inter-node intervals of a phylogeny. Itis assumed that the tree is ultrametric. Note that the function doesnot check that the tree is effectively ultrametric, so if it is not,the returned result may not be meaningful.

Usage

gammaStat(phy)

Arguments

Details

The gamma-statistic is a summary of the information contained in theinter-node intervals of a phylogeny; it follows, under the assumptionthat the clade diversified with constant rates, a normal distributionwith mean zero and standard-deviation unity (Pybus and Harvey2000). Thus, the null hypothesis that the clade diversified withconstant rates may be tested with2*(1 - pnorm(abs(gammaStat(phy)))) for a two-tailed test, or1 - pnorm(abs(gammaStat(phy))) for a one-tailed test, both returningthe corresponding P-value.

Value

a numeric vector of length one.

Author(s)

Emmanuel Paradis

References

Pybus, O. G. and Harvey, P. H. (2000) Testing macro-evolutionarymodels using incomplete molecular phylogenies.Proceedings ofthe Royal Society of London. Series B. Biological Sciences,267, 2267–2272.

See Also

branching.times,ltt.plot,skyline

Read Annotations from GenBank

Description

This function connects to the GenBank database and reads sequence annotationsusing accession number(s) given as argument.

Usage

getAnnotationsGenBank(access.nb, quiet = TRUE)

Arguments

Details

The sequence annotations (a.k.a. feature list) are returned in a dataframe with five or six columns: start, end, type, product, others, andgene (the last being optional). This is the same information that canbe downloaded from NCBI's Web interface by clicking on ‘Send to:’,‘File’, and then selecting ‘Feature Table’ under ‘Format’.

A warning is given if some features are incomplete (this informationis then dropped from the returned object).

A warning is given if some accession numbers are not found on GenBank.

Value

One of the followings: (i) a data frame ifaccess.nb contains asingle accession number; (ii) a list of data frames ifaccess.nb contains several accession numbers, the names are setwithaccess.nb (if some accession numbers are not found onGenBank, the corresponding entries are set toNULL); (iii)NULL if all accession numbers are not found on GenBank.

Author(s)

Emmanuel Paradis

References

https://www.ncbi.nlm.nih.gov/Sequin/table.html (Note: itseems this URL is broken; 2022-01-03)

See Also

read.GenBank,read.gff,DNAbin

Examples

## The 8 sequences of tanagers (Ramphocelus):ref <- c("U15717", "U15718", "U15719", "U15720",         "U15721", "U15722", "U15723", "U15724")## Copy/paste or type the following commands if you## want to try them.## Not run: annot.rampho <- getAnnotationsGenBank(ref)annot.rampho## check all annotations are the same:unique(do.call(rbind, annot.rampho)[, -5])## End(Not run)

Phylogenetic Tree of 193 HIV-1 Sequences

Description

This data set describes an estimated clock-like phylogeny of 193 HIV-1group M sequences sampled in the Democratic Republic of Congo.

Usage

data(hivtree.newick)data(hivtree.table)

Format

hivtree.newick is a string with the tree in Newick format.The data framehivtree.table contains the corresponding internodedistances.

Source

This is a data example from Strimmer and Pybus (2001).

References

Strimmer, K. and Pybus, O. G. (2001) Exploring the demographic historyof DNA sequences using the generalized skyline plot.MolecularBiology and Evolution,18, 2298–2305.

See Also

coalescent.intervals,collapsed.intervals

Calculate Numbers of Phylogenetic Trees

Description

howmanytrees calculates the number of possible phylogenetictrees for a given number of tips.

LargeNumber is a utility function to compute (approximately)large numbers from the powera^b.

Usage

howmanytrees(n, rooted = TRUE, binary = TRUE,             labeled = TRUE, detail = FALSE)LargeNumber(a, b)## S3 method for class 'LargeNumber'print(x, latex = FALSE, digits = 1, ...)

Arguments

Details

In the cases of labeled binary trees, the calculation is done directlyand a single numeric value is returned (or an object of class"LargeNumber").

For multifurcating trees, and bifurcating, rooted, unlabeled trees,the calculation is done iteratively for 1 ton tips. Thus theuser can print all the intermediate values ifdetail = TRUE, oronly a single value ifdetail = FALSE (the default).

For multifurcating trees, ifdetail = TRUE, a matrix isreturned with the number of tips as rows (named from1 ton), and the number of nodes as columns (named from1 ton - 1). For bifurcating, rooted, unlabeled trees, a vector isreturned with names equal to the number of tips (from1 ton).

The number of unlabeled trees (aka tree shapes) can be computed onlyfor the rooted binary cases.

Note that if an infinite value (Inf) is returned this does notmean that there is an infinite number of trees (this cannot be if thenumber of tips is finite), but that the calculation is beyond thelimits of the computer. Only for the cases of rooted, binary, labeledtopologies an approximate number is returned in the form a"LargeNumber" object.

Value

a single numeric value, an object of class"LargeNumber", or inthe case wheredetail = TRUE is used, a named vector ormatrix.

Author(s)

Emmanuel Paradis

References

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

Examples

### Table 3.1 in Felsenstein 2004:for (i in c(1:20, 30, 40, 50))  cat(paste(i, howmanytrees(i), sep = "\t"), sep ="\n")### Table 3.6:howmanytrees(8, binary = FALSE, detail = TRUE)

Graphical Identification of Nodes and Tips

Description

This function allows to identify a clade on a plotted tree by clickingon the plot with the mouse. The tree, specified in the argumentx, must be plotted beforehand.

Usage

## S3 method for class 'phylo'identify(x, nodes = TRUE, tips = FALSE,                  labels = FALSE, quiet = FALSE, ...)

Arguments

Details

By default, the clade is identified by its number as found in the‘edge’ matrix of the tree. Iftips = TRUE, the tips descendingfrom the identified node are returned, possibly together with thenode. Iflabels = TRUE, the labels are returned (if the treehas no node labels, then the node numbered is returned).

The node is identified by the shortest distance where the clickoccurs. If the click occurs close to a tip, the function returns itsinformation.

Value

A list with one or two vectors named"tips" and/or"nodes" with the identification of the tips and/or of thenodes.

Note

This function does not add anything on the plot, but it can be wrappedwith, e.g.,nodelabels (see example), or its results canbe sent to, e.g.,drop.tip.

Author(s)

Emmanuel Paradis

See Also

plot.phylo,nodelabels,identify for the generic function

Examples

## Not run: tr <- rtree(20)f <- function(col) {    o <- identify(tr)    nodelabels(node=o$nodes, pch = 19, col = col)}plot(tr)f("red") # click close to a nodef("green")## End(Not run)

Plot of DNA Sequence Alignement

Description

This function plots an image of an alignment of nucleotide sequences.

Usage

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

Arguments

Details

The idea of this function is to allow flexible plotting and colouringof a nucleotide alignment. By default, the most common bases (a, g, c,t, and n) and alignment gap are plotted using a standard colourscheme.

It is possible to plot only one base specified aswhat with achosen colour: this might be useful to check, for instance, thedistribution of alignment gaps (image(x, "-")) or missing data(see examples).

Author(s)

Emmanuel Paradis, Klaus Schliep

See Also

DNAbin,del.gaps,alex,alview,all.equal.DNAbin,clustal,grid,image.AAbin

Examples

data(woodmouse)image(woodmouse)rug(seg.sites(woodmouse), -0.02, 3, 1)image(woodmouse, "n", "blue") # show missing dataimage(woodmouse, c("g", "c"), "green") # G+Cpar(mfcol = c(2, 2))### barcoding style:for (x in c("a", "g", "c", "t"))    image(woodmouse, x, "black", cex.lab = 0.5, cex.axis = 0.7)par(mfcol = c(1, 1))### zoom on a portion of the data:image(woodmouse[11:15, 1:50], c("a", "n"), c("blue", "grey"))grid(50, 5, col = "black")### see the guanines on a black background:image(woodmouse, "g", "yellow", "black")### Amino acidX <- trans(woodmouse, 2)image(X) # default ape colorsimage(X, scheme="Clustal") # Clustal coloring

Test for Binary Tree

Description

This function tests whether a phylogenetic tree is binary.

Usage

is.binary(phy)## S3 method for class 'phylo'is.binary(phy)## S3 method for class 'multiPhylo'is.binary(phy)## S3 method for class 'tree'is.binary(phy)

Arguments

Details

The test differs whether the tree is rooted or not. An urooted tree isconsidered binary if all its nodes are of degree three (i.e., threeedges connect to each node). A rooted tree is considered binary if allnodes (including the root node) have exactly two descendant nodes, sothat they are of degree three expect the root which is of degree 2.

The test ignores branch lengths. Consider usingdi2multiif you want to treat zero-branch lengths as resulting frommultichotomies.

is.binary.tree is deprecated and will be removed soon:currently it callsis.binary.

Value

a logical vector.

Author(s)

Emmanuel Paradis

See Also

is.rooted,is.ultrametric,multi2di

Examples

is.binary(rtree(10))is.binary(rtree(10, rooted = FALSE))is.binary(stree(10))x <- setNames(rmtree(10, 10), LETTERS[1:10])is.binary(x)

Check Compatibility of Splits

Description

is.compatible is a generic function with a method for the class"bitsplits". It checks whether a set of splits is compatibleusing thearecompatible function.

Usage

is.compatible(obj)## S3 method for class 'bitsplits'is.compatible(obj)arecompatible(x, y, n)

Arguments

Value

TRUE if the splits are compatible,FALSE otherwise.

Author(s)

Andrei Popescu

See Also

as.bitsplits

Is Group Monophyletic

Description

This function tests whether a list of tip labels is monophyletic on a given tree.

Usage

is.monophyletic(phy, tips, reroot = !is.rooted(phy), plot = FALSE, ...)

Arguments

Details

Ifphy is rooted, the test is done on the rooted tree, otherwisethe tree is first unrooted, then arbitrarily rerooted, in order to beindependent on the current position of the root. That is, the testasks iftips could be monophyletic given any favourably rootingofphy.

Ifphy is unrooted the test is done on an unrooted tree, unlessreroot = FALSE is specified.

If tip labels in the listtips are given as characters, they needto be spelled as in the objectphy.

Value

TRUE orFALSE.

Author(s)

Johan Nylanderjnylander@users.sourceforge.net

See Also

which.edge,drop.tip,mrca.

Examples

    ## Test one monophyletic and one paraphyletic group on the bird.orders tree    ## Not run: data("bird.orders")    ## Not run: is.monophyletic(phy = bird.orders, tips = c("Ciconiiformes", "Gruiformes"))    ## Not run: is.monophyletic(bird.orders, c("Passeriformes", "Ciconiiformes", "Gruiformes"))

Test if a Tree is Ultrametric

Description

This function tests whether a tree is ultrametric using the distancesfrom each tip to the root.

Usage

is.ultrametric(phy, ...)## S3 method for class 'phylo'is.ultrametric(phy, tol = .Machine$double.eps^0.5, option = 1, ...)## S3 method for class 'multiPhylo'is.ultrametric(phy, tol = .Machine$double.eps^0.5, option = 1, ...)

Arguments

Details

The test is based on the distances from each tip to the root and acriterion: ifoption = 1, the criterion is the scaled range((max - min/max)), ifoption = 2, the variance is used (thiswas the method used until ape 3.5). The default criterion is invariantto linear changes of the branch lengths.

Value

a logical vector.

Author(s)

Emmanuel Paradis

See Also

is.binary,.Machine

Examples

is.ultrametric(rtree(10))is.ultrametric(rcoal(10))

Plot Multiple Chronograms on the Same Scale

Description

The main argument is a list of (rooted) trees which are plotted on thesame scale.

Usage

kronoviz(x, layout = length(x), horiz = TRUE, ...,         direction = ifelse(horiz, "rightwards", "upwards"), side = 2)

Arguments

Details

The size of the individual plots is proportional to the size of thetrees.

Value

NULL

Author(s)

Emmanuel Paradis, Klaus Schliep

See Also

plot.phylo