Plotting species values at the tips of a tree using "pins"

Earlier this month I had the pleasure of hosting Sonal Singhal for a departmental seminar at UMass-Boston in our Biology Department seminar series.

One of the major projects that she talked about has just come out in the journal Science.

During Sonal’s talk, I noticed this cool phylogeny visualization in which tip data (in their case, tip rates of evolution and phenotypic innovation) are shown as “pins” – lines of different lengths, headed by different colored points.

Naturally, this post is about how to re-create this visualization style in base R graphics.

First, here’s an example using elopomorph eel total body length.

library(phytools)
data(eel.tree)
eel.tree<-ladderize(eel.tree,right=FALSE)
data(eel.data)
eel_size<-setNames(eel.data$Max_TL_cm,rownames(eel.data))
h<-max(nodeHeights(eel.tree))
plotTree(eel.tree,ftype="off",lwd=1,direction="upwards",ylim=c(0,2*h),
  mar=c(0.1,3.1,0.1,0.1))
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
polygon(c(0,max(pp$xx)+1,max(pp$xx)+1,0),
  c(1.1*h,1.1*h,1.9*h,1.9*h),border=FALSE,
  col="#F2F2F2")
lines(range(pp$xx),rep(1.1*h,2),col="black",lty="dotted")
rel.eel_size<-eel_size/max(eel_size)*0.75*h
segments(x0=pp$xx[1:Ntip(eel.tree)],y0=rep(1.1*h,Ntip(eel.tree)),
  x1=pp$xx[1:Ntip(eel.tree)],y1=rel.eel_size[eel.tree$tip.label]+1.1*h)
axis(2,at=1.1*h+max(rel.eel_size)/max(eel_size)*pretty(eel_size),
  labels=pretty(eel_size),las=1,cex.axis=0.6)
cols<-setNames(
  viridisLite::viridis(n=100)[ceiling(100*eel_size/max(eel_size))],
  names(eel_size))
points(pp$xx[1:Ntip(eel.tree)],rel.eel_size[eel.tree$tip.label]+1.1*h,
  pch=21,cex=1.5,bg=cols[eel.tree$tip])

Now let’s try multiple characters at once. These data are for Anolis lizard body size (SVL), and relative hindlimb length (HLL) and tail length (TL).

data(anoletree)
anole.tree<-as.phylo(anoletree)
data(anole.data)
anole_resid<-phyl.resid(anole.tree,x=as.matrix(anole.data[,"SVL",drop=FALSE]),
  Y=as.matrix(anole.data[,c(6,3)]))
anole_data<-cbind(anole_resid$resid,exp(anole.data[,"SVL",drop=FALSE]))
h<-max(nodeHeights(anole.tree))
plotTree(anole.tree,ftype="off",direction="upwards",ylim=c(0,4*h),
  mar=c(0.1,5.1,0.1,0.1),lwd=1)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
titles<-c("relative log(TL)","relative log(HLL)","SVL")
for(i in ncol(anole_data):1){
  d<-max(c(diff(range(anole_data[,i])),max(anole_data[,i])))
  x<-setNames(anole_data[,i]/d*0.8*h,
    rownames(anole.data))
  polygon(c(0,max(pp$xx)+1,max(pp$xx)+1,0),
    (i+c(0.05,0.05,0.95,0.95))*h,border=FALSE,
    col="#F2F2F2")
  hh<-(i+0.1)*h-min(c(0,min(x)))
  lines(range(pp$xx),rep(hh,2),col="black",lty="dotted")
  segments(x0=pp$xx[1:Ntip(anole.tree)],y0=rep(hh,Ntip(anole.tree)),
    x1=pp$xx[1:Ntip(anole.tree)],y1=x[anole.tree$tip.label]+hh)
  labs<-pretty(anole_data[,i])
  labs[!(labs>max(anole_data[,i]))]->labs
  labs[!(labs<min(anole_data[,i]))]->labs
  axis(2,at=hh+max(x)/max(anole_data[,i])*labs,
    labels=labs,las=1,cex.axis=0.6)
  cols<-setNames(
    viridisLite::viridis(n=100)[ceiling(99*((x-min(x))/diff(range(x))))+1],
    names(x))
  points(pp$xx[1:Ntip(anole.tree)],x[anole.tree$tip.label]+hh,
    pch=21,cex=1.2,bg=cols[anole.tree$tip.label])
  mtext(titles[i],2,line=3,at=mean(hh+max(x)/max(anole_data[,i])*labs),
    cex=0.8)
}

That’s it.

Of course, please check out this cool new paper in Science!

Very excited to share a huge project that's been a long time in the making. Interested in squamate macroevolution? Curious about how snakes fit into the broader ecological context of squamates? Looking for a new squamate phylogeny? Check out our new work: https://t.co/xomKFNbr0r pic.twitter.com/EXI99ckObi

— Pascal Title (@pascal_title) February 22, 2024