More on showing a barplot next to a plotted tree

I just added a new function to phytools - plotTree.barplot - to reproduce the visualization shown in a previous recent post.

The idea is to put a barplot next to a plotted tree, but in which the bars line up nicely with the tips.

There already exists a function in phytools to do this kind of visual. It is quite flexible, but it works by just drawing rectangles in the plotting device where we have already drawn a tree. To do this, the user needs to set a scale for the rectangles (the units of measure for our trait are not usually scaled to the edges of the phylogeny), and there is not an obvious way to include an axis.

Here is a demo of that function, using the data & tree that I employed in my prior recent demo:

library(phytools)
plotTree.wBars(eel.tree,bsize,scale=0.5,tip.label=TRUE,fsize=0.7)

The new implementation results in a highly similar plot - but instead of attempting to plot the bars & the tree in the same device, it opens a plot array using par()$mfrow and then prints the tree & bar plot to different subplots. This means that they are no longer required to be rendered to the same scale, and it is very easy to include an axis or other labels on the bar plot.

To accomplish this & still permit the user to pass whatever arguments they wanted to plotTree and barplot, I have made use of the function do.call, which executes a function with a list of arguments passed to it. That means that arguments that the user wants to send to each of these functions needs to be supplied in two lists: args.plotTree (to be passed to plotTree), and args.barplot.

First, here is the function code:

plotTree.barplot<-function(tree,x,args.plotTree=list(),
    args.barplot=list()){
    cw<-reorder(tree)
    args.barplot$height<-x[cw$tip.label]
    args.barplot$plot<-FALSE
    args.barplot$horiz<-TRUE
    args.barplot$axes<-FALSE
    args.barplot$names.arg<-""
    if(is.null(args.barplot$space)) args.barplot$space<-0.7
    if(is.null(args.barplot$mar)) 
        args.barplot$mar<-c(5.1,0,2.1,1.1)
    else args.barplot$mar[2]<-0.1
    args.plotTree$tips<-setNames(do.call(barplot,args.barplot)[,1],
        cw$tip.label)
    args.barplot$plot<-TRUE
    args.barplot$ylim<-range(args.plotTree$tips)
    args.plotTree$tree<-cw
    if(is.null(args.plotTree$mar)) 
        args.plotTree$mar<-c(5.1,1.1,2.1,0)
    else {
        args.plotTree$mar[4]<-0.1
    }
    if(args.plotTree$mar[1]!=args.barplot$mar[1])
        args.plotTree$mar[1]<-args.barplot$mar[1]
    if(args.plotTree$mar[3]!=args.barplot$mar[3])
        args.plotTree$mar[3]<-args.barplot$mar[3]
    if(is.null(args.plotTree$ftype)) args.plotTree$ftype<-"i"
    if(is.null(args.plotTree$lwd)) args.plotTree$lwd<-1
    par(mfrow=c(1,2))
    do.call(plotTree,args.plotTree)
    par(mar=args.barplot$mar)
    obj<-do.call(barplot,args.barplot)
    axis(1)
    if(!is.null(args.barplot$xlab)) title(xlab=args.barplot$xlab)
    else title(xlab="x")
    invisible(obj)
}

Now here is an example:

plotTree.barplot(eel.tree,bsize,list(fsize=0.7),
    list(col="blue",space=1,log="x",
    xlim=c(10,max(bsize)),xlab="length (cm)"))

You get the idea.

Showing a bar plot next to a plotted phylogeny

A phytools user asked me yesterday:

“I'm trying to plot a tree side-by-side with a barplot, however the bars are a bit off from the tip labels, how can I correct this?”

This is a great question, because it is very easy to see how this misalignment results:

library(phytools)
tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  R, E, S, W, X, Q, ...
## 
## Rooted; includes branch lengths.

x

##           A           B           C           D           E           F 
## -0.32084765  0.38196606 -0.12100495 -1.83763522 -0.92862392 -3.79335667 
##           G           H           I           J           K           L 
##  0.72976040  0.08573435 -0.88244919 -2.13005682 -2.75687573 -0.31006456 
##           M           N           O           P           Q           R 
## -1.52206369  1.36749732 -1.61131289  0.21897916 -1.44425879 -0.66191387 
##           S           T           U           V           W           X 
## -1.07884603 -0.99994662  0.49265512 -0.51313930 -0.82487486 -0.69957983 
##           Y           Z 
## -3.32350668 -3.86923149

par(mfrow=c(1,2))
plotTree(cw<-reorder(tree),mar=c(5.1,2.1,2.1,0.1))
par(mar=c(5.1,0.1,2.1,2.1))
barplot(x[cw$tip.label],horiz=TRUE,names.arg="")
title(xlab="log(body size)")

At first, this seems OK - but if we add some grid lines we can see that it is not quite….

par(mfrow=c(1,2))
plotTree(cw<-reorder(tree),mar=c(5.1,2.1,2.1,0.1))
nulo<-sapply(1:Ntip(cw),function(x) abline(h=x,col="grey"))
par(mar=c(5.1,0.1,2.1,2.1))
tmp<-barplot(x[cw$tip.label],horiz=TRUE,names.arg="")
title(xlab="log(body size)")
nulo<-sapply(tmp,function(x) abline(h=x,col="grey"))

The trick is to use the heights of the bars on the barplot to space our tip labels as follows:

cw<-reorder(tree)
tmp<-barplot(x[cw$tip.label],plot=FALSE)
par(mfrow=c(1,2))
plotTree(cw,tips=tmp,mar=c(5.1,1.1,2.1,0.1))
par(mar=c(5.1,0.1,2.1,1.1))
barplot(x[cw$tip.label],horiz=TRUE,axes=FALSE,
    ylim=range(tmp),names.arg="",xlim=range(x))
axis(1)
title(xlab="log(body size)")

and here with our grid lines:

cw<-reorder(tree)
tmp<-barplot(x[cw$tip.label],plot=FALSE)
par(mfrow=c(1,2))
plotTree(cw,tips=tmp,mar=c(5.1,1.1,2.1,0.1))
nulo<-sapply(tmp,function(x) abline(h=x,col="grey"))
par(mar=c(5.1,0.1,2.1,1.1))
barplot(x[cw$tip.label],horiz=TRUE,axes=FALSE,
    ylim=range(tmp),names.arg="",xlim=range(x))
axis(1)
title(xlab="log(body size)")
nulo<-sapply(tmp,function(x) abline(h=x,col="grey"))

Cool. Here it is applied to real data:

eel.tree<-read.tree(file=
    "http://www.phytools.org/SanJuan2016/data/elopomorph.tre")
eel.data<-read.csv(file=
    "http://www.phytools.org/SanJuan2016/data/elopomorph.csv",
    header=TRUE,row.names=1)
bsize<-setNames(eel.data[,2],rownames(eel.data))
cw<-reorder(eel.tree)
tmp<-barplot(bsize[cw$tip.label],plot=FALSE,space=0.7)
par(mfrow=c(1,2))
plotTree(cw,tips=tmp,mar=c(5.1,1.1,2.1,0.1),fsize=0.7,
    ftype="i",lwd=1)
par(mar=c(5.1,0.1,2.1,1.1))
barplot(bsize[cw$tip.label],horiz=TRUE,axes=FALSE,
    ylim=range(tmp),names.arg="",space=0.7)
axis(1)
title(xlab="body size (cm)")

That's it.