A function to graph discrete and continuous traits at the tips of an "arc" or "fan" style plotted tree using phytools

Earlier I posted a hacky solution to the question of how to visualize one or more discrete characters (or discrete & continuous characters) at the tips of a tree plotted in “arc” or “fan” style in R.

Predictably, digging through my blog I discovered some earlier, less hacky solutions to the same visualization function (e.g., here), though, so far as I can recall, I’ve never built any of these solutions as a function for phytools.

I’m not sure I’ll add it to the package, but here is an example of this solution built into a function, plotFanTree.wTraits. I have adapted it so that I can combine both discrete and continuous traits, and returns the color palettes to the user invisibly.

plotFanTree.wTraits<-function(tree,X,type=c("arc","fan"),
  ...){
  X<-if(is.vector(X)) as.matrix(X[tree$tip.label]) else 
    X[tree$tip.label,]
  h<-max(nodeHeights(tree))
  d<-min(ncol(X)*0.07*h,h)
  type<-type[1]
  if(!(type%in%c("arc","fan"))) type<-"fan"
  ftype<-if(hasArg(ftype)) list(...)$ftype else "i"
  fsize<-if(hasArg(fsize)) list(...)$fsize else 0.5
  part<-if(hasArg(part)) list(...)$part else 0.99
  arc_height<-if(hasArg(arc_height)) list(...)$arc_height else 
    0.7
  lwd<-if(hasArg(lwd)) list(...)$lwd else 4
  if(length(lwd)<ncol(X)) lwd<-rep(lwd,ncol(X))
  if(hasArg(colors)) colors<-list(...)$colors
  else {
    colors<-list()
    for(i in 1:ncol(X)){
      if(is.numeric(X[,i])){
        colors[[i]]<-setNames(hcl.colors(n=100),1:100)
      } else {
        if(!is.factor(X[,i])) X[,i]<-as.factor(X[,i])
        colors[[i]]<-setNames(
          palette.colors(n=length(levels(X[,i]))),
          levels(X[,i]))
      }
    }
  }
  tt<-tree
  tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]<-
    tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]+d
  plotTree(tt,type=type,ftype=ftype,fsize=fsize,
    part=part,color="transparent",
    arc_height=arc_height*h/max(nodeHeights(tt)))
  pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
  plotTree(tree,type=type,ftype="off",part=part,
    lwd=1,add=TRUE,xlim=pp$x.lim,ylim=pp$y.lim,
    arc_height=arc_height,ftype="off")
  pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
  par(lend=3)
  for(i in 1:ncol(X)){
    if(is.numeric(X[,i])){
      x_seq<-seq(min(X[,i]),max(X[,i]),length.out=100)
      x_ind<-sapply(X[,i],function(x,y) which.min((x-y)^2),
        y=x_seq)
      colors[[i]]<-colorRampPalette(colors[[i]])(n=100)
      cols<-colors[[i]][x_ind]
    } else {
      cols<-colors[[i]][X[tree$tip.label,i]]  
    }
    for(j in 1:Ntip(tree)){
      start<-if(pp$xx[j]>0) 
        (i-1)*(d/ncol(X))+(2/7)*(d/ncol(X)) else 
          -((i-1)*(d/ncol(X))+(2/7)*(d/ncol(X)))
      end<-if(pp$xx[j]>0) i*d/ncol(X) else -i*d/ncol(X)
      th<-atan(pp$yy[j]/pp$xx[j])
      segments(pp$xx[j]+start*cos(th),
        pp$yy[j]+start*sin(th),
        pp$xx[j]+end*cos(th),
        pp$yy[j]+end*sin(th),
        lwd=lwd[i],
        col=cols[j])
    }
  }
  invisible(colors)
}

Let’s try to apply it to the same liolaemid lizard dataset from Esquerre et al. (2018) that I posted about yesterday.

I’ll start by loading & cleaning up the data.

library(phytools)
data("liolaemid.data")
data("liolaemid.tree")
summary(geiger::name.check(liolaemid.tree,liolaemid.data))

## 1 taxon is present in the data but not the tree:
##     Ctenoblepharys_adspersa
## 
## To see complete list of mis-matched taxa, print object.

liolaemid.data[liolaemid.tree$tip.label,]->liolaemid.data

plotFanTree.wTraits(liolaemid.tree,liolaemid.data[,3:1],
  lwd=10)

Note that the argument lwd is the width of the plotted line segments of our trait – not the tree (which will always be plotted with lwd=1 in our code). There’s no magic number of lwd that will look “right.” We just have to adjust it and figure out what looks good in our specific plotting device.

OK. Let’s try to make this a bit cooler.

colors<-list(
  c("blue","white","red"),
  terrain.colors(n=10),
  setNames(c("#F0EAD6","#DF536B"),c("O","V")))
cols<-plotFanTree.wTraits(liolaemid.tree,
  liolaemid.data[,3:1],lwd=12,colors=colors,ftype="off")
legend(x=0,y=0.7*max(nodeHeights(liolaemid.tree)),
  names(colors[[3]]),lwd=8,col=colors[[3]],
  title="parity mode",bty="n",xjust=0.5,yjust=0.5)
add.color.bar(1.5*max(nodeHeights(liolaemid.tree)),cols[[2]],
  title="maximum altitude (m)",
  lims=range(liolaemid.data[,2]),digits=2,prompt=FALSE,
  x=-0.75*max(nodeHeights(liolaemid.tree)),
  y=0.2*max(nodeHeights(liolaemid.tree)),subtitle="",
  lwd=8,outline=FALSE)
add.color.bar(1.5*max(nodeHeights(liolaemid.tree)),cols[[1]],
  title="environmental temp.",
  lims=range(liolaemid.data[,3]),digits=2,prompt=FALSE,
  x=-0.75*max(nodeHeights(liolaemid.tree)),
  y=-0.15*max(nodeHeights(liolaemid.tree)),subtitle="",
  lwd=8,outline=FALSE)

That’s not bad.