Placing missing taxa and fossils on a phylogenetic tree using continuous traits

In this exercise, we will use a method locate.yeti in the phytools package to attempt to place a missing species, the recently extinct taxonn Anolis roosevelti, into the tree of all anoles using quantitative characters.

First we need to load our packages:

library(phytools)

Next, let's read in our tree & data from file. This is (1) the Anolis tree that we've seen already - but it also includes a larger phenotypic trait dataset for about 20 different characters.

tree<-read.tree("Revell-etal.tree.tre")
## or
tree<-read.tree("http://www.phytools.org/Ilhabela2015/data/Revell-etal.tree.tre")
tree

## 
## Phylogenetic tree with 100 tips and 99 internal nodes.
## 
## Tip labels:
##  ahli, allogus, rubribarbus, imias, sagrei, bremeri, ...
## 
## Rooted; includes branch lengths.

X<-read.csv("Revell-etal.data.csv",header=T,row.names=1)
## or
X<-read.csv("http://www.phytools.org/Ilhabela2015/data/Revell-etal.data.csv",
    header=T,row.names=1)
X<-as.matrix(X)
X<-X[,1:20]
tip<-setdiff(rownames(X),tree$tip.label)
tip

## [1] "roosevelti"

OK, now we are ready to run our estimation procedure in which we will try to place the missing lineage, tip, into our base tree.

mltree<-locate.yeti(tree,X,search="exhaustive",plot=TRUE)

## Optimizing the phylogenetic position of roosevelti using ML. Please wait....

## Done.

mltree$logL

## [1] 3392.184

We have estimated our tree. Let's visualize it and overlay arrows (using the phytools function, add.arrow) pointing to both Anolis roosevelti and the species that we hypothesized a priori would be A. roosevelti's sister taxon.

mltree<-reorder(reorder(mltree,"pruningwise")) ## reorder hack
pca<-phyl.pca(mltree,X)
## flip PC 1 so that it loads positively (not negatively) on size
obj<-contMap(mltree,-pca$S[,1],plot=FALSE)
class(obj)<-"densityMap"
## fix tip label
obj$tree$tip.label[which(obj$tree$tip.label=="pumilis")]<-"pumilus"
## plot
plot(obj,type="fan",lwd=4,legend=0.7,
    leg.txt=c(round(obj$lims[1],3),"PC 1",round(obj$lims[2],3)),
    outline=TRUE)
## add arrows
add.arrow(tree=obj$tree,tip="roosevelti",col="red",lwd=3,hedl=0.06,angle=50)
add.arrow(tree=obj$tree,tip="cuvieri",col="blue",lwd=3,hedl=0.06,angle=50)

Finally, we can do the analyses to test the (null) hypotheses that A. roosevelti is sister to A. cuvieri, or sister-to or nested-within the majority of the other Puerto Rican anoles:

## find ML cuvieri-constraint tree
cuvieri<-which(tree$tip.label=="cuvieri")
cuvieri.tree<-locate.yeti(tree,X,search="exhaustive",constraint=cuvieri)

## Optimizing the phylogenetic position of roosevelti using ML. Please wait....
## Done.

## compute a LR to test the alternative hypothesis that the A. roosevelti
## is *not* sister
LR.cuvieri<-2*(mltree$logL-cuvieri.tree$logL)
## perform simulation for the null distribution of the LR
set.seed(1)
LR.null<-vector()
nsim<-20
obj<-phyl.vcv(X[cuvieri.tree$tip.label,],vcv(cuvieri.tree),lambda=1)
for(i in 1:nsim){
    Xsim<-sim.corrs(tree=cuvieri.tree,vcv=obj$R,anc=obj$alpha[,1])
    mltree.null<-locate.yeti(tree,Xsim,search="exhaustive",plot=FALSE,
        quiet=TRUE)
    cuvieri.tree.null<-locate.yeti(tree,Xsim,search="exhaustive",
        constraint=which(tree$tip.label=="cuvieri"),plot=FALSE,
        quiet=TRUE)
    LR.null[i]<-2*(mltree.null$logL-cuvieri.tree.null$logL)
}
P.cuvieri<-1-mean(LR.cuvieri>=LR.null)
## P-value of the LR test
P.cuvieri

## [1] 0.3

This suggests that we cannot reject the hull hypothesis that A. roosevelti is sister to A. cuvieri.

Challenge question

Test the null hypothesis that Anolis roosevelti is nested-within or sister- to the clade of Anolis lizards consisting of:

sp<-c("cristatellus","cooki","poncensis","gundlachi","pulchellus","krugi",
    "stratulus","evermanni")

Can we reject this placement for Anolis roosevelti?

Written by Liam J. Revell. Last updated July 5, 2015.