Splitting a tree across multiple pages, with node labels & other such things

After posting a simple wrapper to split a plotted tree across plotting devices or pages of a PDF file, I received the following inquiry:

Would there be a way to modify this so that annotations like node labels and axes could be applied to both 1/2s of a split tree at the same time? I've been struggling with a similar problem using plot.phylo() and we can get our tree to look like what we want on one page, but splitting it up while keeping node and branch annotations etc is problematic.

This seems hard - but is actually not too difficult. I can just add an argument fn that is a function containing the functions we want to run for each page of our plotted tree.

First, here is the function:

split.plotTree<-function(tree,splits=NULL,file=NULL,fn,...){
    ef<-0.037037037037
    if(!is.null(file)) pdf(file,width=8.5,height=11)
    if(is.null(splits)) splits<-(floor(0.5*Ntip(tree))+0.5)/Ntip(tree)
    S<-matrix(c(0,splits,splits,1+1/Ntip(tree)),length(splits)+1,2)
    S<-cbind(S[,1]+ef*(S[,2]-S[,1]),S[,2]-ef*(S[,2]-S[,1]))
    for(i in nrow(S):1){
        if(is.null(file)&&i<nrow(S)) par(ask=TRUE)
        plotTree(tree,ylim=Ntip(tree)*S[i,],...)
        fn()
    }
    if(!is.null(file)) oo<-dev.off()
}

Now, just for example, let's create a custom function that adds an axis to each subplot, and includes node labels:

foo<-function(){
    nodelabels()
    axis(1,at=c(0,round(max(nodeHeights(tree)),2)))
}

Finally, let's load phytools, simulate a tree, and try it out:

library(phytools)
tree<-pbtree(n=90)
split.plotTree(tree,splits=c(0.3275,0.683),fn=foo,ftype="i",
    mar=c(3,1,1,1))

That's all there is to it.

Splitting a tree over mutiple plotting devices or pages

An R-sig-phylo subscriber asked the following yesterday:

I'm searching for a way to plot a huge phylogenetic tree to multiple pages in one searchable pdf file from R…. Does anybody know if this is possible and an effective way of doing this?

Well, my first response was that phytools has a function, splitplotTree, that can be used to split a plotted tree across columns in a single plotting window, or across plotting devices. Consequently, I responded:

One option is to use the function splitplotTree in phytools. By default it plots the split tree in two columns, but there is an option to split it into two windows. Then each plotted tree could be saved as a PDF.

On further thought, I realized that there is a better way. plotSimmap in phytools (and, consequently, plotTree, which uses plotTree, which uses plotSimmap internally) allows the user to control the y-limits of the plot. As a result it would be straightforward to write a simple wrapper that split the tree using a moving window over the plotting area. It also allows us to relatively easily split our tree over more than two pages. The following is an example of what that might look like:

split.plotTree<-function(tree,splits=NULL,file=NULL,...){
    ef<-0.037037037037
    if(!is.null(file)) pdf(file,width=8.5,height=11)
    if(is.null(splits)) splits<-(floor(0.5*Ntip(tree))+0.5)/Ntip(tree)
    S<-matrix(c(0,splits,splits,1+1/Ntip(tree)),length(splits)+1,2)
    S<-cbind(S[,1]+ef*(S[,2]-S[,1]),S[,2]-ef*(S[,2]-S[,1]))
    for(i in nrow(S):1){
        if(is.null(file)&&i<nrow(S)) par(ask=TRUE)
        plotTree(tree,ylim=Ntip(tree)*S[i,],...)
    }
    if(!is.null(file)) oo<-dev.off()
}

Note that in order to avoid having the tips & edges near the bottom of each page plot at the top of the following page, I had to trim ~4% off of the top & bottom of each plot. This is because for a given y range, R will create a plotting area that is 4% on each size larger than the area specified.

OK, now let's try it:

library(phytools)
n<-100
## create realistic tip labels
foo<-function(i) paste(sample(LETTERS,1),"._",paste(sample(letters,round(runif(1,min=4,max=8))),
    collapse=""),sep="")
tree<-pbtree(n=n,tip.label=sapply(1:n,foo))
splits<-c(0.255,0.505,0.755)
split.plotTree(tree,splits,ftype="i",mar=rep(1.1,4),fsize=0.9,lwd=1)

The splits mark the proportions of the plotted graph that I want to be displayed on each page.

We can also send the output to a multi-page PDF:

split.plotTree(tree,splits,ftype="i",mar=rep(1,4),file="split.plotTree.pdf",lwd=1)

The output is here: split.plotTree.pdf.

That's it.

Phylogenetic signal under Brownian evolution with bounds

A phytools user recently contacted me with some questions about her analyses and expressed some surprise at having discovered that bounds on the Brownian motion process tended to decrease phylogenetic signal measured using Blomberg et al.'s (2003) K statistic.

K is a measure of phylogenetic signal that is based on a standardized variance ratio - measuring the variance among vs. within clades compared to the ratio expected under Brownian motion. Consequently it has an expected value of 1.0, although the variance for a single simulation condition can be quite high.

Based simply on the logic of what K measures it shouldn't be too difficult to predict what happens to K when bounds are added to the Brownian process. That is, bounds will tend to cause the phenotypic values of species in unrelated clades to be more similar than expected under a pure Brownian process, making variability structured increasingly equally between & among clades. This effect should increase for decreasing bounds, as well as for increasing rate (σ² of the Broanian process), because both will cause evolving lineages to reach the boundaries more often.

This is also pretty easy to demonstrate via simulation. Here I have just simulated under a range of different conditions for the rate & bounds given a single tree, but the results hold generally.

## load packages
library(phytools)
## function to conduct simulation
foo<-function(tree,sig2,range,nsim){
    X<-fastBM(tree,a=range/2,sig2=sig2,bounds=c(0,range),nsim=nsim)
    mean(apply(X,2,phylosig,tree=tree))
}
## simulate tree
tree<-pbtree(n=26,tip.label=LETTERS,scale=1)
## set simulation conditions
sig2<-c(0.1,1:10)
range<-c(0.1,1:10)
nsim<-200
## simulate
K<-sapply(sig2,function(s2,tr,r,n) sapply(r,foo,tree=tr,sig2=s2,
    nsim=n),tr=tree,r=range,n=nsim)
colnames(K)<-sig2
rownames(K)<-range
## plot the result
filled.contour(x=range,y=sig2,K,xlab="range",ylab="rate",
    zlim=c(0,max(K)),color.palette=terrain.colors,
    main="Phylogenetic signal (K) for bounded BM")

Here we can see easily that for increasing rate (vertical axis) and decreasing bounds (horizontal axis), phylogenetic signal decreases - just as we predicted.

Note that, just as we showed in Revell et al. (2008; Syst. Biol.) there is no relationship between σ² and phylogenetic signal, K for unbounded Brownian motion:

foo<-function(tree,sig2,nsim){
    X<-fastBM(tree,sig2=sig2,nsim=nsim)
    apply(X,2,phylosig,tree=tree)
}
sig2<-c(0.001,0.01,0.1,1.0,10,100)
K<-sapply(sig2,foo,tree=tree,nsim=100)
colnames(K)<-sig2
boxplot(K,xlab="rate",ylab="phylogenetic signal (K)",
    main="Phylogenetic signal (K) for unbounded BM")
lines(colMeans(K),lty="dashed")
points(colMeans(K),pch=24)

That's it.