search

xavierdidelot / BactDating

Bayesian inference of ancestral dates on bacterial phylogenetic trees
https://xavierdidelot.github.io/BactDating
MIT License
80 stars 15 forks source link

Title : Error with bactdate : Error in if [...] missing value where TRUE/FALSE needed #52

Closed fredericQC closed 2 years ago

fredericQC commented 2 years ago

Hi!

I used BactDating with some data and it worked well, but I'm having an issue with the tree I'm currently working on.

Here's some info about how I built the tree:

I wanted a tree based on core SNP, without recombinant regions. I have a reference genome and about 70 samples.

Workflow:

tree from FastTree:

(((fr_18:0.000029229,((fr_13:0.0,fr_5:0.0):0.000000005,(fr_6:0.000014614,((Abau_740:0.000000005,(Abau_584:0.000131537,Abau_774:0.000000005)0.915:0.000014614)0.981:0.000073074,((Abau_460:0.000014610,Abau_495:0.000058462)0.986:0.000087696,((PID-0586-AB850-BEN:0.000233878,(((MRSN505383:0.0,MRSN505505:0.0,MRSN505532:0.0):0.000263129,(A10:0.000087690,((((778944:0.0,Abau_501:0.0):0.000014614,(Abau_500:0.0,Abau_502:0.0):0.000000005)0.999:0.000175400,(UK_5:0.000029227,(UK_4:0.000014614,UK_7:0.000000005)0.951:0.000000005)0.995:0.000131550)0.990:0.000116927,(A30:0.000087690,(AMA0341:0.000160783,EGY-R14:0.000292363)0.319:0.000000005)0.999:0.000131543)0.286:0.000000005)0.915:0.000043846)0.976:0.000073072,(1560430:0.000014614,AB32-VUB:0.000014614)1.000:0.000394760)0.278:0.000000006)0.956:0.000043844,((164_Tony:0.000409362,(ACN1-Sample12:0.000131553,Acinetobacter_baumannii_AA-014:0.000102319)0.741:0.000014595)0.888:0.000029237,((MRSN20070:0.000000005,(Ab-Pak-Pesh-01:0.000029228,(Ab-Pak-Pesh-28:0.000058457,Ab-Pak-Pesh-07:0.000131542)0.528:0.000000005)0.958:0.000058458)1.000:0.000175399,((IC8:0.238867666,(((IC2:0.226190705,IC5:0.223374578)1.000:0.025365511,(IC6:0.237299820,(IC4:0.229645992,IC7:0.234270730)1.000:0.019935361)1.000:0.010857958)1.000:0.009943582,(IC9:0.241696712,(IC1:0.188770933,IC3:0.200953945)1.000:0.043849008)1.000:0.016512374)1.000:0.021690933)1.000:0.230468515,(MRSN15151:0.000029229,(MRSN22111:0.000204644,(MRSN636940:0.000014615,MRSN636958:0.000029228)0.979:0.000073077)0.844:0.000000006)0.974:0.000158645)0.972:0.000162862)0.417:0.000000005)0.754:0.000014613)0.940:0.000043848)0.924:0.000043840)0.757:0.000014614)0.875:0.000000005)0.881:0.000000005)0.892:0.000014614,((Reference:0.0,fr_14:0.0,fr_15:0.0,fr_17:0.0,fr_24:0.0,fr_27:0.0,fr_9:0.0):0.000000005,((((fr_12:0.000014614,fr_4:0.000029228)0.789:0.000000005,(fr_1:0.000014614,fr_2:0.000000005)0.861:0.000014614)0.000:0.000000005,fr_16:0.000014614)0.882:0.000014614,fr_25:0.000029229)0.894:0.000000005)0.916:0.000000005)0.853:0.000014614,fr_8:0.000014614,(((fr_19:0.000029228,fr_21:0.000014614)0.890:0.000000005,(fr_20:0.000029228,fr_23:0.000029228)0.902:0.000000005)0.000:0.000000005,((Acinetobacter_baumannii_isolate_K50:0.000000005,fr_10:0.000014614)0.865:0.000014614,(fr_11:0.0,fr_22:0.0,fr_26:0.0,fr_3:0.0,fr_7:0.0):0.000000005)0.000:0.000000005)0.915:0.000000005);

Then I ran the following command in R to get the rooted, outgroup-removed, dated tree with BactDating.:

library(BactDating)
library(ape)

# the following tree (t) is the one written above
t=read.tree('ST158/phylogeny/all_IC_outgroup.clean.core.newick')
t$edge.length=t$edge.length*68692

outgroups=c("IC1","IC2","IC3","IC4","IC5","IC6","IC7","IC8","IC9")

rt=root(t,outgroups,resolve.root=TRUE)

rm_rt=drop.tip(rt,outgroups)

d=c(2008.0,2008.0,2008.0,2008.0,2008.0,2012.0,2008.0,2011.0,2011.0,2008.0,2016.0,2016.0,2016.0,2015.0,NA,2012.0,2012.0,2012.0,2011.0,2011.0,2011.0,2015.0,2012.0,NA,NA,2014.0,NA,NA,2008.0,2013.0,2013.0,2013.0,2013.0,2012.0,2013.0,2016.0,2016.0,2007.0,2008.0,2007.0,2007.0,2008.0,2007.0,2007.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2008.0,2007.0,2008.0,2008.0)

res=bactdate(rm_rt,d)

But I get the following when I try to run bactdate:

> res=bactdate(rm_rt,d)
Error in if (log(runif(1)) < l2 - l + dgamma(mu2, shape = 0.001, scale = 1000,  : 
  missing value where TRUE/FALSE needed

Do you know what's happening here?

Thanks!

xavierdidelot commented 2 years ago

Hello Frederic,

This looks like a nice dataset on which BactDating should work well. Your problem is that you have multifurcations in the tree which is not allowed. So before you run the bactdate command you would need to replace these with binary nodes, using: rm_rt=multi2di(rm_rt)

Best wishes, Xavier

fredericQC commented 2 years ago

It solved the problem! Thanks!

appliedmicrobiologyresearch commented 1 month ago

I am also getting this error message. My data also has many multifurcations, but this solution doesn't seem to help. But the roottotip is giving a nice signal, so I'd love to get BactDating running. Can you help and what information would you need? Thanks! Helena

xavierdidelot commented 1 month ago

Hi Helena,

Do you mean you get this error message even after using multi2di? If so there must be something else in your tree that is not acceptable for BactDating, maybe a branch of negative length for example?

Best wishes, Xavier

appliedmicrobiologyresearch commented 4 weeks ago

Thanks Xavier, this is the tree: (((((ZIB9084:10.215211,(ZIB8559:7.418663,(ZIB9080:4.698386,(ZIB8551:7.57946,ZIB8536:8.13356):3.171639):4.497993):0.000123):1.77989,(((ZIB8528:5.495356,(ZIB8545:1.752863,(((ZIB8565:11.310996,ZIB8570:0.000123):3.540978,((ZIB8572:10.561228,ZIB9076:1.234011):6.216047,ZIB8525:0.000123):1.755952):0.000123,(((ZIB8538:5.47764,(ZIB8544:5.200626,((ZIB8511:0.000123,ZIB8502:3.457654):1.730641,(((ZIB8546:1.863315,((ZIB9095:1.80742,ZIB8560:1.725934):1.714019,((ZIB9059:5.331406,ZIB8542:0.000123):3.531173,((ZIB8562:3.512913,ZIB8569:0.000123):1.662098,((ZIB9065:13.389436,ZIB9056:0.000123):3.515738,ZIB9055:0.000123):3.596584):1.739652):0.000123):1.778895):1.806307,((ZIB8524:3.383465,(ZIB8520:6.91435,(((((ZIB8508:8.661039,(((((ZIB8552:8.438435,ZIB8556:0.000123):1.605677,ZIB8555:0.000123):3.398699,ZIB9068:1.696822):3.429344,(ZIB9090:0.000123,ZIB9094:3.882909):0.000123):3.321362,ZIB9085:1.653697):0.000123):3.436398,(ZIB8516:5.347928,((ZIB9049:7.305974,((ZIB8558:7.418573,((ZIB9089:3.750213,ZIB8547:1.96391):1.754389,(ZIB8522:3.622208,((ZIB8537:12.602034,ZIB9067:3.449689):1.988168,(ZIB8514:7.292768,ZIB8526:0.000123):1.735797):1.831931):0.000123):0.000123):1.78268,ZIB8554:0.000123):0.000123):1.827639,ZIB8500:5.644176):1.694359):0.000123):1.725867,ZIB8509:0.000123):1.680008,((ZIB8513:5.181476,(((ZIB9070:0.000123,ZIB9088:5.300639):3.33318,ZIB8507:0.000123):1.644829,((ZIB8550:3.604687,ZIB8517:1.874589):3.538267,ZIB9074:5.173912):0.000123):0.000123):1.669057,ZIB8510:0.000123):0.000123):1.679951,ZIB8534:0.000123):1.680874):0.000123):1.725032,((ZIB8539:7.026978,ZIB8531:0.000123):3.478508,ZIB8519:0.000123):0.000123):1.783252):1.766322,ZIB8518:0.000123):1.738299):0.000123):1.652032):1.881334,(((ZIB9077:1.802125,ZIB8566:7.280289):1.772585,(((ZIB8527:0.000123,ZIB8530:9.780227):1.97138,ZIB8532:3.796662):1.73164,(((ZIB8577:6.08061,ZIB8563:0.000123):7.320205,ZIB8515:0.000123):1.807267,ZIB9086:3.649018):0.000123):3.684854):1.809135,ZIB8561:0.000123):0.000123):0.000123,((ZIB8529:3.834555,ZIB8573:3.807819):1.665055,((ZIB8568:3.575752,ZIB9069:1.741389):1.739525,ZIB8567:0.000123):3.677617):1.733747):7.419116):0.000123):1.911579):1.568109,ZIB8503:1.817627):3.595802,ZIB8523:3.531688):0.000123):3.577473,ZIB8505:3.416388):0.000123,(ZIB8506:1.858531,ZIB8501:5.259804):1.584795):2.376368,((ZIB8543:7.966907,ZIB9079:1.676918):1.756774,(ZIB8553:3.534481,(ZIB8548:0.000123,(((ZIB9048:1.807364,(ZIB9075:3.954625,ZIB9050:5.603428):1.489329):1.693696,ZIB9057:5.281991):0.000123,((ZIB8535:9.381135,ZIB8540:10.959709):0.000123,(ZIB9092:0.000123,(((ZIB8557:9.230419,ZIB9005:3.392617):2.158717,(ZIB9022:0.000123,ZIB9030:0.000123):0.000123):1.716845,(ZIB9027:3.584503,((ZIB9029:3.577916,(ZIB9007:6.035613,ZIB8512:0.000123):3.490993):3.610029,(ZIB9006:0.000123,((((ZIB9032:3.390037,(ZIB9018:0.000123,(ZIB9012:9.60873,ZIB9021:3.511783):7.583807):1.751557):1.566972,ZIB9023:3.450344):1.802503,(ZIB9043:5.273805,((ZIB9004:3.826974,ZIB9016:3.49855):4.125312,(((((ZIB9054:9.789342,(ZIB9072:6.732906,(ZIB9064:0.000123,ZIB9053:10.930529):5.560166):5.042449):0.000123,ZIB9051:2.22508):2.878528,ZIB9010:6.012818):1.901306,(ZIB9001:6.30488,(ZIB9009:6.160383,(ZIB9008:13.41917,ZIB9003:2.080081):0.000123):2.03167):8.088038):2.394248,ZIB9078:17.420065):0.000123):3.636723):0.000123):0.000123,(ZIB9046:4.037602,ZIB9063:3.442571):1.729443):6.888832):1.688743):1.681314):1.677395):3.543775):1.635865):5.31918):1.704982):1.667829):2.002505):1.237334):0.0;

xavierdidelot commented 4 weeks ago

It seems to work fine with this tree on my side. Could you try:

library(ape)
library(BactDating)
t=read.tree('tree.nwk')
d=runif(Ntip(t),2010,2020)
bactdate(t,d,showProgress = T)

If this works with made up dates on line 4 but not with the actual dates then maybe there is something wrong with your dates.

appliedmicrobiologyresearch commented 2 weeks ago

Aha! Cunning. Yes, there WAS something wrong with my dates. I'm now getting nice consistent results :-) THANKS!