A simple (too simple?) method of graphing a tree with an x-axis split

A phytools user recently asked about graphing a co-phylogenetic plot (“tanglegram”) in which one of the two trees had an axis break. E.g., here is his hand-drawn demo (very helpful):

Well, I’m not doing that yet; however, I did think of a simple (perhaps, too simple) way to create an axis break for a standard, right-facing phylogram plotted alone.

My solution involves using phytools::make.era.map to map discrete, temporal “regimes” onto the tree (and thus create a "simmap" object class), to shrink the “broken” part of the tree, and then to plot the break using a transparent color. (If our graphic device doesn’t support transparency, we can switch to the background color instead.)

Here’s a quick example in which I graph the intermediates as well as the final product.

## load packages
library(phytools)

## load tree
data(salamanders)
plotTree(salamanders,ftype="i",ylim=c(0,Ntip(salamanders)),
  mar=c(2.1,1.1,1.1,1.1))
axis(1,pos=0)
## set start & end points of the break
## (time since the root)
start<-60
end<-140
polygon(c(start,end,end,start),
  par()$usr[c(3,3,4,4)],col=make.transparent("grey",0.5),
  border=FALSE)

The grey box shows the block that I want to slice out.

So far so good. Let’s do it.

## make the "simmap" obect
tt<-make.era.map(salamanders,limits=c(0,start,end))
cols<-setNames(c("black","transparent","black"),1:3)
plot(tt,cols,ftype="i",ylim=c(0,Ntip(salamanders)),
  mar=c(2.1,1.1,1.1,1.1))
axis(1,pos=0)

Now let’s set our “buffer” width for the break. This is how big we want the split in our axis to look in the units of the branch lengths of our tree.

H<-max(nodeHeights(salamanders))
buffer<-0.03*(H-end+start)

Next, I’m going to manually rescale each of our three painted regime eras – the first & the third by 1.0 (so, not at all), and the middle one to have a total width equal to the break we want in the axis of our graph. This can be done as follows.

scale<-setNames(c(1,buffer/(end-start),1),1:3)
tt$maps<-lapply(tt$maps,function(x,scale) 
  x*scale[names(x)],scale=scale)
tt$edge.length<-sapply(tt$maps,sum)

Finally, let’s create our plot. I’m going to use phytools::markChanges to add a bit of a visual element to the break (although I could’ve also done this manually using lines or segments). The hardest part of this is figure out the labels for our x-axis line, particularly if we want to plot it from the present backwards.

plot(tt,cols,mar=c(2.1,1.1,1.1,1.1),ylim=c(0,Ntip(tt)),
  ftype="i",fsize=0.8)
markChanges(tt,setNames(rep("black",3),1:3))
labs<-pretty(c(0,max(nodeHeights(salamanders))),n=8)
labs<-labs[labs<max(nodeHeights(salamanders))]
at<-max(nodeHeights(salamanders))-labs
labs<-labs[-intersect(which(at>start),which(at<end))]
at<-at[-intersect(which(at>start),which(at<end))]
at[at>end]<-at[at>end]-(end-start)+buffer
axis(1,at=at,labels=labs,pos=0,cex.axis=0.8)
polygon(start+c(0,buffer,buffer,0),
  0.5*mean(strheight(LETTERS))*c(-1,-1,1,1),
  col=if(par()$bg=="transparent") "white" else par()$bg,
  border=FALSE)
segments(start+c(0,buffer),0.5*mean(strheight(LETTERS))*c(-1,-1),
  start+c(0,buffer),0.5*mean(strheight(LETTERS))*c(1,1),lwd=2)

Let’s put it in a function.

splitAxisTree<-function(tree,start,end,buffer=0.03,...){
  tt<-make.era.map(tree,limits=c(0,start,end))
  H<-max(nodeHeights(tree))
  buffer<-buffer*(H-end+start)
  scale<-setNames(c(1,buffer/(end-start),1),1:3)
  tt$maps<-lapply(tt$maps,function(x,scale) 
    x*scale[names(x)],scale=scale)
  tt$edge.length<-sapply(tt$maps,sum)
  args<-list(...)
  args$x<-tt
  args$colors<-cols
  args$direction<-"rightwards"
  if(is.null(args$mar)) args$mar<-c(2.1,1.1,1.1,1.1)
  if(is.null(args$ylim)) args$ylim<-c(0,Ntip(tt))
  if(is.null(args$ftype)) args$ftype<-"i"
  do.call(plot,args)
  markChanges(tt,setNames(rep("black",3),1:3))
  ## I'm doing some weird stuff here in case we're
  ## cutting out a lot of our tree!
  labs<-pretty(c(0,max(nodeHeights(tt))),n=8)
  labs<-seq(0,max(nodeHeights(tree)),by=labs[2]-labs[1])
  labs<-labs[labs<max(nodeHeights(tree))]
  at<-max(nodeHeights(tree))-labs
  labs<-labs[-intersect(which(at>start),which(at<end))]
  at<-at[-intersect(which(at>start),which(at<end))]
  at[at>=end]<-at[at>=end]-(end-start)+buffer
  axis(1,at=at,labels=labs,pos=0,cex.axis=0.8)
  polygon(start+c(0,buffer,buffer,0),
    0.5*mean(strheight(LETTERS))*c(-1,-1,1,1),
    col=if(par()$bg=="transparent") "white" else 
      par()$bg,border=FALSE)
  segments(start+c(0,buffer),
    0.5*mean(strheight(LETTERS))*c(-1,-1),
    start+c(0,buffer),
    0.5*mean(strheight(LETTERS))*c(1,1),lwd=2)
}

data(vertebrate.tree)
plotTree(vertebrate.tree,mar=c(3.1,1.1,1.1,1.1),
  ftype="i")
axis(1,at=seq(0,473,by=50))
grid()

splitAxisTree(vertebrate.tree,205,360,buffer=0.02,
  fsize=0.8,lwd=1,mar=c(4.1,1.1,1.1,1.1))
mtext("time (mybp)",side=1,line=2,at=180)

I think that pretty much works. Let’s try it on a coalescent tree.

First, here’s our tree.

set.seed(4)
tree<-rcoal(n=50)
plotTree(tree,lwd=1,ftype="off",
  mar=c(2.1,1.1,1.1,1.1))
grid()
at<-pretty(nodeHeights(tree),n=8)
at<-at[at<max(nodeHeights(tree))]
axis(1,at=max(nodeHeights(tree))-at,labels=at)

h<-max(nodeHeights(tree))
splitAxisTree(tree,0.05,h-0.48,ftype="off",lwd=1)

OK. That’s pretty much what I was going for.