xavierdidelot / TransPhylo

Reconstruction of transmission trees using genomic data
http://xavierdidelot.github.io/TransPhylo/
GNU General Public License v2.0
60 stars 22 forks source link

Error "negative branch lengths" w/o negative branch lengths #7

Closed phiweger closed 6 years ago

phiweger commented 6 years ago

Hi,

This is my input tree in newick format:

(VA13414_h:0.15657,(BK17487_h:0.09975,(VA2520_h:0.97287,VA15976_h:0.38081,(VA21335_h:0.19916,(TP1148_h:0.38815,VA24749_h:0.09725):1.29904):0.35961,(VA22695_h:0.48789,((VA25100_h:0.00000,VA27811_h:0.10130):0.48282,((UR8729_h:0.51740,(TP3995_h:0.28445,VA18593_h:0.26528):0.36710):0.52579,(VA22360_h:0.23612,(VA22374_h:0.00000,(VA6218_h:0.14273,VA8292_h:0.20844):0.35076):0.23612):0.00967):0.14668):0.08447):0.12016):0.31647):0.43191):0.10000;

Then I run

library('TransPhylo')
library(ape)

ptree <- ptreeFromPhylo(read.tree('tree.nwk'), dateLastSample=2013.3)
record <- inferTTree(ptree, dateT=2014)

Error in if (ptree$ptree[ptree$ptree[i, j], 1] - ptree$ptree[i, 1] < 0) stop("The phylogenetic tree contains negative branch lengths!") :
  argument is of length zero

p<-phyloFromPTree(ptree)

Error in tr$edge.length[iedge] <- ptree[ptree[i, 2], 1] - ptree[i, 1] :
  replacement has length zero

However, the tree does not contain negative branch lengths and the format seems to be right compared to the simulated data in the tutorial.

I tried the output from FastTree, RAXML and Treetime, all w/ the same error. What could be the cause?

Thank you!

xavierdidelot commented 6 years ago

Hello,

Indeed your tree does not contain any branch of negative length and so the error message is a bit misleading. The problem though is that there are two branches of length zero, and also at least one multifurcation. You could fix these issues using the following code:

t <- read.tree('tree.nwk')
t <- multi2di(t) # remove multifurcations
t$edge.length <- pmax(t$edge.length,1/365) # use a day as minimum branch length
ptree <- ptreeFromPhylo(t, dateLastSample=2013.3)
record <- inferTTree(ptree, dateT=2014)

Note however that the tree should be a dated tree, so the output of FastTree of raxml is not acceptable, but the output of treetime should be fine.

phiweger commented 6 years ago

Thank you very much for your (ultra-quick) response -- it did work and is now happily fitting the model.