## Friday, January 8, 2016

### 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),
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)
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")
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)