Combining phytools phylo.to.map and contMap plots

A friend & colleague recently asked the following:

I was wondering whether it was possible to combine the output of contMap with phylo.to.map so, that the evolution tree is shown in the gradient used from the continuous trait evolution plot?

This would be a cool feature, for sure. Unfortunately the plotter that is used internally by the method plot.phylo.to.map is not the same as is used by plot.contMap, so this could get complicated.

Fortunately, in R (to paraphrase someone wiser than I) the question is always how not if (or, in this case, whether). So, let's try to figure out the “how.”

Firstly, I had to push a few small updates to phylo.to.map.R.

These weren't strictly necessary, but they'll sure make things easier. As such, to follow along exactly, the reader should update phytools from GitHub using devtools. (With devtools installed this can be done just by running devtools::install_github("liamrevell/phytools") in an interactive R session.)

library(phytools)
packageVersion("phytools")

## [1] '1.0.4'

For the purposes of this example, I'll imagine that I have my phylogeny (Argentina.tree), my lat/long data (Argentina.latlong), and my trait data (Argentina.trait) already in R.

In practice, we'd normally start off by reading these data from an input file or files of some kind.

ls()

## [1] "Argentina.latlong" "Argentina.trait"   "Argentina.tree"    "cMap"             
## [5] "lastPP"            "obj"               "tt"

Argentina.tree

## 
## Phylogenetic tree with 12 tips and 11 internal nodes.
## 
## Tip labels:
##   P.bnqvz, I.vqcxyl, P.nrufw, C.wpndrel, P.xusz, I.gtjyvdriec, ...
## 
## Rooted; includes branch lengths.

Argentina.latlong

##                    lat      long
## P.bnqvz      -54.34253 -67.86878
## I.vqcxyl     -50.58238 -70.52681
## P.nrufw      -47.14639 -70.91579
## C.wpndrel    -47.66503 -67.80395
## P.xusz       -40.92270 -70.00817
## I.gtjyvdriec -43.12693 -67.35014
## G.pzawf      -39.04263 -63.46033
## U.mkflj      -36.19011 -68.90606
## V.ezsji      -35.80113 -61.70992
## C.pyrnvj     -30.93887 -61.77475
## E.suqlrtvh   -31.06853 -65.40523
## G.itjxbykps  -27.50288 -66.18320

Argentina.trait

##      P.bnqvz     I.vqcxyl      P.nrufw    C.wpndrel       P.xusz I.gtjyvdriec 
##   3.43744323   0.43327903  -1.25185738  -7.42113286 -10.57345007  -7.91838325 
##      G.pzawf      U.mkflj      V.ezsji     C.pyrnvj   E.suqlrtvh  G.itjxbykps 
## -10.67775857 -10.65758751  -6.74674621  -0.07399385  -0.53337834  -1.91787840

You can already see that these data are totally made up (simulation code here) – based on the gibberish taxon labels – but I created them to try to emulate reasonable looking tip labels, just to be sure we're properly testing our plotting procedure.

The next step is to simply create our "phylo.to.map" object in the usual way, which we'll do here as follows:

obj<-phylo.to.map(Argentina.tree,Argentina.latlong,
    regions="Argentina",direction="rightwards",plot=FALSE,)

## objective: 6
## objective: 6
## objective: 6
## objective: 6
## objective: 6
## objective: 6
## objective: 6
## objective: 6
## objective: 6
## objective: 6

## objective: 4

obj

## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 12 tips and 11 internal nodes.
## 
## (2) A geographic map with range:
##      -55.03N, -21.8N
##      -73.58W, -53.67W.
## 
## (3) A table containing 12 geographic coordinates (may include
##     more than one set per species).
## 
## If optimized, tree nodes have been rotated to maximize alignment
## with the map when the tree is plotted in a rightwards direction.

Next, I'll create a "contMap" object, but in this case it's very important that I do that using the (potentially rotated) tree from "phylo.to.map" (here obj$tree) rather than my original input phylogeny.

cMap<-contMap(obj$tree,Argentina.trait,plot=FALSE)
cMap

## Object of class "contMap" containing:
## 
## (1) A phylogenetic tree with 12 tips and 11 internal nodes.
## 
## (2) A mapped continuous trait on the range (-10.677759, 3.437443).

Now I'll recolor this object, just because the default color scheme for "contMap" is pretty bad. (Forgive me. I'm colorblind.)

cMap<-setMap(cMap,c("black","red","yellow","white"))

Finally, in one long code chunk, I'll plot my "phylo.to.map" object, then I'll pull out the coordinates (etc.) of the graphed tree, then I'll rescale my "contMap" object to match the total tree length of my plotted phylogeny, then I'll overplot my "contMap" object (first graphing a tree in black to create an outline effect), lining it up the the tree I'd already plotted, finally I'll add a color gradient legend.

Phew!

## plot our phylogeny & map
plot.phylo.to.map(obj,direction="rightwards",ftype="i",
    colors=c("blue","white"),pts=FALSE,psize=1.2,
    map.bg="slategray",fsize=0.8)
## get "last_plot.phylo" object
lastPP<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## rescale our "contMap" tree
tt<-cMap$tree
tt<-rescaleSimmap(tt,depth=diff(range(lastPP$xx)))
## plot a tree using the default fg color for an outline
## effect
plotTree(as.phylo(tt),add=TRUE,lwd=8,ftype="off",
    xlim=lastPP$x.lim-lastPP$x.lim[1],
    ylim=lastPP$y.lim,
    tips=lastPP$yy[1:Ntip(cMap$tree)],
    mar=rep(0.01,4),asp=1,color=par()$fg)
## graph the contMap tree
plot(tt,cMap$cols,add=TRUE,lwd=6,ftype="off",
    xlim=lastPP$x.lim-lastPP$x.lim[1],
    ylim=lastPP$y.lim,
    tips=lastPP$yy[1:Ntip(cMap$tree)],
    mar=rep(0.01,4),asp=1)
## add a color gradient legend
add.color.bar(10,cMap$cols,"phenotype",lims=cMap$lims,
    digits=2,prompt=FALSE,lwd=6,outline=TRUE,x=22,y=-52,
    subtitle="")

Incredible! It works!