## Friday, October 27, 2017

### 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)
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)
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))
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")
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.