Wednesday, September 2, 2015

Bug fix in cophylo

A phytools user (Eliot Miller) reported some strange behavior in the function for co-phylogenetic plotting, cophylo, which can be replicated with the following example:

library(phytools)
tr1<-rtree(n=26,tip.label=LETTERS)
tr2<-rtree(n=26,tip.label=LETTERS)
obj<-cophylo(tr1,tr2)
## Rotating nodes to optimize matching...
## Done.
plot(obj) ## works

plot of chunk unnamed-chunk-1

## read & write to NEXUS tree file:
write.nexus(c(tr1,tr2),file="nexus.trees")
trees<-read.nexus(file="nexus.trees")
obj<-cophylo(trees[[1]],trees[[2]])
## Rotating nodes to optimize matching...
## Done.
plot(obj) ## fails

plot of chunk unnamed-chunk-1

This is because the NEXUS trees use a TRANSLATE tabel and consequently the order of the tip labels is given by a single vector for all trees, and is thus not consistenly 'cladewise' across trees (even if the rows of edge are themselves cladewise.

This can be seen using:

str(trees)
## Class "multiPhylo"
## List of 2
##  $ UNTITLED:List of 3
##   ..$ edge       : int [1:50, 1:2] 27 28 29 29 30 31 31 30 28 32 ...
##   ..$ edge.length: num [1:50] 0.0507 0.6907 0.4726 0.7195 0.6709 ...
##   ..$ Nnode      : int 25
##   ..- attr(*, "class")= chr "phylo"
##   ..- attr(*, "order")= chr "cladewise"
##  $ UNTITLED:List of 3
##   ..$ edge       : int [1:50, 1:2] 27 27 28 29 30 30 29 31 31 32 ...
##   ..$ edge.length: num [1:50] 0.0255 0.5955 0.0601 0.5403 0.294 ...
##   ..$ Nnode      : int 25
##   ..- attr(*, "class")= chr "phylo"
##   ..- attr(*, "order")= chr "cladewise"
##  - attr(*, "TipLabel")= chr [1:26] "T" "Q" "Y" "R" ...

We can fix this using the (classic), if a bit inelegant, read.tree/ write.tree hack - now added to the cophylo function (details here).

source("../phytools/R/cophylo.R")
trees<-read.nexus(file="nexus.trees")
obj<-cophylo(trees[[1]],trees[[2]])
## Rotating nodes to optimize matching...
## Done.
plot(obj) ## works?

plot of chunk unnamed-chunk-3

This fixed version can be installed from GitHub using devtools.

That's it for now.

No comments:

Post a Comment