Alternative density tree plotting method in phytools (that also allows stochastic map trees)

A few of days ago, an R-sig-phylo subscriber asked if there was a way to represent the distribution of trees in a single plot in a manner “graphical style similar to the phytools::fancyTree(type = "phenogram95”)“.

In fact, such a function (densitTree) exists in phangorn; however since the subscriber (Giulio Dalla Riva at the University of Canterbury) added the challenge "bonus point if anybody proposes a solution that works with plotSimmap for painted trees”, I thought it might be worth demonstrating how this can be adapted to work with phytools plotting functions such as plotTree or plotSimmap.

I sent some code to the list, but it is fairly straightforward to write a custom function that does this visualization:

make.transparent<-function(color,alpha){
    RGB<-col2rgb(color)[,1]/255
    rgb(RGB[1],RGB[2],RGB[3],alpha)
}

density.tree<-function(trees,colors="blue",alpha=NULL,method="plotTree",fix.depth=FALSE,...){
    N<-length(trees)
    if(hasArg(use.edge.length)) use.edge.length<-list(...)$use.edge.length
    else use.edge.length<-TRUE
    if(!use.edge.length) trees<-lapply(trees,compute.brlen)
    if(!fix.depth){
        h<-sapply(trees,function(x) max(nodeHeights(x)))
        ii<-order(h,decreasing=TRUE)
        trees<-trees[ii]
        h<-h[ii]
    }
    if(is.null(alpha)) alpha<-1/N
    colors<-setNames(sapply(colors,make.transparent,alpha),names(colors))
    if(method=="plotTree") foo<-plotTree else foo<-plot
    foo(trees[[1]],color=colors,...)
    xlim<-get("last_plot.phylo",envir=.PlotPhyloEnv)$x.lim
    xlim[1]<-xlim[1]+0.03703704*diff(xlim)
    xlim[2]<-xlim[2]-0.03703704*diff(xlim)
    par(fg="transparent")
    for(i in 2:length(trees))
        foo(trees[[i]],
            tips=setNames(1:Ntip(trees[[1]]),trees[[1]]$tip.label),
            color=colors,add=TRUE,
            xlim=if(fix.depth) NULL else xlim-(h[1]-h[i]),...)
    par(fg="black")
}

First, let's use the function to do a basic visualization of topological and branch length uncertainty:

library(phytools)
trees<-read.nexus("Wayqecha.tree")
density.tree(trees,lwd=3,ftype="i",fsize=0.9)

We have all the options of plotTree & plotSimmap at our disposal here, and I find this type of visualization more readable with a non-standard node position which in phytools is called nodes="inner":

density.tree(trees,lwd=3,ftype="i",fsize=0.9,nodes="inner")

In this visualization, trees have variable depth (because the input trees were variable in total depth, but we can also turn this off and plot the trees with a standardized depth:

density.tree(trees,lwd=3,colors="black",ftype="i",fsize=0.9,nodes="inner",fix.depth=TRUE)

Finally, we can also turn off edge lengths entirely for this basic visualization by setting use.edge.length=FALSE:

density.tree(trees,lwd=3,ftype="i",fsize=0.9,nodes="inner",use.edge.length=FALSE)

As noted above, it is fairly straightforward to switch the plotting method used internally by this function to plotSimmap and thus plot stochastic map-style painted trees. For this example I will just use some simulated data:

Q<-matrix(c(-0.02,0.02,0.02,-0.02),2,2)
rownames(Q)<-colnames(Q)<-letters[1:2]
x<-sim.history(trees[[1]],Q)$states

## Done simulation(s).

mtrees<-make.simmap(trees,x,message=FALSE)
mtrees

## 100 phylogenetic trees with mapped discrete characters

density.tree(mtrees,colors=setNames(c("red","blue"),sort(unique(x))),
    nodes="inner",ftype="i",lwd=5,method="plotSimmap")
add.simmap.legend(colors=setNames(c("red","blue"),sort(unique(x))),
    prompt=FALSE,x=par()$usr[1]+8,y=par()$usr[4]-1)

## or
density.tree(mtrees,colors=setNames(c("red","blue"),sort(unique(x))),
    nodes="inner",ftype="i",lwd=5,method="plotSimmap",fix.depth=TRUE)
add.simmap.legend(colors=setNames(c("red","blue"),sort(unique(x))),
    prompt=FALSE,x=par()$usr[1]+1,y=par()$usr[4]-1)

That's it really.

One feature that would be good to include is a preliminary calculation of the tip order in a consensus tree. Presently, the tip order is the order of the deepest tree - which is fine when the trees vary little in topology, as in this case, but will work poorly if that tree happens to be quite topologically different from the rest of the trees in the sample.