## Thursday, March 2, 2017

### Function to plot a tree with error bars on divergence dates in R

Yesterday I posted about plotting error bars around divergence times onto internal nodes of the tree.

The following (`plotTree.errorbars`) is a function to automate this procedure, to be shortly added to phytools.

``````## plot tree with error bars around divergence times at nodes
## written by Liam J. Revell 2017
plotTree.errorbars<-function(tree,CI,...){
args<-list(...)
if(!is.null(args\$gridlines)){
gridlines<-args\$gridlines
args\$gridlines<-NULL
} else gridlines<-TRUE
if(is.null(args\$mar)) args\$mar<-c(4.1,1.1,1.1,1.1)
if(is.null(args\$ftype)) args\$ftype<-"i"
fsize<-if(!is.null(args\$fsize)) args\$fsize else 1
if(is.null(args\$direction)) args\$direction<-"leftwards"
if(!is.null(args\$bar.width)){
bar.width<-args\$bar.width
args\$bar.width<-NULL
} else bar.width<-11
if(!is.null(args\$cex)){
cex<-args\$cex
args\$cex<-NULL
} else cex<-1.2
if(!is.null(args\$bar.col)){
bar.col<-args\$bar.col
args\$bar.col<-NULL
} else bar.col<-"blue"
par(mar=args\$mar)
plot.new()
th<-max(nodeHeights(tree))
h<-max(th,max(CI))
if(is.null(args\$xlim)){
m<-min(min(nodeHeights(tree)),min(CI))
d<-diff(c(m,h))
pp<-par("pin")[1]
sw<-fsize*(max(strwidth(tree\$tip.label,units="inches")))+
1.37*fsize*strwidth("W",units="inches")
alp<-optimize(function(a,d,sw,pp) (a*1.04*d+sw-pp)^2,
d=d,sw=sw,pp=pp,
interval=c(0,1e6))\$minimum
args\$xlim<-if(args\$direction=="leftwards") c(h,m-sw/alp) else
c(m,h+sw/alp)
}
if(is.null(args\$at)) at<-seq(0,h,by=h/5)
else {
at<-args\$at
args\$at<-NULL
}
args\$tree<-tree
do.call(plotTree,args=args)
if(gridlines) abline(v=at,lty="dashed",
col=make.transparent("grey",0.5))
axis(1,at=at,labels=signif(at,3))
obj<-get("last_plot.phylo",envir=.PlotPhyloEnv)
for(i in 1:tree\$Nnode+Ntip(tree))
lines(x=c(CI[i-Ntip(tree),1],CI[i-Ntip(tree),2]),
y=rep(obj\$yy[i],2),lwd=bar.width,lend=0,
col=make.transparent(bar.col,0.4))
points(obj\$xx[1:tree\$Nnode+Ntip(tree)],
obj\$yy[1:tree\$Nnode+Ntip(tree)],pch=19,col=bar.col,
cex=cex)
}
``````

As I have mentioned in an earlier post - one of the most difficult components of plotting a tree in R is leaving enough space in the plotting window for the tip labels. This is because we normally don't know how much space a particular string will occupy (in user units) until we have already opened our plotting device & set our `xlim` and `ylim` values. Consequently, the example will use a tree with tip labels `B` through `Z`, plus one really long tip label, to make sure that enough space has been allocated. BTW, my solution to this problem (above) is adapted from `plot.phylo` in the ape package.

OK, now let's test it out:

``````library(phytools)
tree
``````
``````##
## Phylogenetic tree with 26 tips and 25 internal nodes.
##
## Tip labels:
##  Z, Y, X, W, V, U, ...
##
## Rooted; includes branch lengths.
``````
``````CI
``````
``````##        lower     upper
## 27 2.8414221 2.2855101
## 28 2.4025799 1.7851208
## 29 1.1518398 0.6932472
## 30 1.1016629 0.5915718
## 31 1.5786906 0.9097206
## 32 1.2866289 0.6838407
## 33 0.4025250 0.0000000
## 34 0.7735031 0.2052436
## 35 0.6981268 0.2212514
## 36 2.7675974 2.1440829
## 37 2.3925982 1.8320681
## 38 1.5706797 0.9821414
## 39 1.2980457 0.9041184
## 40 1.1321683 0.5452341
## 41 0.4097613 0.0000000
## 42 0.4676039 0.0000000
## 43 0.7283202 0.2276139
## 44 0.4388847 0.0000000
## 45 1.2948954 0.7909557
## 46 0.6314838 0.1167126
## 47 1.2049871 0.6604358
## 48 1.0396180 0.5733247
## 49 0.5778952 0.1890986
## 50 0.3843571 0.0000000
## 51 1.0729443 0.5512091
``````
``````plotTree.errorbars(tree,CI)
``````

Note that there is a lot we can do to modify this plot. Firstly, any argument of `plotTree` can be passed internally to that function. For instance:

``````plotTree.errorbars(tree,CI,lwd=5,color="navy",lend=0)
``````

However, it is also possible to modify the specify attributes of the plot. For instance:

``````plotTree.errorbars(tree,CI,lwd=1,bar.col="red",bar.width=7,cex=1,
at=seq(0,3,by=0.5))
``````

Pretty neat.

The tree & CIs for this example were simulated as follows:

``````tree<-pbtree(n=26,tip.label=LETTERS[26:1])
tree\$tip.label[26]<-"Some really long tip label."
h<-sapply(1:tree\$Nnode+Ntip(tree),nodeheight,tree=tree)
d<-max(nodeHeights(tree))
CI<-cbind(h-runif(n=length(h),min=.10*d,max=.20*d),h+runif(n=length(h),
min=.05*d,max=.10*d))
CI[CI>max(nodeHeights(tree))]<-max(nodeHeights(tree))
CI<-max(nodeHeights(tree))-CI
rownames(CI)<-1:tree\$Nnode+Ntip(tree)
colnames(CI)<-c("lower","upper")
``````

#### 1 comment:

1. Hi Liam. How can I get a matrix of height_95%_HPD for all nodes from my phylogeny? When I read my nexus file in R doesn't read this values even though this value is in the file. So I'm looking for an easy method to extract these values to assemble this worksheet, but I do not find. Thank you for you help.

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.