R Hackathon 1/Ancestral State Reconstruction
The primary ancestral state reconstruction algorithms in R are accessed through the function "ace" in package "ape". To work through the following example using the Geospiza dataset, make sure that you have installed and loaded ape into your R session and loaded the Geospiza phylogeny and tip data into memory.
library(ape) geotree <- read.nexus("geospiza.nex") geodata <- read.table("geospiza.txt")
Tip data is not available for the outgroup "olivacea", so drop that taxon from the analysis.
geotree <- drop.tip(geotree, "olivacea")
IMPORTANT. The row.names of your dataframe must match the tip.labels of your phylogeny. In the example above, the row.names for geodata will match the tip.labels for geotree after "olivacea" has been culled. However, the individual columns of geodata (e.g. geodata$wingL) do not automatically have the row.names of the whole data table associated with them! If you call a column of the data.table with ace without first dealing with this issue, the analysis will run, but the tip data will be disassociated from the proper tips! There are two workarounds.
One option is to sort the datatable so that the taxa appear in the same order in the table as in the phylogeny. For example:
geodata <- geodata[geotree$tip.label, ]
The other, and better option is to extract the column of data of interest from the overall datatable, creating a new vector, and to transfer the row.names of the datatable to the NAMES of the vector. This is the option least likely to inadverently disassociate the tip data from the tips.
wingL <- geodata$wingL names(wingL) <- row.names(geodata)
The worked examples below assume that you have used the second solution (vector extraction and name assignment).
Reconstructing Ancestral States for Continuous Variables
There are three general options for continuous variables: reconstructions based on maximum likelihood (ML), reconstructions based upon phylogenetic independent contrasts (pic) and reconstructions based on generalized least squares (GLS). Be aware that the generalized least squares algorithms seem to give spurious results near the root of the phylogeny at present.
MLreconstruction <- ace(wingL, geotree, type="continuous", method="ML")
Why am I getting all these warning messages?
ML reconstruction using ace tends to generate a large number of not-a-number error messages. These result when the program calculates the likelihood of particularly poor fits. You can safely ignore these messages.
PHYLOGENETIC INDEPENDENT CONSTRASTS
picreconstruction <- ace(wingL, geotree, type="continuous", method="pic")
GENERALIZED LEAST SQUARES
Reconstructions based upon generalized least squares require specification of a correlation structure for the generalized linear model. There are three basic options for the specification of the correlation structure: 1) corBrownian, which uses a simple Brownian motion model, corMartins, which uses an Ornstein-Uhlenbeck (constrained random-walk) model, and corGrafen, which is a modified Brownian motion model (Grafen 1989) (ADD MORE INFO ON WHAT EXACTLY THE GRAFEN MODEL DOES).
This is the syntax for the simple Brownian model
GLSreconstruction <- ace(wingL, geotree, type="continuous", method="GLS", corStruct = corBrownian(1, geotree))
The syntax for the Ornstein-Uhlenbeck model requires specifying an alpha parameter (here, 0.5).
GLSreconstruction <- ace(wingL, geotree, type="continuous", method="GLS", corStruct = corMartins(0.5, geotree))
The syntax for the Grafen model requires specifying a rho parameter (here, 1) which exponentiates the recalculated branch lengths.
GLSreconstruction <- ace(wingL, geotree, type="continuous", method="GLS", corStruct = corGrafen(1, geotree))
WHAT ABOUT SQUARED CHANGE PARSIMONY?
geotreeones <- compute.brlen(geotree, 1) SQPreconstruction <- ace(wingL, geotreeones, type="continuous", method="ML")
HOW DO I USE THE OUTPUT FROM ACE?
returns the vector of node values for wing length (wingL) that were reconstructed using maximum likelihood, and
wingLfinal <- c(wingL, MLreconstruction$ace)
concatenates the original tip data with the reconstructed node data into a single vector.
HOW DO I PLOT THE OUTPUT FROM ACE? (drawn from a course handout by Gene Hunt)
plot.phylo(geotree) tiplabels(pch = 21, cex=wingL) nodelabels(pch = 21, cex=GLSreconstruction1$ace) ## (pch = 21 is just telling plot.phylo which symbol to use)
HOW DO I FIT A MODEL OF CHARACTER CHANGE TO CONTINUOUS DATA?
While ace returns a rate of Brownian evolution (e.g. MLreconstruction$sigms2), that rate and its associated likelihood vale (e.g. MLreconstruction$loglik) is conditioned on the specific ancestral states reconstructed by ace. More generally applicable model fits using Brownian motion and and the Ornstein-Uhlenbeck model are available in the packages ouch and geiger. Please see here for more information.
Reconstructing Ancestral States for Discrete Variables
char1<-c(1,1,1,1,1,0,1,0,0,1,0,1,1,1) names(char1)<-geotree$tip.label char1