Repository navigation
Different results assignTaxonomy #538
Description
Howdy,
R version 3.5.0 (2018-04-23)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.5
dada2 * 1.8.0 2018-07-31
I have run the following command four times...
tax_silva <- assignTaxonomy(seqtab, "silva_nr_v132_train_set.fa.gz", multithread = TRUE)
and get slightly different results each time. For example,
ASV12456; Bacteria; Proteobacteria; Gammaproteobacteria; Vibrionales; Vibrionaceae; NA
vs.
ASV12456; Bacteria; Proteobacteria; Gammaproteobacteria; NA; NA; NA
Is this to be expected?
I have not been able to identify and rhyme or reason to the patterns. I am trying now to set minBoot manually and call the outputBootstraps flag to have a closer look. But since this is taking ~3hrs a pop, I wanted to ask first.
If you have a quick answer, great, otherwise I am happy to send whatever files you need to reproduce this result. Just let me now what you need and what format would work best:)
Thanks!
Activity
jarrodscott commentedon Aug 8, 2018
I should also add that I ran my whole pipeline three times to make sure there were no underlying differences in the ASVs (e.g., length). The sequence sets I got back were identical--this included the
removeBimeraDenovostep. I am pretty convinced this is related toassignTaxonomyThanks!
benjjneb commentedon Aug 8, 2018
Yes, completely expected.
assignTaxonomyis implementing the naive Bayesian classifier method of Wang et al, 2007. That method uses bootstrapping over subsamples of the kmer-shredded read to assign confidence to each level of classification. That means there is some randomness involved, as the subsamples are randomly chosen and different each time.I believe
assignTaxonomyrespects R's random number generator, so for complete reproducibility you can e.g.set.seed(100)just before runningassignTaxonomy, and it should give the exact same output each time.jarrodscott commentedon Aug 8, 2018
Perfect!!!!!!! Thank you so much for clarifying and for the suggestion. I was just looking through that paper. Sometimes ignorance gets the best of me :) Have a nice day.
soluna1 commentedon May 3, 2019
Hi,
I think that it is better to post here (although this issue is already closed), as my question is similar. I agree that there are some randomness in the function as it is based on a probabilistic model, but I'm trying to understand why the same sequence is assigned to different genera when using different databases. The data are the following:
asv sequences: 130 bp of the v9 18s
database 1: pr2_version_4.11.1_dada2.fasta
database 2: custom-made database with only the v9 section of 18s
Both databases contain a very similar list of taxa, and in both the "problematic-genera" is included.
The code use to run the function "AssignTaxonomy" is the same, set.seed(42) in both cases, and the same parameters :
assignTaxonomy(seqtab, tax_db, taxLevels = c("Kingdom","Supergroup","Division","Class","Order","Family","Genus","Species"), multithread=T, tryRC=T, minBoot=80, outputBootstraps=T)
The problem:
The following asv_sequence "TCCTACCGATTGGATTGGCTGGTGAAGCGTTCGGATTGTGGGTTCGGAGCGGTTCGCCGTTCTGGGCTTGCGAGAAGTTCGTTAAACCCTCCTATCTAGAGGAAGGAGAAGTCGTAACAAGGTTTC"
is assigned to different genera in each case:
for pr2 --> Pseudoscourfieldia
for custom db--> Pycnococcus
in both cases the Bootstrap output is 100%.
This is not a big issue, as the family is the same for both genera, so I can go a level before. But my question is how I'll know when the 100% bootstrapping is really "true" or at least acceptable?
In this case, I would have expected that the 100% boostraping were retrieved at family level, not a genus level. Do you agree with this?
thanks in advance for your collaboration.
Pseudoscourfieldia
Pseudoscourfieldia
benjjneb commentedon May 5, 2019
@soluna1 Could you repost this as a new issue? I think it deserves its own thread.
soluna1 commentedon May 5, 2019
ok, I'll do.