Drawing boxes around clades graphed in a circular (fan or arc) style

In Chapter 13 of my 2022 book with Luke Harmon I provide a simple demo on how to draw a polygon around a clade of a plotted tree.

It goes something like this:

library(phytools)

## set seed
set.seed(99)
## simulate tree
tree<-pbtree(n=80)
tree

## 
## Phylogenetic tree with 80 tips and 79 internal nodes.
## 
## Tip labels:
##   t5, t73, t74, t38, t53, t54, ...
## 
## Rooted; includes branch length(s).

First, let’s draw the tree and figure out which clade we want to box.

plotTree(tree,ftype="off",direction="upwards")
nodelabels(bg="white",cex=0.6)

Let’s draw a polygon around the clade descended from node 124.

## open plotting & calculate tree coordinates
## but don't plot
plotTree(tree,ftype="off",direction="upwards",
  lwd=3,plot=FALSE)
## get coordinates of last plotted phylogeny
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## identify node of interest
node<-124
## find all the descendants that are tips
dd<-getDescendants(tree,node=node)
dd<-dd[dd<=Ntip(tree)]
## find the maximum height of the tree
max_h<-max(nodeHeights(tree))
## identify the upper & lower bounds of our
## box (with a buffer)
h1<-pp$yy[node]-0.02*max_h
h2<-max(pp$yy[dd])+0.02*max_h
## create the coordinates of our polygon
xx<-c(rep(min(pp$xx[dd])-0.5,2),
  rep(max(pp$xx[dd])+0.5,2))
yy<-c(h1,h2,h2,h1)
## graph the polygon
polygon(xx,yy,border=FALSE,col="grey")
## plot tree on top of polygon
plotTree(tree,ftype="off",direction="upwards",
  lwd=3,add=TRUE)
plotTree(tree,ftype="off",direction="upwards",
  lwd=1,add=TRUE,color="white")

In spite of the popularity of this graphing style, I’d never before posted on how to undertake a similar visualization for a “fan” or “arc” style tree. That’s simply because I’d never gone to the trouble of figuring it out! Let’s do that today.

We can start by plotting our tree in this arc style

part<-0.5 ## part of arc to use
plotTree(tree,type="arc",arc_height=1,ftype="off",
  part=part,lwd=1)
nodelabels(bg="white",cex=0.6)

Once again, we want to add our polygon around the clade descended from node 124.

## set part
part<-0.5
## identify node and its descendants
node<-124
dd<-getDescendants(tree,node)
dd<-dd[which(dd<=Ntip(tree))]
## open plotting device (plot=FALSE not available
## for this plotting style)
plotTree(tree,type="arc",arc_height=1,ftype="off",
  part=part,lwd=1,color="transparent")
## get coordinates of last plotted tree
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## compute the open angles of all the tips
theta<-atan(pp$yy[dd]/pp$xx[dd])
theta[pp$xx[dd]<0]<-pi+theta[pp$xx[dd]<0]
theta[pp$xx[dd]>0&pp$yy[dd]<0]<-
  2*pi+theta[pp$xx[dd]>0&pp$yy[dd]<0]
## identify the maximum & minimum open angles
max_theta<-max(theta)+0.5*2*pi/Ntip(tree)*part
min_theta<-min(theta)-0.5*2*pi/Ntip(tree)*part
## choose number of points for the arc
m<-30
## figure out the maximum height (for our buffer)
max_h<-max(sqrt(pp$xx^2+pp$yy^2))
## calculate the inner & outer bounds of box
h1<-sqrt(pp$xx[node]^2+pp$yy[node]^2)-0.02*max_h
h2<-max(sqrt(pp$xx[dd]^2+pp$yy[dd]^2))+0.02*max_h
## get the angles of the m points in our arc
tt<-seq(min_theta,max_theta,length.out=m)
## compute the points in our arc
out_arc.x<-h2*cos(tt)
out_arc.y<-h2*sin(tt)
in_arc.x<-h1*cos(tt)
in_arc.y<-h1*sin(tt)
xx<-c(out_arc.x,in_arc.x[m:1])
yy<-c(out_arc.y,in_arc.y[m:1])
## add the polygon
polygon(xx,yy,lwd=5,col=make.transparent("lightblue",1),
  border=FALSE)
## graph our tree (with an outline)
plotTree(tree,type="arc",arc_height=1,ftype="off",
  add=TRUE,lwd=4,part=part)
plotTree(tree,type="arc",arc_height=1,ftype="off",
  add=TRUE,lwd=2,color="white",part=part)

Let’s see if we can turn this into a general function….

cladebox<-function(node,col="#0000FF40",part=0.5,...){
  if(hasArg(npoints)) npoints<-list(...)$npoints
  else npoints<-30
  pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
  phy<-list(edge=pp$edge,Nnode=pp$Nnode,
    tip.label=1:pp$Ntip)
  class(phy)<-"phylo"
  dd<-getDescendants(phy,node)
  dd<-dd[which(dd<=Ntip(phy))]
  theta<-atan(pp$yy[dd]/pp$xx[dd])
  theta[pp$xx[dd]<0]<-pi+theta[pp$xx[dd]<0]
  theta[pp$xx[dd]>0&pp$yy[dd]<0]<-
    2*pi+theta[pp$xx[dd]>0&pp$yy[dd]<0]
  max_theta<-max(theta)+0.25*2*pi/Ntip(phy)*part
  min_theta<-min(theta)-0.25*2*pi/Ntip(phy)*part
  max_h<-max(sqrt(pp$xx^2+pp$yy^2))
  h1<-sqrt(pp$xx[node]^2+pp$yy[node]^2)-0.02*max_h
  h2<-max(sqrt(pp$xx[dd]^2+pp$yy[dd]^2))+0.02*max_h
  tt<-seq(min_theta,max_theta,length.out=npoints)
  out_arc.x<-h2*cos(tt)
  out_arc.y<-h2*sin(tt)
  in_arc.x<-h1*cos(tt)
  in_arc.y<-h1*sin(tt)
  xx<-c(out_arc.x,in_arc.x[npoints:1])
  yy<-c(out_arc.y,in_arc.y[npoints:1])
  polygon(xx,yy,lwd=5,col=col,border=FALSE)
}

plotTree(tree,type="arc",arc_height=1,ftype="off",
  lwd=1)
nodelabels(bg="white",cex=0.6)
cols<-RColorBrewer::brewer.pal(n=6,name="Pastel1")
cladebox(node=83,col=make.transparent(cols[1],0.4))
cladebox(node=108,col=make.transparent(cols[2],0.4))
cladebox(node=118,col=make.transparent(cols[3],0.4))
cladebox(node=123,col=make.transparent(cols[4],0.4))
cladebox(node=146,col=make.transparent(cols[5],0.4))
cladebox(node=152,col=make.transparent(cols[6],0.4))

OK. That’s not bad. Now let’s make one for Twitter.

par(fg="white",bg="black")
plotTree(tree,type="arc",arc_height=0.5,ftype="off",
  lwd=1,color="transparent")
cols<-RColorBrewer::brewer.pal(n=6,name="Dark2")
cladebox(node=83,col=cols[1])
cladebox(node=108,col=cols[2])
cladebox(node=118,col=cols[3])
cladebox(node=123,col=cols[4])
cladebox(node=146,col=cols[5])
cladebox(node=152,col=cols[6])
plotTree(tree,type="arc",arc_height=0.5,ftype="off",
  lwd=5,color="white",add=TRUE)
plotTree(tree,type="arc",arc_height=0.5,ftype="off",
  lwd=3,color="black",add=TRUE)