phytools: an R package for phylogenetic comparative biology (and other things)

Comparative Biology

Several methods in phylogenetic comparative biology have been implemented in phytools. These cover a wide range of areas including ancestral character estimation (e.g. anc.trend), likelihood-based methods for studying the evolution of character traits over time (e.g. brownie.lite, evol.vcv, fitDiversityModel and phylosig), a Bayesian method for detecting the location of a rate shift in the tree (evol.rate.mcmc), estimation of phylogenetic signal, including with sampling error (phylosig), and various methods for statistical hypothesis testing in a phylogenetic context (e.g. phyl.cca, phyl.pairedttest, phyl.pca and phyl.resid).

Simulation

Several simulation methods are implemented in phytools. These include Brownian motion simulation under various conditions (fastBM), simulation of discrete character evolution (sim.history), simulation of stochastic character maps (make.simmap) and simulation of multiple evolutionary rates (sim.rates), among other functionality (Table 1).

Phylogenetic Inference

A few different phylogenetic inference procedures are implemented in the phytools package. These functions are, in general, highly dependent on calculations and algorithms in the phangorn library of Schliep (2011). Some of the functionality includes matrix representation parsimony supertree estimation (mrp.supertree) and least-squares phylogeny inference (optim.phylo.ls; Table 1).

Graphical Methods

Several graphical methods are implemented in phytools. Among these are projection of a tree into bivariate morphospace (phylomorphospace; Fig. 1), plotting stochastic character maps and histories (plotSimmap), lineage through time plotting with extinct lineages (ltt), animation of Brownian motion and speciation (branching.diffusion), and other functions (Table 1).

Utility Functions

In addition to the aforesaid scientific functions, phytools also includes a number of utility functions for phylogeny input, output and manipulation. These are meant to supplement and complement the existing diverse array of utility functions in the ape and phangorn packages. Several of these functions are listed in Table 1.

Example 1: Detecting the Location of a Rate Shift

In this example, I first simulate a stochastic pure-birth phylogeny; next, I simulate evolutionary change for a single continuously valued character on the phylogeny under two different evolutionary rates in different parts of the tree; I analyse the tree and data using the Bayesian MCMC method for identifying the location of a shift in the evolutionary rate over time (Revell et al. in press); finally, I analyse the MCMC results to estimate the location of the shift and the evolutionary rates tipward and rootward of this point.

First, I loaded the phytools package. This will also load ape and other required packages on first instantiation:

> # load the phytools package (and ape)
> require(phytools)

Loading required package: phytools

Loading required package: ape....

Next, I set the random number seed for reproducibility (here, it is just set to 1):

> set.seed(1) # set the seed

I use the ape function rbdtree to simulate a stochastic pure-birth tree. In this instance, the tree has 91 taxa.

> # simulate a tree (using ape)
> tree<-rbdtree(b=log(50),d = 0,Tmax=1)

Now, for the purposes of simulation, I split the tree at a predetermined position – specified here by the number of the descendant node and the distance along the edge from the root. To do this, I use the phytools function splitTree. It should be noted that the node and edge position used below are only guaranteed to work conditioned on having set the random number seed at 1 (see above), otherwise a different split point should be chosen.

> trees<-splitTree(tree,list(node = 153,bp = 0·09))

Now, I stretch the branches in one part of the tree and reattach the subtrees for simulation. To reattach the two parts of the tree, I use the phytools function paste.tree as follows:

> # stretch the branches in one part of the tree

> trees[[2]]$edge.length<-trees[[2]]$edge.length*10

> # reattach the two subtrees

> sim.tree<-paste.tree(trees[[1]],trees[[2]])

I can then plot the generating tree for simulation (which has its branches stretched to be proportional to the evolutionary rate multiplied by time; Fig. 2), using the phytools function plotTree, and simulate on this stretched tree using phytools function fastBM:

> # plot the generating tree for simulation

> plotTree(sim.tree,fsize = 0·5)

> x<-fastBM(sim.tree) # simulate on the tree

Now, I perform Bayesian MCMC analysis using the phytools function evol.rate.mcmc (Revell et al. in press). This analysis required about 20 min on a Dell i5 650 CPU running at 3·20 GHz.

> # perform MCMC > res<-evol.rate.mcmc(tree,x,ngen = 100000)

Control parameters (set by user or default):

List of 11

$ sig1 : num 4·34

$ sig2 : num 4·34

$ a : num 0·0222

$ sd1 : num 0·867

$ sd2 : num 0·867

$ sda : num 0·00444

$ kloc : num 0·2

$ sdlnr : num 1

$ rand.shift: num 0·05

$ print : num 100

$ sample : num 100

Starting MCMC run:

state sig1 sig2 a node bp likelihood

04·3352 4·3352 0·0222 35 0·345−130·7022

state sig1 sig2 a node bp likelihood

1004·8517 0·9447−0·0248 93 0·044−116·3554

state sig1 sig2 a node bp likelihood

2006·7915 1·2872−0·0205 153 0·057−101·0501

....

Done MCMC run.

The MCMC function first prints the control parameters (which can be set by the user, although above they have been given their default values, see below), and then prints the state of the MCMC chain at a frequency given by the control parameter print (here, every 100 generations; generations after 200 not shown above).

Next, I can estimate the location of the shift point by finding the split in the posterior sample with the smallest summed distance to all the other samples (this is one of multiple possible criteria; see Revell et al. in press). For this analysis, I use the phytools function minSplit and exclude the first 20 000 generations as burn-in:

>est.split<-minSplit(tree,res$mcmc[201:nrow(res$mcmc),])

This analysis took about 6 s to run on the same hardware as described earlier.

Finally, I need to pre-process the posterior sample to get the sampled rates tipward and rootward of the average shift, for each sample (see Revell et al. in press). I do this using the phytools function posterior.evolrate. I can then print the results (estimated shift point and evolutionary rates) to screen:

> pp.result<-posterior.evolrate(tree,est.split,res$mcmc[201:nrow(res$mcmc),], res$tips[201:nrow(res$mcmc)])

> est.sig1<-mean(pp.result[,“sig1”])

> est.sig2 <-mean(pp.result[,“sig2”])

> est.split

$node

[1] 153

$bp

[1] 0·1003489

> est.sig1

[1] 1·010709

> est.sig2

[1] 10·24914

This analysis took about 9 s to complete.

Here, the parameter estimates are very close to the generating shift point of [153, 0·09] and the generating evolutionary rates of and .

It should be noted that in actual practice, the authors should pay much closer attention to the control parameters of the MCMC than is given here, and in particular, to the proposal distribution for the model parameters. More information about function control can be obtained by calling the help file of evol.rate.mcmc:

> ?evol.rate.mcmc

or by referring to Revell et al. (in press). In addition, users should assess convergence and compute effective sample sizes for their samples from the posterior distribution. This can be carried out using the MCMC diagnostics package ‘coda’ (Plummer et al. 2006). Please refer to Revell et al. (in press) for more information about this method.

Example 2: Simulate and Analyse Multi-Rate Brownian Evolution

In this example, I first simulate the character history of a discretely valued character trait with three states evolving on a phylogeny. I then simulate the evolution of a continuous trait with a rate that depends on the value of the discrete trait. Finally, I fit single and multiple rate evolutionary models to the data and tree using the likelihood method of O’Meara et al. (2006).

After loading phytools, I first set the seed (arbitrarily to 10; done here for reproducibility only):

> set.seed(10)

Now, I simulate a stochastic pure-birth tree using ape:

> tree<-rbdtree(b = log(50),d = 0,Tmax = 1)

This tree contains 129 taxa. Next, I simulate a stochastic character history on the tree for a character with three states, A, B, and C, using the phytools function sim.history as follows:

> # this is our transition matrix
> Q<-matrix(c(−2,1,1,1,−2,1,1,1,−2),3,3)

> rownames(Q)<-colnames(Q)<-c(“A”,“B”,“C”)

> mtree<-sim.history(tree,Q)

I can plot the simulated history using the phytools function plotSimmap to see what it looks like:

> # set colors > cols<-c(“red”,“blue”,“green”);
> names(cols)<-rownames(Q)

> # plot tree with labels off
> plotSimmap(mtree,cols,ftype=“off”)

This visualization is shown in Fig. 3.

Next, I simulate continuous character evolution using three different rates using the phytools function sim.rates:

> # set rates
> sig2 <-c(1,10,100)
>names(sig2)<-rownames(Q)

> X<-sim.rates(mtree,sig2) # simulate

Finally, I fit a multi-rate Brownian model using the likelihood method of O’Meara et al. (2006) with the phytools function brownie.lite. This likelihood optimization took about 3 s to run on the same hardware described earlier.

> fit.bm<-brownie.lite(mtree,X,maxit=4000)

> fit.bm

$sig2.single

[1] 47·15693

$a.single

[1]−5·242448

$var.single

[1] 34·47709

$logL1

[1]−329·7582

$k1

[1] 2

$sig2.multiple

B C A

10·4646399 99·9563314 0·8247106

...

$logL.multiple

[1]−287·4271

$k2

[1] 4

$P.chisq

[1] 4·129091e-19

$convergence

[1] “Optimization has converged.”

It should be noted that the order of the three rate regimes in the fitted model is the order in which they are encountered in the tree (Fig. 3), rather than in alphabetical or numerical order. In this case, the fitted parameter estimates (0·82, 10·46, 99·96) are very close to their generating values (1, 10, and 100). For more details on this likelihood method, please refer to O’Meara et al. (2006) or Revell (2008).