Several of the many ways you can plot two trees next to each other while maintaining the same scale

A friend & colleague today sent me an email with the following inquiry:

“I’ve been wrestling with what should be a relatively simple task: I want to plot 2 trees in vertically oriented panels, both on the same time scale. I put together a toy example, attached, using phytools. Clearly, it’s not working. It would seem that even though I set up the plot with the same x/y limits, and use plotTree with add = TRUE, it is not working as I expected to.”

Like many such problems, there’s more than one way to do this. I’m going to run through them, but first I’ll simulate two trees of different depths following my colleague’s example:

library(phytools)
tree1<-pbtree(n=50,scale=80)
tree2<-pbtree(n=50,scale=100)

OK, the first way (possibly the funnest & most elegant) to plot these two trees in a bipanel figure with the same horizontal axis is to plot them backwards (i.e., facing leftwards), but with the axis orientation flipped.

Here’s what that looks like:

par(mfrow=c(2,1))
plotTree(tree1,xlim=c(100,0),direction="leftwards",
  ftype="off",lwd=1,mar=c(2.1,0.6,0.6,0.6))
axis(1,cex.axis=0.8)
plotTree(tree2,xlim=c(100,0),direction="leftwards",
  ftype="off",lwd=1,mar=c(2.1,0.6,0.6,0.6))
axis(1,cex.axis=0.8)

One smallish catch with this approach is that it will make it more difficult to figure out how much space to leave for our labels, should we wish to include them in our plot!

To illustrate how this can be accomplished, I’m going to switch to using two real phylogenies for a change, each time-calibrated in the same units (millions of years), each with a different total depth, and each with different tip labels.

data(primate.tree) ## from Kirk & Kay (2004)
data(whale.tree) ## from Steeman et al. (2009)

Next, let’s pre-plot our trees, but just with their own separate scales. It’s totally fine to use different font sizes, but we must be consistent across this plot and the one to follow.

par(mfrow=c(2,1))
plotTree(primate.tree,direction="leftwards",
  ftype="i",fsize=0.3,lwd=1,mar=c(2.1,0.6,0.6,0.6))
axis(1,cex.axis=0.8)
x1<-get("last_plot.phylo",
  envir=.PlotPhyloEnv)$x.lim[1]
plotTree(whale.tree,direction="leftwards",ftype="i",
  fsize=0.3,lwd=1,mar=c(2.1,0.6,0.6,0.6))
axis(1,cex.axis=0.8)

x2<-get("last_plot.phylo",
  envir=.PlotPhyloEnv)$x.lim[1]

xlim<-c(1,min(c(x1/max(nodeHeights(primate.tree)),
    x2/max(nodeHeights(whale.tree)))))*max(c(
      nodeHeights(primate.tree),
      nodeHeights(whale.tree)))
xlim

## [1]  73.00302 -11.77711

par(mfrow=c(2,1))
plotTree(primate.tree,xlim=xlim,direction="leftwards",
  ftype="i",fsize=0.3,lwd=1,mar=c(2.1,0.6,0.6,0.6))
axis(1,at=seq(0,60,by=20),cex.axis=0.8)
plotTree(whale.tree,xlim=xlim,direction="leftwards",
  ftype="i",fsize=0.3,lwd=1,mar=c(2.1,0.6,0.6,0.6))
axis(1,at=seq(0,60,by=20),cex.axis=0.8)

If we want, we can also apply the same trick to plotting the phylogenies vertically. Here I’ll do that, but keep only one of my axes across the two subplots.

par(mfrow=c(1,2))
plotTree(primate.tree,ylim=xlim,
  direction="downwards",ftype="i",
  fsize=0.3,lwd=1,mar=c(0.1,2.1,1.1,0.1))
mtext("a) primate tree",3,line=0,
  at=0.1*Ntip(primate.tree))
axis(2,at=seq(0,60,by=20),cex.axis=0.8)
plotTree(whale.tree,
  ylim=xlim,direction="downwards",ftype="i",
  fsize=0.3,lwd=1,mar=c(0.1,0.1,1.1,2.1))
mtext("b) whale tree",3,line=0,
  at=0.1*Ntip(whale.tree))

That works.

Way #2 is to add a root edge to the shorter tree, convert this tree into a "phylo" object with unbranching nodes, map two separate “regimes” onto this tree, and the plot the root edge in a transparent color. This sounds complicated, but is actually pretty straightforward.

First, here with our simulated trees.

par(mfrow=c(2,1))
## tree1 is my shorter tree
re_tree1<-tree1
re_tree1$root.edge<-max(nodeHeights(tree2))-
  max(nodeHeights(tree1))
re_tree1<-rootedge.to.singleton(re_tree1)
re_tree1<-paintBranches(re_tree1,
  edge=Ntip(re_tree1)+2,state="0")
cols<-setNames(c("transparent","black"),0:1)
plot(re_tree1,cols,ftype="off",lwd=1,
  mar=c(2.1,0.6,0.6,0.6))
axis(1,cex.axis=0.8)
## tree2 is my longer tree
plotTree(tree2,ftype="off",lwd=1,
  mar=c(2.1,0.6,0.6,0.6))
axis(1,cex.axis=0.8)

This gets more complicated if tip labels are involved, but we can still use basically the same solution.

par(mfrow=c(2,1))
## primate.tree is my longer tree
plotTree(primate.tree,ftype="i",fsize=0.3,lwd=1,
  mar=c(2.1,0.6,0.6,0.6),plot=FALSE)
x1<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim[2]
## whale.tree is my shorter tree
re_whale.tree<-whale.tree
re_whale.tree$root.edge<-max(nodeHeights(primate.tree))-
  max(nodeHeights(whale.tree))
re_whale.tree<-rootedge.to.singleton(re_whale.tree)
re_whale.tree<-paintBranches(re_whale.tree,
  edge=Ntip(re_whale.tree)+2,state="0")
cols<-setNames(c("transparent","black"),0:1)
plotTree(re_whale.tree,ftype="i",fsize=0.3,lwd=1,
  mar=c(2.1,0.6,0.6,0.6),plot=FALSE,add=TRUE)
x2<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim[2]
plotTree(primate.tree,ftype="i",fsize=0.3,lwd=1,
  mar=c(2.1,0.6,0.6,0.6),xlim=c(0,max(x1,x2)),add=TRUE)
axis(1,at=max(nodeHeights(primate.tree))-seq(0,60,by=20),
  seq(0,60,by=20),cex.axis=0.8)
plot(re_whale.tree,cols,ftype="i",fsize=0.3,lwd=1,
  mar=c(2.1,0.6,0.6,0.6),xlim=c(0,max(x1,x2)))
axis(1,at=max(nodeHeights(re_whale.tree))-seq(0,60,by=20),
  seq(0,60,by=20),cex.axis=0.8)

Finally, way #3 involves first graphing the longer tree, then figuring out how the x limits need to be “fudged” so that the long and the short tree end at the same place. I’m going to demonstrate this just with our unlabeled trees, just to keep the complicated relatively simple.

par(mfrow=c(2,1))
## graph the shorter tree without plotting it
plotTree(tree1,plot=FALSE,ftype="off",
  mar=c(2.1,0.6,0.6,0.6))
x_diff<-diff(get("last_plot.phylo",
  envir=.PlotPhyloEnv)$x.lim)
## graph the shorter tree
plotTree(tree1,ftype="off",mar=c(2.1,0.6,0.6,0.6),
  xlim=c(x_diff-max(nodeHeights(tree2)),
    max(nodeHeights(tree1))),add=TRUE,lwd=1)
axis(1,at=max(nodeHeights(tree1))-seq(0,100,by=20),
  seq(0,100,by=20))
## graph the longer tree
plotTree(tree2,ftype="off",mar=c(2.1,0.6,0.6,0.6),
  lwd=1)
axis(1,at=max(nodeHeights(tree2))-seq(0,100,by=20),
  seq(0,100,by=20))

That’s it, folks.