X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bayesian.cpp;h=e5ab89e6216426e2cd9ca1c8a7eca16db203cd00;hb=9013e13ecfb2fda3c2664a76f76cc99b8c7fa74c;hp=22eab72b60574679f61eca1f02990a3178aecd75;hpb=4a2d841cb97fb02351022efe9d7068b1dc212bf9;p=mothur.git diff --git a/bayesian.cpp b/bayesian.cpp index 22eab72..e5ab89e 100644 --- a/bayesian.cpp +++ b/bayesian.cpp @@ -17,9 +17,11 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { try { /************calculate the probablity that each word will be in a specific taxonomy*************/ - 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"; + string tfileroot = tfile.substr(0,tfile.find_last_of(".")+1); + string tempfileroot = getRootName(getSimpleName(tempFile)); + string phyloTreeName = tfileroot + "tree.train"; + string probFileName = tfileroot + tempfileroot + char('0'+ kmerSize) + "mer.prob"; + string probFileName2 = tfileroot + tempfileroot + char('0'+ kmerSize) + "mer.numNonZero"; ofstream out; ofstream out2; @@ -48,125 +50,134 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { //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); + //prevents errors caused by creating shortcut files if you had an error in the sanity check. + if (m->control_pressed) { remove(phyloTreeName.c_str()); remove(probFileName.c_str()); remove(probFileName2.c_str()); } + else{ + 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; - #ifdef USE_MPI - } - #endif - + //initialze probabilities + wordGenusProb.resize(numKmers); - //for each word - for (int i = 0; i < numKmers; i++) { - if (m->control_pressed) { break; } + for (int j = 0; j < wordGenusProb.size(); j++) { wordGenusProb[j].resize(genusNodes.size()); } + + ofstream out; + ofstream out2; #ifdef USE_MPI + int pid; MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are if (pid == 0) { #endif - out << i << '\t'; - #ifdef USE_MPI - } - #endif + openOutputFile(probFileName, out); - vector seqsWithWordi = database->getSequencesWithKmer(i); + out << numKmers << endl; - map count; - for (int k = 0; k < genusNodes.size(); k++) { count[genusNodes[k]] = 0; } - - //for each sequence with that word - for (int j = 0; j < seqsWithWordi.size(); j++) { - int temp = phyloTree->getIndex(names[seqsWithWordi[j]]); - count[temp]++; //increment count of seq in this genus who have this word - } + openOutputFile(probFileName2, out2); - //probabilityInTemplate = (# of seqs with that word in template + 0.05) / (total number of seqs in template + 1); - float probabilityInTemplate = (seqsWithWordi.size() + 0.50) / (float) (names.size() + 1); + #ifdef USE_MPI + } + #endif + - int numNotZero = 0; - 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) { - #ifdef USE_MPI - MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are - if (pid == 0) { - #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 - out << k << '\t' << wordGenusProb[i][k] << '\t'; + if (pid == 0) { + #endif + + out << i << '\t'; + + #ifdef USE_MPI + } + #endif + + 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 + for (int j = 0; j < seqsWithWordi.size(); j++) { + int temp = phyloTree->getIndex(names[seqsWithWordi[j]]); + count[temp]++; //increment count of seq in this genus who have this word + } + + //probabilityInTemplate = (# of seqs with that word in template + 0.05) / (total number of seqs in template + 1); + float probabilityInTemplate = (seqsWithWordi.size() + 0.50) / (float) (names.size() + 1); + + int numNotZero = 0; + 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) { + #ifdef USE_MPI + int pid; + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are - #ifdef USE_MPI - } - #endif + if (pid == 0) { + #endif + + out << k << '\t' << wordGenusProb[i][k] << '\t'; + + #ifdef USE_MPI + } + #endif - numNotZero++; + 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 << endl; - out2 << probabilityInTemplate << '\t' << numNotZero << endl; + 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); } - - #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); } m->mothurOut("DONE."); m->mothurOutEndLine(); @@ -213,7 +224,7 @@ string Bayesian::getTaxonomy(Sequence* seq) { int index = getMostProbableTaxonomy(queryKmers); if (m->control_pressed) { return tax; } - +//cout << seq->getName() << '\t' << index << endl; //bootstrap - to set confidenceScore int numToSelect = queryKmers.size() / 8; tax = bootstrapResults(queryKmers, index, numToSelect); @@ -249,10 +260,10 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { //get taxonomy int newTax = getMostProbableTaxonomy(temp); - TaxNode taxonomy = phyloTree->get(newTax); - + TaxNode taxonomyTemp = phyloTree->get(newTax); + //add to confidence results - while (taxonomy.level != 0) { //while you are not at the root + while (taxonomyTemp.level != 0) { //while you are not at the root itBoot2 = confidenceScores.find(newTax); //is this a classification we already have a count on @@ -262,8 +273,8 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { confidenceScores[newTax]++; } - newTax = taxonomy.parent; - taxonomy = phyloTree->get(taxonomy.parent); + newTax = taxonomyTemp.parent; + taxonomyTemp = phyloTree->get(newTax); } } @@ -362,14 +373,16 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string #ifdef USE_MPI - int pid, num, num2; - vector positions; - vector positions2; + int pid, num, num2, processors; + 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 + MPI_Comm_size(MPI_COMM_WORLD, &processors); + int tag = 2001; char inFileName[1024]; strcpy(inFileName, inNumName.c_str()); @@ -382,26 +395,24 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string 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 + for(int i = 1; i < processors; i++) { + MPI_Send(&num, 1, MPI_INT, i, tag, MPI_COMM_WORLD); + MPI_Send(&positions[0], (num+1), MPI_LONG, i, tag, MPI_COMM_WORLD); + + MPI_Send(&num2, 1, MPI_INT, i, tag, MPI_COMM_WORLD); + MPI_Send(&positions2[0], (num2+1), MPI_LONG, i, tag, MPI_COMM_WORLD); + } }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_Recv(&num, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status); + positions.resize(num+1); + MPI_Recv(&positions[0], (num+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status); - 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 - + MPI_Recv(&num2, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status); + positions2.resize(num2+1); + MPI_Recv(&positions2[0], (num2+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status); } //read numKmers @@ -473,6 +484,7 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string } MPI_File_close(&inMPI2); + MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case #else in >> numKmers; gobble(in);