# Difference between revisions of "R Hackathon 1/Ancestral State Reconstruction"

(→Reconstructing Ancestral States for Discrete Variables) |
(→Reconstructing Ancestral States for Discrete Variables) |
||

Line 131: | Line 131: | ||

Test ER against ARD for Geospiza data. | Test ER against ARD for Geospiza data. | ||

+ | |||

+ | P value | ||

+ | |||

+ | 1-pchisq(2*(abs(ERreconstruction$loglik) - abs(ARDreconstruction$loglik)), 1) | ||

+ | |||

+ | Just the absolute difference in log likelihoods. | ||

+ | |||

+ | (abs(ERreconstruction$loglik) - abs(ARDreconstruction$loglik) | ||

How do I define a custom transition model. | How do I define a custom transition model. |

## Revision as of 13:26, 13 December 2007

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.

**MAXIMUM LIKELIHOOD**

This syntax will reconstruct the ancestral states for the variable "wingL" (extracted from geodata) using the Brownian motion-based maximum likelihood (ML) estimator of Schluter et. al. (1997). This is the default method.

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**

This syntax will reconstruct the ancestral states for the variable "wingL" (extracted from geodata) using Felsenstein's (1985) phylogenetic independent contrasts (pic). This is also a Brownian-motion based estimator, but it only takes descendants of each node into account in reconstructing the state at that node. More basal nodes are ignored.

picreconstruction <- ace(wingL, geotree, type="continuous", method="pic")

**GENERALIZED LEAST SQUARES**

The generalized least squares algorithm in ace seems to given spurious results near the root of the phylogeny, at least at present. Specifically, the root node is always reconstructed at possessing state zero, and internal nodes near the root may receive reconstructed values out of the range of values observed among the tips.

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?**

Squared change parsimony is mathematically equivalent to a special case of Schluter et. al.'s maximum likelihood method in which branch length information is ignored. Thus, to perform a squared change parsimony reconstruction, set all branch lengths equal to 1 and calculate a maximum likelihood solution.

geotreeones <- compute.brlen(geotree, 1) SQPreconstruction <- ace(wingL, geotreeones, type="continuous", method="ML")

**HOW DO I USE THE OUTPUT FROM ACE?**

Each of the objects (for example: MLreconstruction) that we reconstructed above is a list of several elements. These are loglik (the log-likelihood of the most likely reconstruction), ace (the vector of reconstructed node values), sigma2 (a two element vector of the rate estimate and its standard error) and CI95, the 95% confidence intervals around the vector of reconstructed node values. Thus

MLreconstruction$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 can scale symbols at tips and nodes of your tree to your character data. Try this to visualize the evolution of wing length (wingL) across the Geospiza phylogeny.

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)

Note that the blue internal nodes are all similar in size to the yellow tip labels. This is because there is little variation in wing length across the dataset.

**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

Set up space.

Reload tree.

Set up new characters.

char1<-c(1,1,1,1,1,0,1,0,0,1,0,1,1,1) names(char1)<-geotree$tip.label

Three different options, or set your own matrix.

Run these using the geospiza data.

ERreconstruction <- ace(char1, geotree, type="discrete", model="ARD") SYMreconstruction <- ace(char1, geotree, type="discrete", model="SYM") ARDreconstruction <- ace(char1, geotree, type="discrete", model="ARD")

For Two states, ER = SYM. Explain why, and how three states differs,

Test ER against ARD for Geospiza data.

P value

1-pchisq(2*(abs(ERreconstruction$loglik) - abs(ARDreconstruction$loglik)), 1)

Just the absolute difference in log likelihoods.

(abs(ERreconstruction$loglik) - abs(ARDreconstruction$loglik)

How do I define a custom transition model.

Defining a custom four-parameter matrix for a three state character.

transitions <- matrix (c(0, 1, 2, 1, 0, 3, 2, 4, 0), nrow=3)

Matrix should look like.

[,1] [,2] [,3] [1,] 0 1 2 [2,] 1 0 4 [3,] 2 3 0