X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=bayesian.cpp;h=8278afb32e1d2c028a83ffdece17ed9910a78e20;hp=49be4af57ff66f46ac04912535d7918518cc75d2;hb=d1c97b8c04bb75faca1e76ffad60b37a4d789d3d;hpb=90708fe9701e3827e477c82fb3652539c3bf2a0d diff --git a/bayesian.cpp b/bayesian.cpp index 49be4af..8278afb 100644 --- a/bayesian.cpp +++ b/bayesian.cpp @@ -145,6 +145,8 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { //for each word for (int i = 0; i < numKmers; i++) { + //m->mothurOut("[DEBUG]: kmer = " + toString(i) + "\n"); + if (m->control_pressed) { break; } #ifdef USE_MPI @@ -239,7 +241,9 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { } } + if (m->debug) { m->mothurOut("[DEBUG]: about to generateWordPairDiffArr\n"); } generateWordPairDiffArr(); + if (m->debug) { m->mothurOut("[DEBUG]: done generateWordPairDiffArr\n"); } //save probabilities if (rdb->save) { rdb->wordGenusProb = wordGenusProb; rdb->WordPairDiffArr = WordPairDiffArr; } @@ -255,9 +259,8 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { /**************************************************************************************************/ Bayesian::~Bayesian() { try { - - delete phyloTree; - if (database != NULL) { delete database; } + if (phyloTree != NULL) { delete phyloTree; } + if (database != NULL) { delete database; } } catch(exception& e) { m->errorOut(e, "Bayesian", "~Bayesian"); @@ -299,7 +302,7 @@ string Bayesian::getTaxonomy(Sequence* seq) { } } - if (queryKmers.size() == 0) { m->mothurOut(seq->getName() + "is bad."); m->mothurOutEndLine(); simpleTax = "unknown;"; return "unknown;"; } + if (queryKmers.size() == 0) { m->mothurOut(seq->getName() + " is bad. It has no kmers of length " + toString(kmerSize) + "."); m->mothurOutEndLine(); simpleTax = "unknown;"; return "unknown;"; } int index = getMostProbableTaxonomy(queryKmers); @@ -309,7 +312,11 @@ string Bayesian::getTaxonomy(Sequence* seq) { //bootstrap - to set confidenceScore int numToSelect = queryKmers.size() / 8; + if (m->debug) { m->mothurOut(seq->getName() + "\t"); } + tax = bootstrapResults(queryKmers, index, numToSelect); + + if (m->debug) { m->mothurOut("\n"); } return tax; } @@ -375,6 +382,7 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { int seqTaxIndex = tax; TaxNode seqTax = phyloTree->get(tax); + while (seqTax.level != 0) { //while you are not at the root itBoot2 = confidenceScores.find(seqTaxIndex); //is this a classification we already have a count on @@ -384,11 +392,13 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { confidence = itBoot2->second; } + if (m->debug) { m->mothurOut(seqTax.name + "(" + toString(((confidence/(float)iters) * 100)) + ");"); } + if (((confidence/(float)iters) * 100) >= confidenceThreshold) { confidenceTax = seqTax.name + "(" + toString(((confidence/(float)iters) * 100)) + ");" + confidenceTax; simpleTax = seqTax.name + ";" + simpleTax; } - + seqTaxIndex = seqTax.parent; seqTax = phyloTree->get(seqTax.parent); }