Plotting facing "contMap" objects with the same trait range & color palette

A phytools user recently asked the following question:

“Dear Liam. I am using the command contMap to do an ancestral reconstruction of a continuous trait. The command is working well. My question is that I want to compare two different trees, one for the data for females and another one for the data of males (facing one another). My problem is that the range of the trait for females is much narrower than the range of the trait for males. So, for example, the highest value for females is 0.6 (which is the brightest red color in the tree) while for males the highest value is 0.9 (being the brightest red) and the 0.6 value in males turns out to be yellowish. I would like to manually set the trait value from 0 to 1 in both trees so that the 0.6 value is the same color in both.”

This is actually surprisingly easy to accomplish using phytools::contMap. Here I’ll show how using the 1st & 3rd principal components scores from a morphological trait analysis of cordylid lizards by Broeckhoven et al. (2016). (I picked these because the 1st & 2nd had very similar numerical ranges!)

First, let’s see what happens by default.

## load packages
library(phytools)

## load tree and data
data(cordylid.tree)
data(cordylid.data)

## pull out PCs 1 & 2
pc1<-setNames(cordylid.data$pPC1,rownames(cordylid.data))
pc3<-setNames(cordylid.data$pPC3,rownames(cordylid.data))

## compute contMaps & update color palette
pc1_cmap<-contMap(cordylid.tree,pc1,plot=FALSE)
pc1_cmap<-setMap(pc1_cmap,viridisLite::viridis(n=20))
pc3_cmap<-contMap(cordylid.tree,pc3,plot=FALSE)
pc3_cmap<-setMap(pc3_cmap,viridisLite::viridis(n=20))

## graph our two contMaps
par(mfrow=c(1,2))
plot(pc1_cmap,leg.txt="PC 1")
plot(pc3_cmap,leg.txt="PC 3")

We can see that even though our two different principle components have substantively different numerical ranges, the range on the color palette is (by design) identical.

Sometimes, as per the user’s comment, we might want to force a particular numerical range on the "contMap" object. This is (surprisingly) quite easily done in phytools::contMap with the argument lims.

In the code below, I’ll do this & also illustrate how to create a facing graph with the tip labels and a separately plotted legend below. This is actually much harder than fixing the trait limits to be equivalent between the two plots.

## fix range
x_range<-range(c(pc1,pc3))
## compute contMaps & update color palette
pc1_cmap<-contMap(cordylid.tree,pc1,plot=FALSE,
  lims=x_range)
pc1_cmap<-setMap(pc1_cmap,viridisLite::viridis(n=20))
pc3_cmap<-contMap(cordylid.tree,pc3,plot=FALSE,
  lims=x_range)
pc3_cmap<-setMap(pc3_cmap,viridisLite::viridis(n=20))

## use layout instead of par(mfrow)
mat<-matrix(c(1,3,2,4,4,4),byrow=TRUE,2,3)
## this may need to be adjusted based on size of
## plotting device
layout(mat,widths=c(0.43,0.14,0.43),heights=c(0.9,0.1))
## graph contMaps
plot(pc1_cmap,ftype="off",legend=FALSE,mar=c(0.3,0.3,2.1,0.3))
mtext("a) PC 1",adj=0.05,line=-0.5,cex=0.9)
plot(pc3_cmap,ftype="off",legend=FALSE,
  direction="leftwards",mar=c(0.3,0.3,2.1,0.3))
mtext("b) PC 3",adj=0.05,line=-0.5,cex=0.9)
## identify aspect ratio in user/inches
xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
xasp<-diff(xlim)/par()$pin[1]
## get ylims for tip labels
ylim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$y.lim
## add tip labels back in
plot(NA,axes=FALSE,bty="n",xlim=c(-1,1),ylim=ylim)
text(x=0,y=1:Ntip(pc1_cmap$tree),
  gsub("_"," ",pc1_cmap$tree$tip.label),
  font=3)
## open plot for legend
par(mar=rep(0,4))
plot(NA,axes=FALSE,bty="n",xlim=c(-1,1),ylim=c(-1,1))
## compute aspect ratio of new plot
new.xasp<-2/par()$pin[1]
## add color bar
add.color.bar(leg=new.xasp/xasp,cols=pc1_cmap$cols,
  title="PC score",lims=pc1_cmap$lims,digits=2,
  prompt=FALSE,lwd=6,x=-0.5*(new.xasp/xasp),y=0,
  subtitle="length = 1.0 time units")

That’s it folks!