X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bayesian.cpp;h=00467ff0feb0a5df82314685520943376582278e;hb=905cc2b0bd18c5ce611b048d785e93859865a5ea;hp=5272b2d1761285f934e9e8e1a474fd27f1edc328;hpb=a5afca18544555fba2d9c3670ad1f8574916b0a0;p=mothur.git diff --git a/bayesian.cpp b/bayesian.cpp index 5272b2d..00467ff 100644 --- a/bayesian.cpp +++ b/bayesian.cpp @@ -9,56 +9,100 @@ #include "bayesian.h" #include "kmer.hpp" +#include "phylosummary.h" /**************************************************************************************************/ -Bayesian::Bayesian(string tfile, string tempFile, string method, int ksize, int cutoff, bool p) : -Classify(tfile, tempFile, method, ksize, 0.0, 0.0, 0.0, 0.0), kmerSize(ksize), confidenceThreshold(cutoff), probs(p) { +Bayesian::Bayesian(string tfile, string tempFile, string method, int ksize, int cutoff, int i) : +Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { try { - numKmers = database->getMaxKmer() + 1; - - //initialze probabilities - wordGenusProb.resize(numKmers); - - genusNodes = phyloTree->getGenusNodes(); - - for (int j = 0; j < wordGenusProb.size(); j++) { wordGenusProb[j].resize(genusNodes.size()); } - - //reset counts because we are on a new word - for (int j = 0; j < genusNodes.size(); j++) { - TaxNode temp = phyloTree->get(genusNodes[j]); - genusTotals.push_back(temp.accessions.size()); - } - - /************calculate the probablity that each word will be in a specific taxonomy*************/ - ofstream out; - string probFileName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer.prob"; - ifstream probFileTest(probFileName.c_str()); + string phyloTreeName = tfile.substr(0,tfile.find_last_of(".")+1) + "tree.train"; + string probFileName = tfile.substr(0,tfile.find_last_of(".")+1) + tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer.prob"; + string probFileName2 = tfile.substr(0,tfile.find_last_of(".")+1) + tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer.numNonZero"; + ofstream out; ofstream out2; - string probFileName2 = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer.numNonZero"; + + ifstream phyloTreeTest(phyloTreeName.c_str()); ifstream probFileTest2(probFileName2.c_str()); + ifstream probFileTest(probFileName.c_str()); int start = time(NULL); - if(probFileTest && probFileTest2){ - mothurOut("Reading template probabilities... "); cout.flush(); - readProbFile(probFileTest, probFileTest2); + if(probFileTest && probFileTest2 && phyloTreeTest){ + m->mothurOut("Reading template taxonomy... "); cout.flush(); + + phyloTree = new PhyloTree(phyloTreeTest, phyloTreeName); + + m->mothurOut("DONE."); m->mothurOutEndLine(); + + genusNodes = phyloTree->getGenusNodes(); + genusTotals = phyloTree->getGenusTotals(); + + m->mothurOut("Reading template probabilities... "); cout.flush(); + readProbFile(probFileTest, probFileTest2, probFileName, probFileName2); + }else{ - mothurOut("Calculating template probabilities... "); cout.flush(); + + //create search database and names vector + generateDatabaseAndNames(tfile, tempFile, method, ksize, 0.0, 0.0, 0.0, 0.0); + + genusNodes = phyloTree->getGenusNodes(); + genusTotals = phyloTree->getGenusTotals(); + + m->mothurOut("Calculating template taxonomy tree... "); cout.flush(); + + phyloTree->printTreeNodes(phyloTreeName); + + m->mothurOut("DONE."); m->mothurOutEndLine(); + + m->mothurOut("Calculating template probabilities... "); cout.flush(); + + numKmers = database->getMaxKmer() + 1; + + //initialze probabilities + wordGenusProb.resize(numKmers); + + for (int j = 0; j < wordGenusProb.size(); j++) { wordGenusProb[j].resize(genusNodes.size()); } + + + #ifdef USE_MPI + int pid; + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + + if (pid == 0) { + #endif ofstream out; openOutputFile(probFileName, out); + out << numKmers << endl; + ofstream out2; openOutputFile(probFileName2, out2); + #ifdef USE_MPI + } + #endif + + //for each word for (int i = 0; i < numKmers; i++) { + if (m->control_pressed) { break; } + #ifdef USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + + if (pid == 0) { + #endif + out << i << '\t'; + #ifdef USE_MPI + } + #endif + vector seqsWithWordi = database->getSequencesWithKmer(i); map count; @@ -77,22 +121,59 @@ Classify(tfile, tempFile, method, ksize, 0.0, 0.0, 0.0, 0.0), kmerSize(ksize), c for (int k = 0; k < genusNodes.size(); k++) { //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)); - if (count[genusNodes[k]] != 0) { out << k << '\t' << wordGenusProb[i][k] << '\t'; numNotZero++; } + if (count[genusNodes[k]] != 0) { + #ifdef USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + if (pid == 0) { + #endif + + out << k << '\t' << wordGenusProb[i][k] << '\t'; + + #ifdef USE_MPI + } + #endif + + numNotZero++; + } } + + #ifdef USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + if (pid == 0) { + #endif + out << endl; out2 << probabilityInTemplate << '\t' << numNotZero << endl; + + #ifdef USE_MPI + } + #endif } + #ifdef USE_MPI + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + if (pid == 0) { + #endif + out.close(); out2.close(); + + #ifdef USE_MPI + } + #endif + + //read in new phylotree with less info. - its faster + ifstream phyloTreeTest(phyloTreeName.c_str()); + delete phyloTree; + + phyloTree = new PhyloTree(phyloTreeTest, phyloTreeName); } - - - mothurOut("DONE."); mothurOutEndLine(); - mothurOut("It took " + toString(time(NULL) - start) + " seconds get probabilities. "); mothurOutEndLine(); + + m->mothurOut("DONE."); m->mothurOutEndLine(); + m->mothurOut("It took " + toString(time(NULL) - start) + " seconds get probabilities. "); m->mothurOutEndLine(); } catch(exception& e) { - errorOut(e, "Bayesian", "getTaxonomy"); + m->errorOut(e, "Bayesian", "Bayesian"); exit(1); } } @@ -105,33 +186,30 @@ string Bayesian::getTaxonomy(Sequence* seq) { //get words contained in query //getKmerString returns a string where the index in the string is hte kmer number //and the character at that index can be converted to be the number of times that kmer was seen + string queryKmerString = kmer.getKmerString(seq->getUnaligned()); + vector queryKmers; for (int i = 0; i < queryKmerString.length(); i++) { if (queryKmerString[i] != '!') { //this kmer is in the query queryKmers.push_back(i); } } - + + if (queryKmers.size() == 0) { m->mothurOut(seq->getName() + "is bad."); m->mothurOutEndLine(); return "bad seq"; } + int index = getMostProbableTaxonomy(queryKmers); + + if (m->control_pressed) { return tax; } //bootstrap - to set confidenceScore - if (probs) { - int numToSelect = queryKmers.size() / 8; - tax = bootstrapResults(queryKmers, index, numToSelect); - }else{ - TaxNode seqTax = phyloTree->get(index); - while (seqTax.level != 0) { //while you are not at the root - tax = seqTax.name + ";" + tax; - seqTax = phyloTree->get(seqTax.parent); - } - simpleTax = tax; - } - + int numToSelect = queryKmers.size() / 8; + tax = bootstrapResults(queryKmers, index, numToSelect); + return tax; } catch(exception& e) { - errorOut(e, "Bayesian", "getTaxonomy"); + m->errorOut(e, "Bayesian", "getTaxonomy"); exit(1); } } @@ -151,7 +229,9 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { map::iterator itBoot2; map::iterator itConvert; - for (int i = 0; i < 100; i++) { + for (int i = 0; i < iters; i++) { + if (m->control_pressed) { return "control"; } + vector temp; for (int j = 0; j < numToSelect; j++) { @@ -175,9 +255,10 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { }else{ confidenceScores[taxonomy.level][taxonomy.name]++; } - + taxonomy = phyloTree->get(taxonomy.parent); } + } string confidenceTax = ""; @@ -185,7 +266,7 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { TaxNode seqTax = phyloTree->get(tax); while (seqTax.level != 0) { //while you are not at the root - + itBoot2 = confidenceScores[seqTax.level].find(seqTax.name); //is this a classification we already have a count on int confidence = 0; @@ -193,31 +274,31 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { confidence = confidenceScores[seqTax.level][seqTax.name]; } - if (confidence >= confidenceThreshold) { - confidenceTax = seqTax.name + "(" + toString(confidence) + ");" + confidenceTax; + if (((confidence/(float)iters) * 100) >= confidenceThreshold) { + confidenceTax = seqTax.name + "(" + toString(((confidence/(float)iters) * 100)) + ");" + confidenceTax; simpleTax = seqTax.name + ";" + simpleTax; } seqTax = phyloTree->get(seqTax.parent); } + if (confidenceTax == "") { confidenceTax = "unclassified;"; simpleTax = "unclassified;"; } return confidenceTax; } catch(exception& e) { - errorOut(e, "Bayesian", "bootstrapResults"); + m->errorOut(e, "Bayesian", "bootstrapResults"); exit(1); } } /**************************************************************************************************/ int Bayesian::getMostProbableTaxonomy(vector queryKmer) { try { - int indexofGenus; + int indexofGenus = 0; double maxProbability = -1000000.0; //find taxonomy with highest probability that this sequence is from it for (int k = 0; k < genusNodes.size(); k++) { - //for each taxonomy calc its probability double prob = 1.0; for (int i = 0; i < queryKmer.size(); i++) { @@ -234,7 +315,7 @@ int Bayesian::getMostProbableTaxonomy(vector queryKmer) { return indexofGenus; } catch(exception& e) { - errorOut(e, "Bayesian", "getMostProbableTaxonomy"); + m->errorOut(e, "Bayesian", "getMostProbableTaxonomy"); exit(1); } } @@ -261,46 +342,170 @@ map Bayesian::parseTaxMap(string newTax) { } catch(exception& e) { - errorOut(e, "Bayesian", "parseTax"); + m->errorOut(e, "Bayesian", "parseTax"); exit(1); } } /**************************************************************************************************/ -void Bayesian::readProbFile(ifstream& in, ifstream& inNum) { +void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string inNumName) { try{ - int kmer, name, count; count = 0; - vector num; num.resize(numKmers); - float prob; - vector zeroCountProb; zeroCountProb.resize(numKmers); - - while (inNum) { - inNum >> zeroCountProb[count] >> num[count]; - count++; - gobble(inNum); - } - inNum.close(); + #ifdef USE_MPI + + int pid, num, num2; + vector positions; + vector positions2; + + MPI_Status status; + MPI_File inMPI; + MPI_File inMPI2; + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + + char inFileName[1024]; + strcpy(inFileName, inNumName.c_str()); + + char inFileName2[1024]; + strcpy(inFileName2, inName.c_str()); + + MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI); //comm, filename, mode, info, filepointer + MPI_File_open(MPI_COMM_WORLD, inFileName2, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI2); //comm, filename, mode, info, filepointer + + if (pid == 0) { + positions = setFilePosEachLine(inNumName, num); + + //send file positions to all processes + MPI_Bcast(&num, 1, MPI_INT, 0, MPI_COMM_WORLD); //send numSeqs + MPI_Bcast(&positions[0], (num+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos + + positions2 = setFilePosEachLine(inName, num2); + + //send file positions to all processes + MPI_Bcast(&num2, 1, MPI_INT, 0, MPI_COMM_WORLD); //send numSeqs + MPI_Bcast(&positions2[0], (num2+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos + + }else{ + MPI_Bcast(&num, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs + positions.resize(num); + MPI_Bcast(&positions[0], (num+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions + + MPI_Bcast(&num2, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs + positions2.resize(num2); + MPI_Bcast(&positions2[0], (num2+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions + + } - while(in) { - in >> kmer; + //read numKmers + int length = positions2[1] - positions2[0]; + char* buf = new char[length]; + + MPI_File_read_at(inMPI2, positions2[0], buf, length, MPI_CHAR, &status); + + string tempBuf = buf; + if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); } + delete buf; + + istringstream iss (tempBuf,istringstream::in); + iss >> numKmers; + + //initialze probabilities + wordGenusProb.resize(numKmers); - //set them all to zero value - for (int i = 0; i < genusNodes.size(); i++) { - wordGenusProb[kmer][i] = log(zeroCountProb[kmer] / (float) (genusTotals[i]+1)); + for (int j = 0; j < wordGenusProb.size(); j++) { wordGenusProb[j].resize(genusNodes.size()); } + + int kmer, name; + vector numbers; numbers.resize(numKmers); + float prob; + vector zeroCountProb; zeroCountProb.resize(numKmers); + + //read file + for(int i=0;i length) { tempBuf = tempBuf.substr(0, length); } + delete buf4; + + istringstream iss (tempBuf,istringstream::in); + iss >> zeroCountProb[i] >> numbers[i]; } - //get probs for nonzero values - for (int i = 0; i < num[kmer]; i++) { - in >> name >> prob; - wordGenusProb[kmer][name] = prob; + MPI_File_close(&inMPI); + + for(int i=1;i length) { tempBuf = tempBuf.substr(0, length); } + delete buf4; + + istringstream iss (tempBuf,istringstream::in); + + iss >> kmer; + + //set them all to zero value + for (int i = 0; i < genusNodes.size(); i++) { + wordGenusProb[kmer][i] = log(zeroCountProb[kmer] / (float) (genusTotals[i]+1)); + } + + //get probs for nonzero values + for (int i = 0; i < numbers[kmer]; i++) { + iss >> name >> prob; + wordGenusProb[kmer][name] = prob; + } + } + MPI_File_close(&inMPI2); + #else + + in >> numKmers; gobble(in); - gobble(in); - } - in.close(); + //initialze probabilities + wordGenusProb.resize(numKmers); + + for (int j = 0; j < wordGenusProb.size(); j++) { wordGenusProb[j].resize(genusNodes.size()); } + + int kmer, name, count; count = 0; + vector num; num.resize(numKmers); + float prob; + vector zeroCountProb; zeroCountProb.resize(numKmers); + + while (inNum) { + inNum >> zeroCountProb[count] >> num[count]; + count++; + gobble(inNum); + } + inNum.close(); + + while(in) { + in >> kmer; + + //set them all to zero value + for (int i = 0; i < genusNodes.size(); i++) { + wordGenusProb[kmer][i] = log(zeroCountProb[kmer] / (float) (genusTotals[i]+1)); + } + + //get probs for nonzero values + for (int i = 0; i < num[kmer]; i++) { + in >> name >> prob; + wordGenusProb[kmer][name] = prob; + } + + gobble(in); + } + in.close(); + + #endif } catch(exception& e) { - errorOut(e, "Bayesian", "readProbFile"); + m->errorOut(e, "Bayesian", "readProbFile"); exit(1); } }