Co-phylogenetic plot with geological colors

A couple of days ago, I received the following phytools user email:

“I wonder whether it would be possible to use geo.legend along with cophylo. Alternatively, if you know of any workaround to add a geological time scale to a cophylo plot object. Any help would be very much appreciated.”

Well, in R the question is never if but always how (to paraphrase Simone Blomberg) – so let's see how.

For this example, I'll use a single estimated tree obtained from the Laurasiatherian dataset in the phangorn package.

I'm going to set the divergence date separating marsupial & placental mammals at 165 mybp based on http://www.timetree.org/, and then (crudely) time-calibrate the tree using chronos.

library(phytools)
tree<-read.nexus("LaurasiatherianML.nex")
tree<-drop.tip(root(tree,"Platypus"),"Platypus")
tree<-chronos(tree,calibration=makeChronosCalib(tree,age.min=165,
    age.max=165))

## 
## Setting initial dates...
## Fitting in progress... get a first set of estimates
##          (Penalised) log-lik = -17.00735 
## Optimising rates... dates... -17.00735 
## 
## log-Lik = -17.00713 
## PHIIC = 302.01

Here's what this tree looks like with the default geological timescale overlain.

plotTree(tree,ylim=c(-10,Ntip(tree)),fsize=0.6)
geo.legend()

Now what I'm going to do is subsample the tips of this tree – just so that I have some incongruence between my left & right trees in the cophylo plot.

In a more realistic use-case, perhaps we'd have a time-calibrated tree for host & parasite species, and then table of associations between the two.

left.tree<-drop.tip(tree,sample(tree$tip.label,6))
right.tree<-drop.tip(tree,sample(tree$tip.label,10))

Next, I'm going to create my plot. To do this, I'll first create & graph my "cophylo" object.

Then, I'll pull the "last_plot.cophylo" object that plot.cophylo creates and assign the left "cophylo" tree to a new object called "last_plot.phylo" that'll be used internally by geo.legend.

Then I'll make my geological colors palette using geo.palette. Actually, I could use a different set of time periods or colors if I wanted to – I'd just have to create a new table containing them.

Next, I'll rescale the time periods in my table to match those of my new plotting space.

This is critical because when we plot a "cophylo" object the space is no longer in the units of the edge lengths of our original tree!

Finally, I'll add the geological period legend using geo.legend – and then simply repeat the same thing for my other other plotted tree.

## plot cophylo
plot(cophylo(left.tree,right.tree),ylim=c(-0.15,1),fsize=0.6,
    link.type="curved",link.lty="solid",
    link.col=make.transparent("blue",0.25),link.lwd=2,
    pts=FALSE,part=0.47)

## Rotating nodes to optimize matching...
## Done.

## pull out the last plotted cophylo
obj<-get("last_plot.cophylo",envir=.PlotPhyloEnv)
## reset some attributes of the last plotted cophylo
tmp<-obj$left
tmp$direction<-FALSE
tmp$Ntip<-1
## assign it to last_plot.phylo
assign("last_plot.phylo",tmp,envir=.PlotPhyloEnv)
## extract geological color palette using geo.palette
pal<-geo.palette()
## find the time period that matches the depth of the 
## original tree
ii<-which(pal$period[,1]>max(nodeHeights(left.tree)))[1]
## subsample the geological palette matrix
pal$period[ii,1]<-max(nodeHeights(left.tree))
pal$colors<-pal$colors[1:ii]
## find the x-limits of the tree
XLIM<-range(tmp$xx)
## rescale to these new limits
new.pal<-XLIM[2]-pal$period/max(nodeHeights(left.tree))*diff(XLIM)
## add legend
geo.legend(leg=new.pal,colors=pal$colors,cex=0.7)
## repeat for the other tree
tmp<-obj$right
tmp$direction<-FALSE
tmp$Ntip<-1
assign("last_plot.phylo",tmp,envir=.PlotPhyloEnv)
pal<-geo.palette()
ii<-which(pal$period[,1]>max(nodeHeights(right.tree)))[1]
pal$period[ii,1]<-max(nodeHeights(right.tree))
pal$colors<-pal$colors[1:ii]
XLIM<-range(tmp$xx)
new.pal<-pal$period/max(nodeHeights(tree))*diff(XLIM)+XLIM[1]
geo.legend(leg=new.pal,colors=pal$colors,cex=0.7)

That's it.