X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=classify.cpp;h=36179f471da4ac5d40b563e609feb2f8d3d32e6c;hb=1a968f34ae2d2680eaf189a197d1a21b8dfd6c03;hp=605c1f8f3c6470310a649cebfcc090141eb7d70d;hpb=63e089e0b3aad1741bab60119ed7ccc784dce347;p=mothur.git diff --git a/classify.cpp b/classify.cpp index 605c1f8..36179f4 100644 --- a/classify.cpp +++ b/classify.cpp @@ -12,113 +12,388 @@ #include "kmerdb.hpp" #include "suffixdb.hpp" #include "blastdb.hpp" +#include "distancedb.hpp" +#include "referencedb.h" /**************************************************************************************************/ - -Classify::Classify(string tfile, string tempFile, string method, int kmerSize, int gapOpen, int gapExtend, int match, int misMatch) : taxFile(tfile), templateFile(tempFile) { - try { +void Classify::generateDatabaseAndNames(string tfile, string tempFile, string method, int kmerSize, float gapOpen, float gapExtend, float match, float misMatch) { + try { + ReferenceDB* rdb = ReferenceDB::getInstance(); + + if (tfile == "saved") { tfile = rdb->getSavedTaxonomy(); } - readTaxonomy(taxFile); + taxFile = tfile; int numSeqs = 0; - //need to know number of template seqs for suffixdb - if (method == "suffix") { - ifstream inFASTA; - openInputFile(tempFile, inFASTA); - numSeqs = count(istreambuf_iterator(inFASTA),istreambuf_iterator(), '>'); - inFASTA.close(); - } + + if (tempFile == "saved") { + int start = time(NULL); + m->mothurOutEndLine(); m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory."); m->mothurOutEndLine(); + + numSeqs = rdb->referenceSeqs.size(); + templateFile = rdb->getSavedReference(); + tempFile = rdb->getSavedReference(); + + bool needToGenerate = true; + string kmerDBName; + if(method == "kmer") { + database = new KmerDB(tempFile, kmerSize); + + kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer"; + ifstream kmerFileTest(kmerDBName.c_str()); + if(kmerFileTest){ + bool GoodFile = m->checkReleaseVersion(kmerFileTest, m->getVersion()); + if (GoodFile) { needToGenerate = false; } + } + } + else if(method == "suffix") { database = new SuffixDB(numSeqs); } + else if(method == "blast") { database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch, "", threadID); } + else if(method == "distance") { database = new DistanceDB(); } + else { + m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8."); + m->mothurOutEndLine(); + database = new KmerDB(tempFile, 8); + } + + if (needToGenerate) { + for (int k = 0; k < rdb->referenceSeqs.size(); k++) { + Sequence temp(rdb->referenceSeqs[k].getName(), rdb->referenceSeqs[k].getAligned()); + names.push_back(temp.getName()); + database->addSequence(temp); + } + if ((method == "kmer") && (!shortcuts)) {;} //don't print + else {database->generateDB(); } + }else if ((method == "kmer") && (!needToGenerate)) { + ifstream kmerFileTest(kmerDBName.c_str()); + database->readKmerDB(kmerFileTest); + + for (int k = 0; k < rdb->referenceSeqs.size(); k++) { + names.push_back(rdb->referenceSeqs[k].getName()); + } + } + + database->setNumSeqs(numSeqs); + + m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(rdb->referenceSeqs.size()) + " sequences and generate the search databases.");m->mothurOutEndLine(); + + }else { + + templateFile = tempFile; + + int start = time(NULL); + + m->mothurOut("Generating search database... "); cout.flush(); + #ifdef USE_MPI + int pid, processors; + vector positions; + int tag = 2001; + + MPI_Status status; + MPI_File inMPI; + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + MPI_Comm_size(MPI_COMM_WORLD, &processors); - mothurOut("Generating search database... "); cout.flush(); + //char* inFileName = new char[tempFile.length()]; + //memcpy(inFileName, tempFile.c_str(), tempFile.length()); - bool needToGenerate = true; - string kmerDBName; - if(method == "kmer") { - database = new KmerDB(tempFile, kmerSize); - - kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer"; - ifstream kmerFileTest(kmerDBName.c_str()); - if(kmerFileTest){ needToGenerate = false; } - } - else if(method == "suffix") { database = new SuffixDB(numSeqs); } - else if(method == "blast") { database = new BlastDB(gapOpen, gapExtend, match, misMatch); } - else { - mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8."); - mothurOutEndLine(); - database = new KmerDB(tempFile, 8); - } - - if (needToGenerate) { - ifstream fastaFile; - openInputFile(tempFile, fastaFile); - - while (!fastaFile.eof()) { - Sequence temp(fastaFile); - gobble(fastaFile); + char inFileName[1024]; + strcpy(inFileName, tempFile.c_str()); + + MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI); //comm, filename, mode, info, filepointer + //delete inFileName; + + if (pid == 0) { //only one process needs to scan file + positions = m->setFilePosFasta(tempFile, numSeqs); //fills MPIPos, returns numSeqs + + //send file positions to all processes + for(int i = 1; i < processors; i++) { + MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD); + MPI_Send(&positions[0], (numSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD); + } + }else{ + MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status); + positions.resize(numSeqs+1); + MPI_Recv(&positions[0], (numSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status); + } + + //create database + if(method == "kmer") { database = new KmerDB(tempFile, kmerSize); } + else if(method == "suffix") { database = new SuffixDB(numSeqs); } + else if(method == "blast") { database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch, "", pid); } + else if(method == "distance") { database = new DistanceDB(); } + else { + m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8."); m->mothurOutEndLine(); + database = new KmerDB(tempFile, 8); + } + + //read file + for(int i=0;i length) { tempBuf = tempBuf.substr(0, length); } + delete buf4; + istringstream iss (tempBuf,istringstream::in); + + Sequence temp(iss); + if (temp.getName() != "") { + if (rdb->save) { rdb->referenceSeqs.push_back(temp); } + names.push_back(temp.getName()); + database->addSequence(temp); + } + } + + database->generateDB(); + MPI_File_close(&inMPI); + MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case + #else - names.push_back(temp.getName()); - database->addSequence(temp); + //need to know number of template seqs for suffixdb + if (method == "suffix") { + ifstream inFASTA; + m->openInputFile(tempFile, inFASTA); + m->getNumSeqs(inFASTA, numSeqs); + inFASTA.close(); } - fastaFile.close(); - database->generateDB(); + bool needToGenerate = true; + string kmerDBName; + if(method == "kmer") { + database = new KmerDB(tempFile, kmerSize); + + kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer"; + ifstream kmerFileTest(kmerDBName.c_str()); + if(kmerFileTest){ + bool GoodFile = m->checkReleaseVersion(kmerFileTest, m->getVersion()); + if (GoodFile) { needToGenerate = false; } + } + } + else if(method == "suffix") { database = new SuffixDB(numSeqs); } + else if(method == "blast") { database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch, "", threadID); } + else if(method == "distance") { database = new DistanceDB(); } + else { + m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8."); + m->mothurOutEndLine(); + database = new KmerDB(tempFile, 8); + } - }else if ((method == "kmer") && (!needToGenerate)) { - ifstream kmerFileTest(kmerDBName.c_str()); - database->readKmerDB(kmerFileTest); + if (needToGenerate) { + ifstream fastaFile; + m->openInputFile(tempFile, fastaFile); + + while (!fastaFile.eof()) { + Sequence temp(fastaFile); + m->gobble(fastaFile); + + if (rdb->save) { rdb->referenceSeqs.push_back(temp); } + + names.push_back(temp.getName()); + + database->addSequence(temp); + } + fastaFile.close(); + + if ((method == "kmer") && (!shortcuts)) {;} //don't print + else {database->generateDB(); } + + }else if ((method == "kmer") && (!needToGenerate)) { + ifstream kmerFileTest(kmerDBName.c_str()); + database->readKmerDB(kmerFileTest); - ifstream fastaFile; - openInputFile(tempFile, fastaFile); + ifstream fastaFile; + m->openInputFile(tempFile, fastaFile); + + while (!fastaFile.eof()) { + Sequence temp(fastaFile); + m->gobble(fastaFile); + + if (rdb->save) { rdb->referenceSeqs.push_back(temp); } + names.push_back(temp.getName()); + } + fastaFile.close(); + } + #endif + + database->setNumSeqs(names.size()); + + m->mothurOut("DONE."); m->mothurOutEndLine(); + m->mothurOut("It took " + toString(time(NULL) - start) + " seconds generate search database. "); m->mothurOutEndLine(); + } + + readTaxonomy(taxFile); + + //sanity check + bool okay = phyloTree->ErrorCheck(names); + + if (!okay) { m->control_pressed = true; } + } + catch(exception& e) { + m->errorOut(e, "Classify", "generateDatabaseAndNames"); + exit(1); + } +} +/**************************************************************************************************/ +Classify::Classify() { m = MothurOut::getInstance(); database = NULL; phyloTree=NULL; flipped=false; } +/**************************************************************************************************/ + +int Classify::readTaxonomy(string file) { + try { + + phyloTree = new PhyloTree(); + string name, taxInfo; - while (!fastaFile.eof()) { - Sequence temp(fastaFile); - gobble(fastaFile); + m->mothurOutEndLine(); + m->mothurOut("Reading in the " + file + " taxonomy...\t"); cout.flush(); + if (m->debug) { m->mothurOut("[DEBUG]: Taxonomies read in...\n"); } + +#ifdef USE_MPI + int pid, num, processors; + vector positions; + int tag = 2001; + + MPI_Status status; + MPI_File inMPI; + MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are + MPI_Comm_size(MPI_COMM_WORLD, &processors); + + char inFileName[1024]; + strcpy(inFileName, file.c_str()); + + MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI); //comm, filename, mode, info, filepointer + //delete inFileName; + + if (pid == 0) { + positions = m->setFilePosEachLine(file, num); - names.push_back(temp.getName()); + //send file positions to all processes + 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); } - fastaFile.close(); + }else{ + 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); } + + //read file + for(int i=0;i length) { tempBuf = tempBuf.substr(0, length); } + delete buf4; + + istringstream iss (tempBuf,istringstream::in); + iss >> name; m->gobble(iss); + iss >> taxInfo; + if (m->debug) { m->mothurOut("[DEBUG]: name = " + name + " tax = " + taxInfo + "\n"); } + if (m->inUsersGroups(name, names)) { + taxonomy[name] = taxInfo; + phyloTree->addSeqToTree(name, taxInfo); + }else { + m->mothurOut("[WARNING]: " + name + " is in your taxonomy file and not in your reference file, ignoring.\n"); + } + } - database->setNumSeqs(names.size()); + MPI_File_close(&inMPI); + MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case +#else + + taxonomy.clear(); + m->readTax(file, taxonomy); + map tempTaxonomy; + for (map::iterator itTax = taxonomy.begin(); itTax != taxonomy.end(); itTax++) { + if (m->inUsersGroups(itTax->first, names)) { + phyloTree->addSeqToTree(itTax->first, itTax->second); + tempTaxonomy[itTax->first] = itTax->second; + }else { + m->mothurOut("[WARNING]: " + itTax->first + " is in your taxonomy file and not in your reference file, ignoring.\n"); + } + } + taxonomy = tempTaxonomy; +#endif + phyloTree->assignHeirarchyIDs(0); - mothurOut("DONE."); mothurOutEndLine(); - + phyloTree->setUp(file); + + m->mothurOut("DONE."); + m->mothurOutEndLine(); cout.flush(); + + return phyloTree->getNumSeqs(); + } catch(exception& e) { - errorOut(e, "Classify", "Classify"); + m->errorOut(e, "Classify", "readTaxonomy"); exit(1); } } /**************************************************************************************************/ -void Classify::readTaxonomy(string file) { +vector Classify::parseTax(string tax) { try { - - ifstream inTax; - openInputFile(file, inTax); - - mothurOutEndLine(); - mothurOut("Reading in the " + file + " taxonomy...\t"); cout.flush(); + vector taxons; - string name, taxInfo; - //read template seqs and save - while (!inTax.eof()) { - inTax >> name >> taxInfo; + tax = tax.substr(0, tax.length()-1); //get rid of last ';' + + //parse taxonomy + string individual; + while (tax.find_first_of(';') != -1) { + individual = tax.substr(0,tax.find_first_of(';')); + tax = tax.substr(tax.find_first_of(';')+1, tax.length()); + taxons.push_back(individual); - taxonomy[name] = taxInfo; - - gobble(inTax); } - inTax.close(); - - mothurOut("DONE."); - mothurOutEndLine(); cout.flush(); - + //get last one + taxons.push_back(tax); + + return taxons; } catch(exception& e) { - errorOut(e, "Classify", "readTaxonomy"); + m->errorOut(e, "Classify", "parseTax"); exit(1); } } /**************************************************************************************************/ +double Classify::getLogExpSum(vector probabilities, int& maxIndex){ + try { + // http://jblevins.org/notes/log-sum-exp + + double maxProb = probabilities[0]; + maxIndex = 0; + + int numProbs = (int)probabilities.size(); + + for(int i=1;i= maxProb){ + maxProb = probabilities[i]; + maxIndex = i; + } + } + + double probSum = 0.0000; + + for(int i=0;ierrorOut(e, "Classify", "getLogExpSum"); + exit(1); + } +} + +/**************************************************************************************************/ +