Neat new function to color branches by trait value

I had a user request yesterday for a function that would plot edge colors by probability - for example, if we have a value on [0, 1] for each edge in the tree - how do we use a built-in color palette in R to plot the tree using these edges as our colors? I think that the idea is somewhat akin to what contMap or densityMap accomplish; however if all we want is to plot each edge a different color, the 'ape' function plot.phylo is just fine at doing that.

I wrote a new function, plotBranchbyTrait, which automates this - calling plot.phylo internally; but the central piece of code used by the function is relatively simple. It merely takes our probabilities or reconstructed trait values along branches of the tree and translates them to colors from a palette as follows (using, in this case, a blue→red color map):

cols<-rainbow(1000,start=0.7,end=0) # blue->red
if(is.null(xlims)) xlims<-range(x)+c(-tol,tol)
breaks<-0:1000/1000*(xlims[2]-xlims[1])+xlims[1]
whichColor<-function(p,cols,breaks){
  i<-1
  while(p>=breaks[i]&&p>breaks[i+1]) i<-i+1
  cols[i]
}
colors<-sapply(x,whichColor,cols=cols,breaks=breaks)

What we come out with is a vector of colors to use as input in the argument plot.phylo.

I also have another internally called function, add.color.map, which takes user input for location and then plots a translation map (that also serves as a scale bar, as in, for instance, densityMap). I'd never done this before, but getting an interactively user supplied position for the color map was pretty easy. Here is the code snippet that I used:

cat("Click where you want to draw the bar\n")
x<-unlist(locator(1))
y<-x[2]
x<-x[1]

This is based on something similar in the ape function add.scale.bar.

One warning to the user - the function hangs before printing the "Click where you want to draw the bar" and waits for you to supply the location before the prompt! Not sure how to fix this.

The function has three different modes. mode="edges" uses the user supplied trait value for all the edges in the tree. The input values should be in the row order of tree$edge. mode="tips" just takes tip values and it reconstructs ancestral states by using fastAnc to compute the ancestral node states, and the averages the tipward and rootward states for each edge to get the plotted color. mode="nodes" does the same averaging, but the tip and node values are user supplied.

Here's a quick and dirty demo:

> source("plotBranchbyTrait.R") # load source
> # simplest use: ancestral BM values
> tree<-pbtree(n=40)
> x<-fastBM(tree)
> plotBranchbyTrait(tree,x,method="tips")
Click where you want to draw the bar

> # now let's say we have probabilities for each edge
> pp<-fastBM(tree,bounds=c(0,1),internal=TRUE)
> p<-rowMeans(matrix(pp[tree$edge],nrow(tree$edge),2))
> plotBranchbyTrait(tree,p,xlims=c(0,1),palette="gray")
Click where you want to draw the bar

Other methods in this function should be fairly self explanatory, I hope.

That's it.