Learning about the structure of the "phylo" object class through a phylogenetic Hanukkah menorah

Several years ago I posted a phylogenetic Christmas tree, adapted from one cleverly created by Gustavo Ballen.

This year, Christmas has passed – but Hanukkah lasts until the 2nd of January, 2025. One of the symbols of Hanukkah is the the menorah. As a gentile, I profess to be ignorant of its cultural & religious significance – however, it does look suspiciously like a phylogeny!

Here, for example, is a menorah for sale on Amazon:

To recreate this as a "phylo" object, I’m first going to re-draw it as follows.

Then, I’ll number the tips 1 through 9:

Finally, I’ll number the nodes 10 through 14, as follows:

Having drawn it out, I can easily re-write my menorah “tree” as an edge matrix, as follows.

edge<-matrix(c(
  10,11,
  11,1,
  11,12,
  11,9,
  12,2,
  12,13,
  12,8,
  13,3,
  13,14,
  13,7,
  14,4,
  14,5,
  14,6),13,2,byrow=TRUE)

To create a complete "phylo" object, we just need to combine this with a vector of tip labels, and a numeric value giving the number of nodes, as follows:

phy<-list(
  edge=edge,
  tip.label=1:9,
  Nnode=5
)

Then assign a class attribute:

class(phy)<-c("menorah","phylo")

(There’s no "menorah" class methods. I just did that for fun.)

Let’s load phytools and see if we’re on the right track.

library(phytools)

phy

## 
## Phylogenetic tree with 9 tips and 5 internal nodes.
## 
## Tip labels:
##   1, 2, 3, 4, 5, 6, ...
## 
## Rooted; no branch length.

par(mar=rep(2.1,4))
plot(phy,direction="upwards",show.tip.label=FALSE)
nodelabels()
tiplabels()

(We should be too alarmed that our “stem” branch is not plotted. Since we didn’t assign it a length, ape::plot.phylo uses compute.brlen internally, which assigns a length of zero to unbranching stem edges.)

Let’s assign our own edge lengths to (more or less) match the menorah from Amazon. These will take the form of a vector in the order of our rows of phy$edge, as follows.

phy$edge

##       [,1] [,2]
##  [1,]   10   11
##  [2,]   11    1
##  [3,]   11   12
##  [4,]   11    9
##  [5,]   12    2
##  [6,]   12   13
##  [7,]   12    8
##  [8,]   13    3
##  [9,]   13   14
## [10,]   13    7
## [11,]   14    4
## [12,]   14    5
## [13,]   14    6

phy$edge.length<-c(0.2,0.5,0.1,0.5,0.4,0.1,0.4,0.3,0.1,0.3,
  0.2,0.3,0.2)
plotTree(phy,direction="upwards",lwd=4)

Great. We’re totally on the right track. I think we can do better. How about this:

sigmoidPhylogram(phy,direction="upwards",lwd=20)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
sigmoidPhylogram(phy,direction="upwards",lwd=16,
  color="#C0C0C0",add=TRUE,ftype="off",ylim=pp$y.lim)

Finally, let’s get fancy. (Edit: updated to reflect the correct lighting sequence for the 8th night of Hanukkah, as I understand it.)

sigmoidPhylogram(phy,direction="upwards",lwd=14,
  ylim=c(0,0.92),ftype="off",color="blue",m1=0.06,b=5)
pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
## star of David
polygon(x=pp$xx[5]+c(0,0.5,-0.5),y=pp$yy[5]-c(0,0.17,0.17),
  border="transparent",col="blue")
polygon(x=pp$xx[5]+c(0,0.5,-0.5),y=pp$yy[5]-c(0.22,0.07,0.07),
  border="transparent",col="blue")
polygon(x=pp$xx[5]+c(0,0.5,-0.5),y=pp$yy[5]-c(0,0.17,0.17),
  border="grey")
polygon(x=pp$xx[5]+c(0,0.5,-0.5),y=pp$yy[5]-c(0.22,0.07,0.07),
  border="grey")
## candle bases
for(i in 1:Ntip(phy)){
  polygon(x=pp$xx[i]+c(-0.15,0.15,0.15,-0.15),
    y=pp$yy[i]+c(0,0,0.05,0.05),col="#C0C0C0")
}
## flames
x1<-c(0,1,2,3)/30
x2<-c(0,0.5,1.5,2.5,3)/30
y1<-(c(2,2.5,2.1,2.3)-1.6)/4
y2<-(c(1.2,0.9,0.8,1.5,2.3)-1.6)/4
f1<-smooth.spline(x1,y1,df=4)
left<-predict(f1,x=seq(0,max(x1),length.out=20))
f2<-smooth.spline(x2,y2,df=4)
right<-predict(f2,x=seq(0,max(x1),length.out=20))
xx<-c(left$y,right$y[20:1])
yy<-c(left$x,right$x[20:1])
lighting_sequence<-c(5,1,2,3,4,6,7,8,9)
for(i in lighting_sequence){
  polygon(x=pp$xx[i]+xx,y=pp$yy[i]+0.05+yy,
    col="gold")
}
## base
polygon(x=c(3.75,6.25,5.5,5.25,4.75,4.5),
  y=c(-0.03,-0.03,0.025,0.1,0.1,0.025),
  border="transparent",col="blue")

(Just don’t mess with the aspect ratio, or all bets are off!)

As a side lesson on how S3 classes work, let’s combine all of this into a class method for the object class "menorah".

plot.menorah<-function(phy,...){
  phy<-if(is.null(attr(phy,"order"))) reorder(phy) else phy
  ## rescale our tree & add unbranching root (if necessary)
  if(sum(phy$edge[,1]==(Ntip(phy)+1))>1){
    phy$edge.length<-phy$edge.length/max(nodeHeights(phy))*0.6
    phy$root.edge<-0.2
    phy<-rootedge.to.singleton(phy)
  } else {
    phy$edge.length<-phy$edge.length/max(nodeHeights(phy))*0.8
  }
  ## plot in sigmoid style
  sigmoidPhylogram(phy,direction="upwards",lwd=14,
    ylim=c(0,0.92),ftype="off",color="blue",m1=0.06,b=5)
  pp<-get("last_plot.phylo",envir=.PlotPhyloEnv)
  ## add candle bases
  for(i in 1:Ntip(phy)){
    polygon(x=pp$xx[i]+c(-0.15,0.15,0.15,-0.15),
      y=pp$yy[i]+c(0,0,0.05,0.05),col="#C0C0C0")
  }
  ## add flames
  x1<-c(0,1,2,3)/30
  x2<-c(0,0.5,1.5,2.5,3)/30
  y1<-(c(2,2.5,2.1,2.3)-1.6)/4
  y2<-(c(1.2,0.9,0.8,1.5,2.3)-1.6)/4
  f1<-smooth.spline(x1,y1,df=4)
  left<-predict(f1,x=seq(0,max(x1),length.out=20))
  f2<-smooth.spline(x2,y2,df=4)
  right<-predict(f2,x=seq(0,max(x1),length.out=20))
  xx<-c(left$y,right$y[20:1])
  yy<-c(left$x,right$x[20:1])
  for(i in 1:Ntip(phy)){
    polygon(x=pp$xx[i]+xx,y=pp$yy[i]+0.05+yy,
      col="gold")
  }
  ## base
  polygon(x=c(-1.25,1.25,0.5,0.25,-0.25,-0.5)+pp$xx[Ntip(phy)+1],
    y=c(-0.03,-0.03,0.025,0.1,0.1,0.025),
    border="transparent",col="blue")
}

Now let’s load a random tree, update the class, and test it. (Once again, I haven’t done anything to ensure this is insensitive to our aspect ratio, so results may vary!)

data(sunfish.tree)
sunfish.tree<-as.phylo(sunfish.tree)
class(sunfish.tree)<-c("menorah",class(sunfish.tree))
plot(sunfish.tree)

Happy Hanukkah!