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.

Stochastic mapping when tip states are uncertain - revisited

Stimulated by a user inquiry, here is a quick re-hash of ancestral character estimation using the phytools function make.simmap when tip states are uncertain.

Firstly, let's simulate some data with the property of uncertainty in the tip values - that is, we do not know the states for some extant taxa in the tree. The way this is expressed is as a prior probability distribution on the states for those tips that we do not know. For instance, if our character has two states "a" & "b", and we are completely ignorant of the values of certain species in the tree, then we might say that the prior probability of being in each of states "a" or "b" was 0.5.

## load phytools
library(phytools)
## simulate stochastic pure-birth tree
tree<-pbtree(n=26,tip.label=LETTERS,scale=1)
## generate character transition matrix
Q<-matrix(c(-1,1,1,-1),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
Q

##    a  b
## a -1  1
## b  1 -1

## simulate character
x<-sim.history(tree,Q)$states

## Done simulation(s).

x

##   A   B   C   D   E   F   G   H   I   J   K   L   M   N   O   P   Q   R 
## "b" "b" "b" "b" "b" "b" "b" "b" "b" "b" "b" "a" "a" "a" "b" "b" "b" "b" 
##   S   T   U   V   W   X   Y   Z 
## "b" "b" "b" "a" "a" "a" "a" "a"

Next, we have to add some uncertainty. Arbitrarily, let's say that we do not know the states of five of the twenty-six taxa in our tree. For the simulation, let's choose the five taxa with missing data at random:

## first encode our original data as a matrix
x<-to.matrix(x,seq=letters[1:2])
x

##   a b
## A 0 1
## B 0 1
## C 0 1
## D 0 1
## E 0 1
## F 0 1
## G 0 1
## H 0 1
## I 0 1
## J 0 1
## K 0 1
## L 1 0
## M 1 0
## N 1 0
## O 0 1
## P 0 1
## Q 0 1
## R 0 1
## S 0 1
## T 0 1
## U 0 1
## V 1 0
## W 1 0
## X 1 0
## Y 1 0
## Z 1 0

## now set some of these taxa to be uncertain
x[sample(1:26,5),]<-rep(0.5,2)
x

##     a   b
## A 0.0 1.0
## B 0.0 1.0
## C 0.0 1.0
## D 0.0 1.0
## E 0.0 1.0
## F 0.5 0.5
## G 0.0 1.0
## H 0.5 0.5
## I 0.0 1.0
## J 0.0 1.0
## K 0.0 1.0
## L 1.0 0.0
## M 1.0 0.0
## N 1.0 0.0
## O 0.0 1.0
## P 0.0 1.0
## Q 0.0 1.0
## R 0.5 0.5
## S 0.0 1.0
## T 0.5 0.5
## U 0.0 1.0
## V 1.0 0.0
## W 1.0 0.0
## X 0.5 0.5
## Y 1.0 0.0
## Z 1.0 0.0

Finally, let's generate some stochastic character maps using the phytools function make.simmap:

trees<-make.simmap(tree,x,nsim=100)

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##            a          b
## a -0.9951156  0.9951156
## b  0.9951156 -0.9951156

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##   a   b 
## 0.5 0.5

## Done.

trees

## 100 phylogenetic trees

If we want to obtain the posterior probabilities at nodes or tips, we can use the function describe.simmap:

obj<-describe.simmap(trees)
obj

## 100 trees with a mapped discrete character with states:
##  a, b 
## 
## trees have 7.6 changes between states on average
## 
## changes are of the following types:
##       a,b  b,a
## x->y 3.55 4.05
## 
## mean total time spent in each state is:
##              a         b    total
## raw  2.4524900 5.0739916 7.526482
## prop 0.3258481 0.6741519 1.000000

plot(obj)

The matrices of posterior probabilities at nodes and tips are stored in the two matrices obj$ace and obj$tips, respectively. Check it out:

obj$ace

##       a    b
## 27 0.35 0.65
## 28 0.07 0.93
## 29 0.05 0.95
## 30 0.05 0.95
## 31 0.13 0.87
## 32 0.11 0.89
## 33 0.00 1.00
## 34 0.00 1.00
## 35 0.00 1.00
## 36 0.00 1.00
## 37 0.00 1.00
## 38 0.00 1.00
## 39 0.00 1.00
## 40 0.00 1.00
## 41 0.59 0.41
## 42 1.00 0.00
## 43 0.00 1.00
## 44 0.07 0.93
## 45 0.49 0.51
## 46 0.27 0.73
## 47 0.88 0.12
## 48 0.89 0.11
## 49 0.83 0.17
## 50 0.95 0.05
## 51 1.00 0.00

obj$tips

##      a    b
## A 0.00 1.00
## B 0.00 1.00
## C 0.00 1.00
## D 0.00 1.00
## E 0.00 1.00
## F 0.09 0.91
## G 0.00 1.00
## H 0.01 0.99
## I 0.00 1.00
## J 0.00 1.00
## K 0.00 1.00
## L 1.00 0.00
## M 1.00 0.00
## N 1.00 0.00
## O 0.00 1.00
## P 0.00 1.00
## Q 0.00 1.00
## R 0.20 0.80
## S 0.00 1.00
## T 0.42 0.58
## U 0.00 1.00
## V 1.00 0.00
## W 1.00 0.00
## X 0.86 0.14
## Y 1.00 0.00
## Z 1.00 0.00

Of course, the posterior probabilities at the tips will normally be different than the prior probabilities, as we would expect.

Another possibility, of course, for visualizing the variability across maps at nodes & tips is the phytools function densityMap.

densityMap(trees,outline=TRUE)

## sorry - this might take a while; please be patient

Finally, if we didn't know that the function describe.simmap existed, it would not be too difficult to get these same values using the function getStates from phytools, as follows:

## first nodes
X<-getStates(trees)
levs<-sort(unique(as.vector(X)))
ace<-t(apply(X,1,function(x,l) summary(factor(x,levels=l)),l=levs))/ncol(X)
ace

##       a    b
## 27 0.35 0.65
## 28 0.07 0.93
## 29 0.05 0.95
## 30 0.05 0.95
## 31 0.13 0.87
## 32 0.11 0.89
## 33 0.00 1.00
## 34 0.00 1.00
## 35 0.00 1.00
## 36 0.00 1.00
## 37 0.00 1.00
## 38 0.00 1.00
## 39 0.00 1.00
## 40 0.00 1.00
## 41 0.59 0.41
## 42 1.00 0.00
## 43 0.00 1.00
## 44 0.07 0.93
## 45 0.49 0.51
## 46 0.27 0.73
## 47 0.88 0.12
## 48 0.89 0.11
## 49 0.83 0.17
## 50 0.95 0.05
## 51 1.00 0.00

## now tips
X<-getStates(trees,"tips")
levs<-sort(unique(as.vector(X)))
tips<-t(apply(X,1,function(x,l) summary(factor(x,levels=l)),l=levs))/ncol(X)
tips

##      a    b
## A 0.00 1.00
## B 0.00 1.00
## C 0.00 1.00
## D 0.00 1.00
## E 0.00 1.00
## F 0.09 0.91
## G 0.00 1.00
## H 0.01 0.99
## I 0.00 1.00
## J 0.00 1.00
## K 0.00 1.00
## L 1.00 0.00
## M 1.00 0.00
## N 1.00 0.00
## O 0.00 1.00
## P 0.00 1.00
## Q 0.00 1.00
## R 0.20 0.80
## S 0.00 1.00
## T 0.42 0.58
## U 0.00 1.00
## V 1.00 0.00
## W 1.00 0.00
## X 0.86 0.14
## Y 1.00 0.00
## Z 1.00 0.00

That's it.