X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bayesian.cpp;h=9ae0777912ba21ad231540eed5658c09298b4ece;hb=1e8d08e96f4fe99604a6b3502568de464bf60891;hp=ae96a220ae2646c59a45a36152ed2d0b099c0c64;hpb=2df35fdeea85f574630d75b11fb5b08c39aec31a;p=mothur.git diff --git a/bayesian.cpp b/bayesian.cpp index ae96a22..9ae0777 100644 --- a/bayesian.cpp +++ b/bayesian.cpp @@ -15,10 +15,10 @@ Bayesian::Bayesian(string tfile, string tempFile, string method, int ksize, int cutoff, int i) : Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { try { - + /************calculate the probablity that each word will be in a specific taxonomy*************/ string tfileroot = tfile.substr(0,tfile.find_last_of(".")+1); - string tempfileroot = getRootName(getSimpleName(tempFile)); + string tempfileroot = m->getRootName(m->getSimpleName(tempFile)); string phyloTreeName = tfileroot + "tree.train"; string phyloTreeSumName = tfileroot + "tree.sum"; string probFileName = tfileroot + tempfileroot + char('0'+ kmerSize) + "mer.prob"; @@ -34,7 +34,13 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { int start = time(NULL); - if(probFileTest && probFileTest2 && phyloTreeTest && probFileTest3){ + //if they are there make sure they were created after this release date + bool FilesGood = false; + if(probFileTest && probFileTest2 && phyloTreeTest && probFileTest3){ + FilesGood = checkReleaseDate(probFileTest, probFileTest2, phyloTreeTest, probFileTest3); + } + + if(probFileTest && probFileTest2 && phyloTreeTest && probFileTest3 && FilesGood){ m->mothurOut("Reading template taxonomy... "); cout.flush(); phyloTree = new PhyloTree(phyloTreeTest, phyloTreeName); @@ -70,9 +76,9 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { //initialze probabilities wordGenusProb.resize(numKmers); - + //cout << numKmers << '\t' << genusNodes.size() << endl; for (int j = 0; j < wordGenusProb.size(); j++) { wordGenusProb[j].resize(genusNodes.size()); } - + //cout << numKmers << '\t' << genusNodes.size() << endl; ofstream out; ofstream out2; @@ -84,11 +90,17 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { #endif - openOutputFile(probFileName, out); + m->openOutputFile(probFileName, out); + + //output mothur version + out << "#" << m->getVersion() << endl; out << numKmers << endl; - openOutputFile(probFileName2, out2); + m->openOutputFile(probFileName2, out2); + + //output mothur version + out2 << "#" << m->getVersion() << endl; #ifdef USE_MPI } @@ -193,6 +205,7 @@ Classify(), kmerSize(ksize), confidenceThreshold(cutoff), iters(i) { /**************************************************************************************************/ Bayesian::~Bayesian() { try { + delete phyloTree; if (database != NULL) { delete database; } } @@ -211,7 +224,7 @@ 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; @@ -223,14 +236,16 @@ string Bayesian::getTaxonomy(Sequence* seq) { 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; } -//cout << seq->getName() << '\t' << index << endl; + //bootstrap - to set confidenceScore int numToSelect = queryKmers.size() / 8; + tax = bootstrapResults(queryKmers, index, numToSelect); - + return tax; } catch(exception& e) { @@ -243,6 +258,17 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { try { map confidenceScores; + + //initialize confidences to 0 + int seqIndex = tax; + TaxNode seq = phyloTree->get(tax); + confidenceScores[tax] = 0; + + while (seq.level != 0) { //while you are not at the root + seqIndex = seq.parent; + confidenceScores[seqIndex] = 0; + seq = phyloTree->get(seq.parent); + } map::iterator itBoot; map::iterator itBoot2; @@ -252,7 +278,6 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { if (m->control_pressed) { return "control"; } vector temp; - for (int j = 0; j < numToSelect; j++) { int index = int(rand() % kmers.size()); @@ -262,17 +287,15 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { //get taxonomy int newTax = getMostProbableTaxonomy(temp); + //int newTax = 1; TaxNode taxonomyTemp = phyloTree->get(newTax); - + //add to confidence results 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 - if (itBoot2 == confidenceScores.end()) { //not already in confidence scores - confidenceScores[newTax] = 1; - }else{ - confidenceScores[newTax]++; + if (itBoot2 != confidenceScores.end()) { //this is a classification we need a confidence for + (itBoot2->second)++; } newTax = taxonomyTemp.parent; @@ -293,7 +316,7 @@ string Bayesian::bootstrapResults(vector kmers, int tax, int numToSelect) { int confidence = 0; if (itBoot2 != confidenceScores.end()) { //already in confidence scores - confidence = confidenceScores[seqTaxIndex]; + confidence = itBoot2->second; } if (((confidence/(float)iters) * 100) >= confidenceThreshold) { @@ -327,14 +350,14 @@ int Bayesian::getMostProbableTaxonomy(vector queryKmer) { for (int i = 0; i < queryKmer.size(); i++) { prob += wordGenusProb[queryKmer[i]][k]; } - + //is this the taxonomy with the greatest probability? if (prob > maxProbability) { indexofGenus = genusNodes[k]; maxProbability = prob; } } - +// cout << phyloTree->get(indexofGenus).name << '\t' << maxProbability << endl; return indexofGenus; } catch(exception& e) { @@ -396,8 +419,8 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string 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); - positions2 = setFilePosEachLine(inName, num2); + positions = m->setFilePosEachLine(inNumName, num); + positions2 = m->setFilePosEachLine(inName, num2); for(int i = 1; i < processors; i++) { MPI_Send(&num, 1, MPI_INT, i, tag, MPI_COMM_WORLD); @@ -416,12 +439,19 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string positions2.resize(num2+1); MPI_Recv(&positions2[0], (num2+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status); } - - //read numKmers + + //read version int length = positions2[1] - positions2[0]; + char* buf5 = new char[length]; + + MPI_File_read_at(inMPI2, positions2[0], buf5, length, MPI_CHAR, &status); + delete buf5; + + //read numKmers + length = positions2[2] - positions2[1]; char* buf = new char[length]; - MPI_File_read_at(inMPI2, positions2[0], buf, length, MPI_CHAR, &status); + MPI_File_read_at(inMPI2, positions2[1], buf, length, MPI_CHAR, &status); string tempBuf = buf; if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); } @@ -438,10 +468,17 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string int kmer, name; vector numbers; numbers.resize(numKmers); float prob; - vector zeroCountProb; zeroCountProb.resize(numKmers); + vector zeroCountProb; zeroCountProb.resize(numKmers); + + //read version + length = positions[1] - positions[0]; + char* buf6 = new char[length]; + MPI_File_read_at(inMPI2, positions[0], buf6, length, MPI_CHAR, &status); + delete buf6; + //read file - for(int i=0;i> numKmers; gobble(in); + //read version + string line = m->getline(in); m->gobble(in); + + in >> numKmers; m->gobble(in); //initialze probabilities wordGenusProb.resize(numKmers); @@ -500,11 +539,14 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string vector num; num.resize(numKmers); float prob; vector zeroCountProb; zeroCountProb.resize(numKmers); - + + //read version + string line2 = m->getline(inNum); m->gobble(inNum); + while (inNum) { inNum >> zeroCountProb[count] >> num[count]; count++; - gobble(inNum); + m->gobble(inNum); } inNum.close(); @@ -522,7 +564,7 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string wordGenusProb[kmer][name] = prob; } - gobble(in); + m->gobble(in); } in.close(); @@ -534,6 +576,61 @@ void Bayesian::readProbFile(ifstream& in, ifstream& inNum, string inName, string } } /**************************************************************************************************/ +bool Bayesian::checkReleaseDate(ifstream& file1, ifstream& file2, ifstream& file3, ifstream& file4) { + try { + + bool good = true; + + vector lines; + lines.push_back(m->getline(file1)); + lines.push_back(m->getline(file2)); + lines.push_back(m->getline(file3)); + lines.push_back(m->getline(file4)); + + //before we added this check + if ((lines[0][0] != '#') || (lines[1][0] != '#') || (lines[2][0] != '#') || (lines[3][0] != '#')) { good = false; } + else { + //rip off # + for (int i = 0; i < lines.size(); i++) { lines[i] = lines[i].substr(1); } + + //get mothurs current version + string version = m->getVersion(); + + vector versionVector; + m->splitAtChar(version, versionVector, '.'); + + //check each files version + for (int i = 0; i < lines.size(); i++) { + vector linesVector; + m->splitAtChar(lines[i], linesVector, '.'); + + if (versionVector.size() != linesVector.size()) { good = false; break; } + else { + for (int j = 0; j < versionVector.size(); j++) { + int num1, num2; + convert(versionVector[j], num1); + convert(linesVector[j], num2); + + //if mothurs version is newer than this files version, then we want to remake it + if (num1 > num2) { good = false; break; } + } + } + + if (!good) { break; } + } + } + + if (!good) { file1.close(); file2.close(); file3.close(); file4.close(); } + else { file1.seekg(0); file2.seekg(0); file3.seekg(0); file4.seekg(0); } + + return good; + } + catch(exception& e) { + m->errorOut(e, "Bayesian", "checkReleaseDate"); + exit(1); + } +} +/**************************************************************************************************/