Plotting a "contMap" tree with a custom, segmented color gradient

Recently, my student Kristin Winchell asked if it was possible to define a custom, segmented color gradient for a contMap plot. The idea is to plot estimated liabilities from an analysis conducted under the threshold model, and use a color gradient in which we have defined certain segments to correspond to the regions between each estimated threshold.

We can do this - but it's not completely straightforward. The strategy that I used employed was to generate custom gradients using colorRampPalette and then populate the corresponding sections of the color map gradient of our "contMap" object. For example:

library(phytools)
library(RColorBrewer)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  A, B, C, D, E, F, ...
## 
## Rooted; includes branch lengths.

x

##           A           B           C           D           E           F 
## -2.58294501 -0.87787706 -0.79231783 -1.45633035 -2.21125431 -1.18664645 
##           G           H           I           J           K           L 
## -0.73050079 -1.45626005 -1.64320288 -1.30336800 -1.77872106 -0.50475681 
##           M           N           O           P           Q           R 
##  1.31066282 -0.04524848  1.59883780  2.67919886  1.26184776  0.93059768 
##           S           T           U           V           W           X 
##  1.96968190  1.85420614  0.92747280 -0.81124266 -1.74906958 -1.72499681 
##           Y           Z 
##  0.19247578  0.17086918

obj<-contMap(tree,x,plot=FALSE)
## our thresholds
threshold<-c(0,1)
## our palette for the segments
colors<-brewer.pal(n=7,name="YlGnBu")
## our new colors
lims<-obj$lims
colvals<-seq(lims[1],lims[2],by=diff(lims)/(length(obj$cols)-1))
ii<-which((colvals-threshold[1])^2==min((colvals-threshold[1])^2))
newcols<-colorRampPalette(colors[1:3])(ii)
jj<-which((colvals-threshold[2])^2==min((colvals-threshold[2])^2))
newcols<-c(newcols,colorRampPalette(colors[3:5])(jj-ii))
newcols<-c(newcols,colorRampPalette(colors[5:7])(length(colvals)-jj))
obj$cols[]<-newcols
## desired lengend width
leg.width<-8
## our plot
h<-max(nodeHeights(obj$tree))
plot(obj,legend=FALSE,xlim=c(-0.25*h,1.1*h),lwd=6)
add.color.bar(Ntip(obj$tree)-3,obj$cols,title="liability",
    lims=NULL,digits=3,direction="upwards",
    subtitle="",lwd=leg.width,x=-0.2*h,y=2,prompt=FALSE)
LWD<-diff(par()$usr[1:2])/dev.size("px")[1]
lines(x=rep(-0.2*h+LWD*leg.width/2,2),y=c(2,Ntip(obj$tree)-1))
ticks<-c(lims[1],threshold,lims[2])
nticks<-length(ticks)
Y<-cbind((ticks-min(ticks))/diff(range(ticks))*(Ntip(obj$tree)-3)+2,
    (ticks-min(ticks))/diff(range(ticks))*(Ntip(obj$tree)-3)+2)
X<-cbind(rep(-0.2*h+LWD*leg.width/2,nticks),
    rep(-0.2*h+LWD*leg.width/2+0.02*h,nticks))
for(i in 1:nrow(Y)) lines(X[i,],Y[i,])
text(x=X[,2],y=Y[,2],round(ticks,3),pos=4,cex=0.8)

We can also select completely divergent color gradients for each segment. For example:

leg.width<-8
## our palette for the segments
colors<-list(c("yellow","green"),
    c("green","blue","purple"),
    c("purple","red"))
## our new colors
lims<-obj$lims
colvals<-seq(lims[1],lims[2],by=diff(lims)/(length(obj$cols)-1))
ii<-which((colvals-threshold[1])^2==min((colvals-threshold[1])^2))
newcols<-colorRampPalette(colors[[1]])(ii)
jj<-which((colvals-threshold[2])^2==min((colvals-threshold[2])^2))
newcols<-c(newcols,colorRampPalette(colors[[2]])(jj-ii))
newcols<-c(newcols,colorRampPalette(colors[[3]])(length(colvals)-jj))
obj$cols[]<-newcols
## our plot
plot(obj,legend=FALSE,xlim=c(-0.25*h,1.1*h),lwd=6)
add.color.bar(Ntip(obj$tree)-3,obj$cols,title="liability",
    lims=NULL,digits=3,direction="upwards",
    subtitle="",lwd=leg.width,x=-0.2*h,y=2,prompt=FALSE)
LWD<-diff(par()$usr[1:2])/dev.size("px")[1]
lines(x=rep(-0.2*h+LWD*leg.width/2,2),y=c(2,Ntip(obj$tree)-1))
ticks<-c(lims[1],threshold,lims[2])
nticks<-length(ticks)
Y<-cbind((ticks-min(ticks))/diff(range(ticks))*(Ntip(obj$tree)-3)+2,
    (ticks-min(ticks))/diff(range(ticks))*(Ntip(obj$tree)-3)+2)
X<-cbind(rep(-0.2*h+LWD*leg.width/2,nticks),
    rep(-0.2*h+LWD*leg.width/2+0.02*h,nticks))
for(i in 1:nrow(Y)) lines(X[i,],Y[i,])
text(x=X[,2],y=Y[,2],round(ticks,3),pos=4,cex=0.8)

It is no problem to plot this in a fan style if desired. We just have to move the legend:

h<-max(nodeHeights(obj$tree))
plot(obj,type="fan",legend=FALSE,xlim=c(-1.5*h,1.1*h))
add.color.bar(2*h,obj$cols,title="liability",
    lims=NULL,digits=3,direction="upwards",
    subtitle="",lwd=leg.width,x=-1.4*h,y=-h,prompt=FALSE)
LWD<-diff(par()$usr[1:2])/dev.size("px")[1]
lines(x=rep(-1.4*h+LWD*leg.width/2,2),y=c(-h,h))
ticks<-c(lims[1],threshold,lims[2])
nticks<-length(ticks)
Y<-cbind((ticks-min(ticks))/diff(range(ticks))*2*h-h,
    (ticks-min(ticks))/diff(range(ticks))*2*h-h)
X<-cbind(rep(-1.4*h+LWD*leg.width/2,nticks),
    rep(-1.4*h+LWD*leg.width/2+0.04*h,nticks))
for(i in 1:nrow(Y)) lines(X[i,],Y[i,])
text(x=X[,2],y=Y[,2],round(ticks,3),pos=4,cex=0.8)

Something like that.

I'm going to simulate a larger dataset to make it look cooler:

bigger.tree<-pbtree(n=150,scale=1)
y<-fastBM(bigger.tree,a=1)
obj<-contMap(bigger.tree,y,plot=FALSE)
leg.width<-8
colors<-list(c("yellow","green"),
    c("green","blue","purple"),
    c("purple","red"))
lims<-obj$lims
colvals<-seq(lims[1],lims[2],by=diff(lims)/(length(obj$cols)-1))
ii<-which((colvals-threshold[1])^2==min((colvals-threshold[1])^2))
newcols<-colorRampPalette(colors[[1]])(ii)
jj<-which((colvals-threshold[2])^2==min((colvals-threshold[2])^2))
newcols<-c(newcols,colorRampPalette(colors[[2]])(jj-ii))
newcols<-c(newcols,colorRampPalette(colors[[3]])(length(colvals)-jj))
obj$cols[]<-newcols
h<-max(nodeHeights(obj$tree))
plot(obj,type="fan",legend=FALSE,xlim=c(-1.3*h,h),ftype="off")
add.color.bar(2*h,obj$cols,title="liability",
    lims=NULL,digits=3,direction="upwards",
    subtitle="",lwd=leg.width,x=-1.2*h,y=-h,prompt=FALSE)
LWD<-diff(par()$usr[1:2])/dev.size("px")[1]
lines(x=rep(-1.2*h+LWD*leg.width/2,2),y=c(-h,h))
ticks<-c(lims[1],threshold,lims[2])
nticks<-length(ticks)
Y<-cbind((ticks-min(ticks))/diff(range(ticks))*2*h-h,
    (ticks-min(ticks))/diff(range(ticks))*2*h-h)
X<-cbind(rep(-1.2*h+LWD*leg.width/2,nticks),
    rep(-1.2*h+LWD*leg.width/2+0.04*h,nticks))
for(i in 1:nrow(Y)) lines(X[i,],Y[i,])
text(x=X[,2],y=Y[,2],round(ticks,3),pos=4,cex=0.8)

That's OK.