X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=bayesian.cpp;h=6f65965bfb912e983d4de12dcbaa3dfe307e23e1;hp=eca63b1a9168613a92232777e520efcaa586f328;hb=df7e3ff9f68ef157b0328a2d353c3258c5d45d89;hpb=d6c0a11d1cecfac18b323285e7ffadb7f58e848f diff --git a/bayesian.cpp b/bayesian.cpp index eca63b1..6f65965 100644 --- a/bayesian.cpp +++ b/bayesian.cpp @@ -12,13 +12,14 @@ #include "phylosummary.h" #include "referencedb.h" /**************************************************************************************************/ -Bayesian::Bayesian(string tfile, string tempFile, string method, int ksize, int cutoff, int i, int tid, bool f) : +Bayesian::Bayesian(string tfile, string tempFile, string method, int ksize, int cutoff, int i, int tid, bool f, bool sh) : Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { try { ReferenceDB* rdb = ReferenceDB::getInstance(); threadID = tid; flip = f; + shortcuts = sh; string baseName = tempFile; if (baseName == "saved") { baseName = rdb->getSavedReference(); } @@ -27,7 +28,7 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { if (baseTName == "saved") { baseTName = rdb->getSavedTaxonomy(); } /************calculate the probablity that each word will be in a specific taxonomy*************/ - string tfileroot = baseTName.substr(0,baseTName.find_last_of(".")+1); + string tfileroot = m->getFullPathName(baseTName.substr(0,baseTName.find_last_of(".")+1)); string tempfileroot = m->getRootName(m->getSimpleName(baseName)); string phyloTreeName = tfileroot + "tree.train"; string phyloTreeSumName = tfileroot + "tree.sum"; @@ -63,7 +64,7 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { } saveIn.close(); } - + if(probFileTest && probFileTest2 && phyloTreeTest && probFileTest3 && FilesGood){ if (tempFile == "saved") { m->mothurOutEndLine(); m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory."); m->mothurOutEndLine(); } @@ -113,7 +114,7 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { WordPairDiffArr.resize(numKmers); for (int j = 0; j < wordGenusProb.size(); j++) { wordGenusProb[j].resize(genusNodes.size()); } - ofstream out; + ofstream out; ofstream out2; #ifdef USE_MPI @@ -124,23 +125,24 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { #endif - m->openOutputFile(probFileName, out); + if (shortcuts) { + m->openOutputFile(probFileName, out); - //output mothur version - out << "#" << m->getVersion() << endl; + //output mothur version + out << "#" << m->getVersion() << endl; - out << numKmers << endl; + out << numKmers << endl; - m->openOutputFile(probFileName2, out2); + m->openOutputFile(probFileName2, out2); - //output mothur version - out2 << "#" << m->getVersion() << endl; + //output mothur version + out2 << "#" << m->getVersion() << endl; + } #ifdef USE_MPI } #endif - //for each word for (int i = 0; i < numKmers; i++) { if (m->control_pressed) { break; } @@ -151,7 +153,7 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { if (pid == 0) { #endif - out << i << '\t'; + if (shortcuts) { out << i << '\t'; } #ifdef USE_MPI } @@ -159,12 +161,10 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { vector seqsWithWordi = database->getSequencesWithKmer(i); - map count; - for (int k = 0; k < genusNodes.size(); k++) { count[genusNodes[k]] = 0; } - //for each sequence with that word + vector count; count.resize(genusNodes.size(), 0); for (int j = 0; j < seqsWithWordi.size(); j++) { - int temp = phyloTree->getIndex(names[seqsWithWordi[j]]); + int temp = phyloTree->getGenusIndex(names[seqsWithWordi[j]]); count[temp]++; //increment count of seq in this genus who have this word } @@ -178,9 +178,9 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { //probabilityInThisTaxonomy = (# of seqs with that word in this taxonomy + probabilityInTemplate) / (total number of seqs in this taxonomy + 1); - wordGenusProb[i][k] = log((count[genusNodes[k]] + probabilityInTemplate) / (float) (genusTotals[k] + 1)); + wordGenusProb[i][k] = log((count[k] + probabilityInTemplate) / (float) (genusTotals[k] + 1)); - if (count[genusNodes[k]] != 0) { + if (count[k] != 0) { #ifdef USE_MPI int pid; MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are @@ -188,7 +188,7 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { if (pid == 0) { #endif - out << k << '\t' << wordGenusProb[i][k] << '\t' ; + if (shortcuts) { out << k << '\t' << wordGenusProb[i][k] << '\t' ; } #ifdef USE_MPI } @@ -204,8 +204,10 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { if (pid == 0) { #endif - out << endl; - out2 << probabilityInTemplate << '\t' << numNotZero << '\t' << log(probabilityInTemplate) << endl; + if (shortcuts) { + out << endl; + out2 << probabilityInTemplate << '\t' << numNotZero << '\t' << log(probabilityInTemplate) << endl; + } #ifdef USE_MPI } @@ -218,9 +220,10 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { if (pid == 0) { #endif - out.close(); - out2.close(); - + if (shortcuts) { + out.close(); + out2.close(); + } #ifdef USE_MPI } #endif @@ -230,7 +233,7 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { delete phyloTree; phyloTree = new PhyloTree(phyloTreeTest, phyloTreeName); - + //save probabilities if (rdb->save) { rdb->wordGenusProb = wordGenusProb; rdb->WordPairDiffArr = WordPairDiffArr; } } @@ -252,9 +255,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"); @@ -296,7 +298,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); @@ -306,7 +308,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; } @@ -372,6 +378,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 @@ -381,11 +388,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); } @@ -605,7 +614,7 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string istringstream iss (tempBuf,istringstream::in); float probTemp; iss >> zeroCountProb[i] >> numbers[i] >> probTemp; - WordPairDiffArr[i].prob = tempProb; + WordPairDiffArr[i].prob = probTemp; }