Sunday, February 15, 2015

Perhaps the last (wishful thinking?) bug fix for reroot function for trees with node labels

I have just posted yet another bug fix for the phytools function reroot for trees containing node labels. This bug (reported by a helpful user) had the effect of causing one of the edge lengths descended from the new root to be wrong.

All are basically related to the same underlying issue which is that the ape function drop.tip(...,trim.internal=FALSE) will label the stem of a trimmed clade with the tip label "NA" when node labels are absent, but the node label itself when present. All of the internals used by reroot were originally built around the assumption that the trimmed edge would be labeled "NA". Unfortunately, this created issues for multiple internals - which I have been peeling back & fixing one-by-one.

This update is also in the latest (non-CRAN) phytools version which, as always, can be obtained from the phytools webpage.

Here is a quick demo of the bug & its fix.

First reroot behaving properly for a tree without node labels:

library(phytools)
## simulate a tree
set.seed(1)
tree<-rtree(n=26,tip.label=LETTERS)
plotTree(tree)
nodelabels()

plot of chunk unnamed-chunk-1

## here with a tree without node labels
t1<-reroot(tree,31,0.5*tree$edge.length[which(tree$edge[,2]==31)])
plotTree(t1)

plot of chunk unnamed-chunk-1

## we re-rooting in the middle of an edge so
## these should be equal
t1$edge.length[which(t1$edge[,1]==(Ntip(t1)+1))]
## [1] 0.3661569 0.3661569

OK, now here it is misbehaving. The tree is re-rooted in the right spot, but the edge lengths descended from the root are wrong:

## add node labels
tree$node.label<-paste("n",1:tree$Nnode+Ntip(tree),sep="")
plotTree(tree)
nodelabels(tree$node.label)

plot of chunk unnamed-chunk-2

t2<-reroot(tree,31,0.5*tree$edge.length[which(tree$edge[,2]==31)])
plotTree(t2) ## look at the edge lengths
nodelabels(t2$node.label)

plot of chunk unnamed-chunk-2

## these should be equal, but they are not
t2$edge.length[which(t1$edge[,1]==(Ntip(t2)+1))]
## [1] 2.2862548 0.3661569

OK, now let's load the fix and try again:

detach("package:phytools", unload=TRUE)
install.packages("phytools_0.4-49.tar.gz",type="source")
## Installing package into 'C:/Users/Liam/Documents/R/win-library/3.1'
## (as 'lib' is unspecified)
## inferring 'repos = NULL' from 'pkgs'
library(phytools)
t2<-reroot(tree,31,0.5*tree$edge.length[which(tree$edge[,2]==31)])
plotTree(t2)
nodelabels(t2$node.label)

plot of chunk unnamed-chunk-3

## these should be equal
t2$edge.length[which(t1$edge[,1]==(Ntip(t2)+1))]
## [1] 0.3661569 0.3661569

That's it.

Wednesday, February 11, 2015

Yet another bug fix for reroot for an input tree containing node labels

Recently, a phytools user correctly reported a bug in the phytools function reroot to re-root a tree along an edge when the input tree contained node labels. I thought I had identified it's cause and fixed it (details here). Specifically, the problem is caused by the fact that the function drop.tip(...,trim.internal=FALSE) will leave behind "NA" as a tip label for the trimmed clade when node labels are absent; however it will instead use the node labels if they are present. I fixed the (primarily internal) phytools function splitTree for this issue; however I failed to realize that I also needed to update the function drop.clade. I have now done this and posted the updated code here.

Let's try it with the old version & the new one:

library(phytools)
set.seed(1) ## set seed for reproducibility
tree<-rtree(n=26)
tree$tip.label<-LETTERS
plotTree(tree)

plot of chunk unnamed-chunk-1

node<-37
position=0.5*tree$edge.length[which(tree$edge[,2]==node)]
plotTree(reroot(tree,node,position)) ## works fine

plot of chunk unnamed-chunk-1

## now try with node labels
tree$node.label<-paste("n",1:tree$Nnode,sep="")
rerooted.tree<-reroot(tree,node,position) ## breaks
## Error in if (newroot == ROOT) {: argument is of length zero

Now let's update to the latest phytools version (not on CRAN) and re-attempt:

detach("package:phytools",unload=TRUE)
install.packages("phytools_0.4-48.tar.gz",type="source",repos=NULL)
## Installing package into 'C:/Users/Liam/Documents/R/win-library/3.1'
## (as 'lib' is unspecified)
library(phytools)
rerooted.tree<-reroot(tree,node,position)
plotTree(rerooted.tree)
nodelabels(rerooted.tree$node.label)

plot of chunk unnamed-chunk-2

The most recent working version of phytools can always be downloaded here.

That's it.

Tuesday, February 10, 2015

Bug fix for reroot for trees with node labels

A phytools user last night very helpfully reported that the phytools function reroot, which reroots the tree along an edge, throws an error if the input tree has node labels. Indeed, this seems to be the case. Here is the code supplied by the user (slightly modified) demonstrating the bug:

library(phytools)
set.seed(1) ## just for reproducibility
##  no node labels
tr<-rtree(12)
plotTree(tr)

plot of chunk unnamed-chunk-1

nod<-fastMRCA(tr,"t6","t10") ## arbitrarily
rtr<-reroot(tr,node.number=nod,position=0.5*tr$edge.length[which(tr$edge[,2]==nod)])
plotTree(rtr) ## works fine

plot of chunk unnamed-chunk-1

## now try the same with node labels
tr$node.label<-paste("n",1:tr$Nnode,sep="")
plotTree(tr)
nodelabels(tr$node.label)

plot of chunk unnamed-chunk-1

rtr<-reroot(tr, node.number=nod, position=0.5*tr$edge.length[which(tr$edge[,2]==nod)])
## Error in if (newroot == ROOT) {: argument is of length zero

This morning I was able to get to the bottom of this bug. Basically, the function reroot uses the phytools (mostly internal, but in the namespace) function splitTree to first split the tree at the desired point for re-rooting, then reroots the tree using the split point in the “parent” subtree as the new root. It identifies this split point in the parent tree because it has the tip label of "NA". Unfortunately, if the input tree has node labels, then the new node will be labeled not with "NA", but with the node label of the input tree.

I have now posted a new version of this function with this bug fixed. I accomplished this by checking for tips labeled first with "NA" and then, if none are found, with any of the node labels of the input tree. This will work fine unless any of the node labels are the same as tip labels (or "NA"). This probably should be avoided anyway.

Here is how it works:

detach("package:phytools",unload=TRUE)
## install new phytools with bug fix
install.packages("phytools_0.4-47.tar.gz",type="source",repos=NULL)
## Installing package into 'C:/Users/Liam/Documents/R/win-library/3.1'
## (as 'lib' is unspecified)
library(phytools)
packageVersion("phytools")
## [1] '0.4.47'
rtr<-reroot(tr, node.number=nod, position=0.5*tr$edge.length[which(tr$edge[,2]==nod)])
plotTree(rtr)
nodelabels(rtr$node.label)

plot of chunk unnamed-chunk-2

Great!

Monday, February 9, 2015

Bug fix in plotTree.singletons to permit multifurcating nodes

A phytools user correctly reports that the phytools function plotTree.singletons (which, as the function name suggests, plots trees containing singleton nodes) breaks if more than two edges emerge from a single node.

For instance:

library(phytools)
text<-"(A:0.1,(S:0.3)B:0.2,(C:0.3,D:0.4)E:0.5)F;"
tree<-read.tree(text=text)
tree ## doesn't work at all
## 
## Phylogenetic tree with 4 tips and 3 internal nodes.
## 
## Tip labels:
## [1] "A" "S" "B" "" 
## Node labels:
## [1] ""  "D" "C"
## 
## Rooted; includes branch lengths.
tree<-read.newick(text=text)
tree ## read in correctly
## 
## Phylogenetic tree with 4 tips and 3 internal nodes.
## 
## Tip labels:
## [1] "A" "S" "C" "D"
## Node labels:
## [1] "F" "B" "E"
## 
## Unrooted; includes branch lengths.
plotTree.singletons(tree)
## Error in xy.coords(x, y): 'x' and 'y' lengths differ

plot of chunk unnamed-chunk-1

Ok, the user very helpfully identified a fix in the code to address this issue. I also seem to have found a simpler fix. The code is available here and will also be in the next version of phytools.

Let's plot with the fixed code:

source("map.to.singleton.R")
plotTree.singletons(tree)

plot of chunk unnamed-chunk-2

OK, that's it.