Graphing a discrete trait and a text label for each tip of a plotted tree

Today a phytools user emailed me with the following question.

“I wanted to ask you for help with something (I am desperate). I am trying to create a phylogenetic tree where the tip labels have multiple columns, which have labels themselves. The first column is the species name; the second is a dot indicating the value of a discrete character trait; the third is the IUCN status. Here is a drawing showing what I mean:”

OK – this may look like a scribble, but it’s exactly what I like to see when someone asks me for help with graphing. Kudos!

Actually, this isn’t hard to do at all. By the looks of their plot, in addition to the species names, it looks like they’d like to graph \(\lambda\) (I assume the rate of population growth, with \(\lambda = 1\)), but perhaps with all \(\lambda > 1\) shown in one color & \(\lambda < 1\) in another, and then IUCN threat status (i.e., EX, EW, CR, EN, VU, NT, and LC).

Before I continue, let me add that this is exactly the type of custom plot that Luke Harmon & I focus on in Chapter 13 of our recent book with Princeton University Press.

To illustrate how one might go about making this kind of plot, I’ll use a real tree – this phylogeny of 26 species of Plethodon salamanders from Highton & Larson (1979).

library(phytools)
data(salamanders)
salamanders$tip.label

##  [1] "neomexicanus" "larselli"     "vandykei"     "idahoensis"   "vehiculum"   
##  [6] "dunni"        "elongatus"    "stormi"       "fourchensis"  "ouachitae"   
## [11] "caddoensis"   "yonahlossee"  "jordani"      "glutinosus"   "punctatus"   
## [16] "wehrlei"      "websteri"     "dorsalis"     "welleri"      "serratus"    
## [21] "hubrichti"    "cinereus"     "nettingi"     "hoffmani"     "shenandoah"  
## [26] "richmondi"

To continue, let’s populate some IUCN values. (Don’t take my word on it, but I believe these to be the true IUCN classification of these 26 species!)

iucn<-setNames(
  c("NT","NT","LC","LC","LC",
    "LC","NT","EN","NT","NT",
    "NT","LC","NT","LC","NT",
    "LC","LC","LC","EN","LC",
    "VU","LC","NT","LC","VU",
    "LC"),salamanders$tip.label)

Now, I’m just going to generate some random values of \(\lambda\) between \(\lambda = 0.95\) (declining) and \(\lambda = 1.05\) (increasing), without presuming that an endangered species is in decline or a species of least concern is growing. (Indeed, this might not even be expected to be the case – for example, if conservation or recovery efforts were focused on vulnerable or endangered species.)

lambda<-setNames(
  runif(n=Ntip(salamanders),0.95,1.05),
  salamanders$tip.label)
lambda

## neomexicanus     larselli     vandykei   idahoensis    vehiculum        dunni    elongatus 
##    0.9565825    0.9887832    1.0361898    1.0020955    0.9516036    0.9502447    1.0105711 
##       stormi  fourchensis    ouachitae   caddoensis  yonahlossee      jordani   glutinosus 
##    1.0274635    1.0125997    1.0154203    0.9846055    0.9588364    1.0381764    1.0457046 
##    punctatus      wehrlei     websteri     dorsalis      welleri     serratus    hubrichti 
##    1.0411562    0.9987080    0.9933846    1.0023392    0.9525084    0.9857701    1.0076362 
##     cinereus     nettingi     hoffmani   shenandoah    richmondi 
##    1.0432823    0.9975689    1.0282926    1.0166452    1.0169282

(Please keep in mind – these are random!)

Now we’re pretty much ready to create our plot!

To start, however, we might graph our tree without all the trait info – just to gauge how much space we’re going to need for it!

plotTree(salamanders, ftype="i")

Now let’s use get("last_plot.phylo",envir=.PlotPhyloEnv) to pull down the information about this plot, including the x and y user unit dimensions. (This is the most important part for our purposes.)

pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
pp$x.lim

## [1]   0.0000 507.4171

pp$y.lim

## [1]  1 26

Now let’s simply re-plot our tree, updating these values, and add all of our desired elements on top of it. Since I have to do this in a single code chunk, I’ll do my best to explain in comments!

## graph the tree, updating xlim & ylim
plotTree(salamanders,ftype="off",
  xlim=c(pp$x.lim[1],pp$x.lim[2]*1.15),
  ylim=c(pp$y.lim[1],pp$y.lim[2]*1.15))
## pull the last plotted phylo info
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## compute the string width in user units of our
## labels
sw<-sapply(salamanders$tip.label,strwidth,font=3)
## right align & re-graph our labels (this just 
## looks nicer)
par(lend=1)
for(i in 1:Ntip(salamanders)){
  text(x=pp$xx[i]+1.05*max(sw),y=pp$yy[i],
    salamanders$tip.label[i],cex=0.9,font=3,
    pos=2)
  segments(pp$xx[i],pp$yy[i],pp$xx[i]+max(sw)-sw[i],
    pp$yy[i],lty="dotted")
}
## sort our lambda values & discretize them as 
## > 1 (increasing) or < 1 (decreasing)
lambda<-lambda[salamanders$tip.label]
foo<-function(x) if(x>1) "white" else "black"
lam_col<-sapply(lambda,foo)
## graph the levels using points (this could 
## also be done more scientifically with e.g.,
## plotrix::draw.circle)
## the specific coordinates might need adjustment
points(pp$xx[1:Ntip(salamanders)]+
    1.05*max(sw)+0.02*max(nodeHeights(salamanders)),
  pp$yy[1:Ntip(salamanders)],pch=21,bg=lam_col,
  cex=2.2)
## now add IUCN status as text. Once again, the 
## font size & offset might need adjustment 
## depending on the case
text(pp$xx[1:Ntip(salamanders)]+
    1.05*max(sw)+0.08*max(nodeHeights(salamanders)),
  pp$yy[1:Ntip(salamanders)],
  iucn[salamanders$tip.label])
## finall, column headers & legend
text(pp$xx[1:Ntip(salamanders)]+
    max(sw)+0.01*max(nodeHeights(salamanders)),
  max(pp$yy)+0.75,expression(paste(lambda," (growth)")),
  srt=80,pos=4)
text(pp$xx[1:Ntip(salamanders)]+
    max(sw)+0.07*max(nodeHeights(salamanders)),
  max(pp$yy)+0.75,"IUCN status",srt=80,
  pos=4)
legend("topleft",c("decreasing","increasing"),pch=21,
  pt.bg=c("black","white"),pt.cex=2,bty="n",inset=0.05)

I mean, we could’ve gotten fancier – but that’s not bad!