Tuesday, June 7, 2016

Fix to contMap for short edges descended from the root node

A phytools user, Bruno Vilela, recently reported the bug that the phytools function contMap can fail for very short internal edges. This tends to be the case when an edge has a length that is less than 1/res × the tree height, in which res is an argument of contMap specifying the graininess of the color gradient to be plotted along edges.

[Actually, I believe I already fixed this bug except for edges descended from the root node. This is because for those edges the function will attempt to create a list of mappings along edges, but create a vector instead. Silly!]

I have previously suggested collapsing very short internal edges or increasing res; however Bruno did one better by providing an incredibly simple solution, which was to create a list of mappings along edges before iterating across the edges of the tree to compute the mappings. Perhaps obvious in retrospect, please keep in mind that I wrote the first version of this function in 2012 when I was quite an R programming novice. (I'm not sure what I'd call myself now. Slightly more than an a novice? I guess in an era when Donald Trump is running for present I'm a world-class R programming expert!)

Here's a quick demo of the problem, along with a demonstration that Bruno's fix does indeed restore function:

library(phytools)
tree
## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  t5, t26, t18, t17, t19, t22, ...
## 
## Rooted; includes branch lengths.
x
##          t5         t26         t18         t17         t19         t22 
##  0.20063042  0.90298749  2.28555569  2.11331867  1.38012423  0.01110355 
##         t15         t12          t6          t3          t2         t24 
##  0.31029532 -0.31270146  0.36006396  0.86328160 -0.16281977 -1.88665386 
##          t4          t9         t25         t11          t8         t21 
## -2.40507772 -2.69645998 -1.57485299  1.30619095  1.96062919 -0.91760839 
##          t1         t14          t7         t16         t20         t23 
## -1.45220994 -2.58959187 -1.44919866  0.93419977  0.19530914  0.06853757 
##         t13         t10 
##  1.79467390  0.29857041
plotTree(tree) ## short node from the root

plot of chunk unnamed-chunk-1

obj<-contMap(tree,x,plot=FALSE)
## Error in tree$maps[[i]] <- XX[, 2] - XX[, 1]: more elements supplied than there are to replace
## increase res
obj<-contMap(tree,x,plot=FALSE,res=1000)
plot(obj)

plot of chunk unnamed-chunk-1

Now let's load the new version from source off GitHub:

source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/contMap.R")
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/densityMap.R")
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/read.simmap.R")
obj<-contMap(tree,x,plot=FALSE)
plot(obj)

plot of chunk unnamed-chunk-2

This fix also allows us to set the value of res to any arbitrarily low value. Although the resultant plot is grainier, this might nonetheless be useful for users working with very large trees where computation is an issue.

contMap(tree,x,res=100)

plot of chunk unnamed-chunk-3

contMap(tree,x,res=10)

plot of chunk unnamed-chunk-3

Thanks Bruno!

No comments:

Post a Comment