Saturday, May 30, 2015

Bug fix for ltt with node labels

Today, Luke Harmon & Josef Uyeda reported a bug in the phytools function ltt, which does lineage-through-time plots. The difference between ltt and equivalent functions in other packages is that ltt permits trees to include lineages that end prior to the end of the tree (for instance, a tree with some lineages that have gone extinct). For example:

library(phytools)
tree<-pbtree(n=100,b=1.6,d=0.6)
plotTree(tree,ftype="off")

plot of chunk unnamed-chunk-1

obj<-ltt(tree)

plot of chunk unnamed-chunk-2

The reported bug is for conditions when the tree includes node labels. The source of this issue is that for trees in which all tips are contemporaneous (i.e., ultrametric trees) ltt uses the ape function branching.times internally and assumes that the names of the vector returned by branching.times are the node indices from the "phylo" object. They are if the tree lacks node labels; however if node labels are present, then they will instead contain the node labels of the tree. This creates an error when ltt looks for the node indices & doesn't find them! Here's a demo:

tree<-pbtree(n=26,tip.label=LETTERS,scale=10)
plotTree(tree)

plot of chunk unnamed-chunk-3

## no node labels
ltt(tree)

plot of chunk unnamed-chunk-4

## $ltt
##    27 31 39 32 28 40 33 48 35 45 49 37 51 41 46 43 50 47 29 38 44 30 34 36 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 
## 42    
## 26 26 
## 
## $times
##                        27           31           39           32 
## 0.000000e+00 9.999983e-10 1.594837e+00 1.613968e+00 1.809402e+00 
##           28           40           33           48           35 
## 2.384807e+00 4.478021e+00 5.274344e+00 5.812117e+00 6.371105e+00 
##           45           49           37           51           41 
## 6.434691e+00 7.038966e+00 7.336978e+00 7.372788e+00 7.762271e+00 
##           46           43           50           47           29 
## 8.524592e+00 8.967709e+00 9.034180e+00 9.082955e+00 9.457377e+00 
##           38           44           30           34           36 
## 9.480731e+00 9.539301e+00 9.748171e+00 9.763473e+00 9.872351e+00 
##           42              
## 9.946383e+00 1.000000e+01 
## 
## $gamma
## [1] 0.8561059
## 
## $p
## [1] 0.3919392
## add node labels
tree$node.label<-letters[1:25]
plotTree(tree)
nodelabels(tree$node.label)

plot of chunk unnamed-chunk-5

ltt(tree)
## Warning in ltt(tree): NAs introduced by coercion
## Error in 2:n: NA/NaN argument

plot of chunk unnamed-chunk-6

My bug fix, which is a little hacky to be fair, I first copy the node labels from the tree into a vector with names corresponding to the node indices, then I strip them, then I replace the names in the object returned by ltt with the original labels. The code is here. The following is a demo in which I load the function from source:

source("http://www.phytools.org/ltt/v1.0/ltt.R")
ltt(tree)

plot of chunk unnamed-chunk-7

## $ltt
##     a  e  m  f  b  n  g  v  i  s  w  k  y  o  t  q  x  u  c  l  r  d  h  j 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 
##  p    
## 26 26 
## 
## $times
##                         a            e            m            f 
## 0.000000e+00 9.999983e-10 1.594837e+00 1.613968e+00 1.809402e+00 
##            b            n            g            v            i 
## 2.384807e+00 4.478021e+00 5.274344e+00 5.812117e+00 6.371105e+00 
##            s            w            k            y            o 
## 6.434691e+00 7.038966e+00 7.336978e+00 7.372788e+00 7.762271e+00 
##            t            q            x            u            c 
## 8.524592e+00 8.967709e+00 9.034180e+00 9.082955e+00 9.457377e+00 
##            l            r            d            h            j 
## 9.480731e+00 9.539301e+00 9.748171e+00 9.763473e+00 9.872351e+00 
##            p              
## 9.946383e+00 1.000000e+01 
## 
## $gamma
## [1] 0.8561059
## 
## $p
## [1] 0.3919392

This update should be in the next version of phytools.

That's it.

No comments:

Post a Comment

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