classify.cpp

   1 /*
   2  *  classify.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 11/3/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "classify.h"
  11 #include "sequence.hpp"
  12 #include "kmerdb.hpp"
  13 #include "suffixdb.hpp"
  14 #include "blastdb.hpp"
  15 #include "distancedb.hpp"
  16 #include "referencedb.h"
  17
  18 /**************************************************************************************************/
  19 void Classify::generateDatabaseAndNames(string tfile, string tempFile, string method, int kmerSize, float gapOpen, float gapExtend, float match, float misMatch)  {
  20         try {
  21                 ReferenceDB* rdb = ReferenceDB::getInstance();
  22
  23                 if (tfile == "saved") { tfile = rdb->getSavedTaxonomy(); }
  24
  25                 taxFile = tfile;
  26                 readTaxonomy(taxFile);
  27                 int numSeqs = 0;
  28
  29                 if (tempFile == "saved") {
  30                         int start = time(NULL);
  31                         m->mothurOutEndLine();  m->mothurOut("Using sequences from " + rdb->getSavedReference() + " that are saved in memory.");        m->mothurOutEndLine();
  32
  33                         numSeqs = rdb->referenceSeqs.size();
  34                         templateFile = rdb->getSavedReference();
  35                         tempFile = rdb->getSavedReference();
  36
  37                         bool needToGenerate = true;
  38                         string kmerDBName;
  39                         if(method == "kmer")                    {
  40                                 database = new KmerDB(tempFile, kmerSize);
  41
  42                                 kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
  43                                 ifstream kmerFileTest(kmerDBName.c_str());
  44                                 if(kmerFileTest){
  45                                         bool GoodFile = m->checkReleaseVersion(kmerFileTest, m->getVersion());
  46                                         if (GoodFile) {  needToGenerate = false;        }
  47                                 }
  48                         }
  49                         else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
  50                         else if(method == "blast")              {       database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch, "", threadID);     }
  51                         else if(method == "distance")   {       database = new DistanceDB();    }
  52                         else {
  53                                 m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
  54                                 m->mothurOutEndLine();
  55                                 database = new KmerDB(tempFile, 8);
  56                         }
  57
  58                         if (needToGenerate) {
  59                                 for (int k = 0; k < rdb->referenceSeqs.size(); k++) {
  60                                         Sequence temp(rdb->referenceSeqs[k].getName(), rdb->referenceSeqs[k].getAligned());
  61                                         names.push_back(temp.getName());
  62                                         database->addSequence(temp);
  63                                 }
  64                                 database->generateDB();
  65                         }else if ((method == "kmer") && (!needToGenerate)) {
  66                                 ifstream kmerFileTest(kmerDBName.c_str());
  67                                 database->readKmerDB(kmerFileTest);
  68
  69                                 for (int k = 0; k < rdb->referenceSeqs.size(); k++) {
  70                                         names.push_back(rdb->referenceSeqs[k].getName());
  71                                 }
  72                         }
  73
  74                         database->setNumSeqs(numSeqs);
  75
  76                         //sanity check
  77                         bool okay = phyloTree->ErrorCheck(names);
  78
  79                         if (!okay) { m->control_pressed = true; }
  80
  81                         m->mothurOut("It took " + toString(time(NULL) - start) + " to load " + toString(rdb->referenceSeqs.size()) + " sequences and generate the search databases.");m->mothurOutEndLine();
  82
  83                 }else {
  84
  85                         templateFile = tempFile;
  86
  87                         int start = time(NULL);
  88
  89                         m->mothurOut("Generating search database...    "); cout.flush();
  90         #ifdef USE_MPI
  91                                 int pid, processors;
  92                                 vector<unsigned long long> positions;
  93                                 int tag = 2001;
  94
  95                                 MPI_Status status;
  96                                 MPI_File inMPI;
  97                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
  98                                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
  99
 100                                 //char* inFileName = new char[tempFile.length()];
 101                                 //memcpy(inFileName, tempFile.c_str(), tempFile.length());
 102
 103                                 char inFileName[1024];
 104                                 strcpy(inFileName, tempFile.c_str());
 105
 106                                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 107                                 //delete inFileName;
 108
 109                                 if (pid == 0) { //only one process needs to scan file
 110                                         positions = m->setFilePosFasta(tempFile, numSeqs); //fills MPIPos, returns numSeqs
 111
 112                                         //send file positions to all processes
 113                                         for(int i = 1; i < processors; i++) {
 114                                                 MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 115                                                 MPI_Send(&positions[0], (numSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
 116                                         }
 117                                 }else{
 118                                         MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 119                                         positions.resize(numSeqs+1);
 120                                         MPI_Recv(&positions[0], (numSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
 121                                 }
 122
 123                                 //create database
 124                                 if(method == "kmer")                    {       database = new KmerDB(tempFile, kmerSize);                      }
 125                                 else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
 126                                 else if(method == "blast")              {       database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch, "", pid);  }
 127                                 else if(method == "distance")   {       database = new DistanceDB();    }
 128                                 else {
 129                                         m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8."); m->mothurOutEndLine();
 130                                         database = new KmerDB(tempFile, 8);
 131                                 }
 132
 133                                 //read file
 134                                 for(int i=0;i<numSeqs;i++){
 135                                         //read next sequence
 136                                         int length = positions[i+1] - positions[i];
 137                                         char* buf4 = new char[length];
 138                                         MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
 139
 140                                         string tempBuf = buf4;
 141                                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
 142                                         delete buf4;
 143                                         istringstream iss (tempBuf,istringstream::in);
 144
 145                                         Sequence temp(iss);
 146                                         if (temp.getName() != "") {
 147                                                 if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
 148                                                 names.push_back(temp.getName());
 149                                                 database->addSequence(temp);
 150                                         }
 151                                 }
 152
 153                                 database->generateDB();
 154                                 MPI_File_close(&inMPI);
 155                                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 156                 #else
 157
 158                         //need to know number of template seqs for suffixdb
 159                         if (method == "suffix") {
 160                                 ifstream inFASTA;
 161                                 m->openInputFile(tempFile, inFASTA);
 162                                 m->getNumSeqs(inFASTA, numSeqs);
 163                                 inFASTA.close();
 164                         }
 165
 166                         bool needToGenerate = true;
 167                         string kmerDBName;
 168                         if(method == "kmer")                    {
 169                                 database = new KmerDB(tempFile, kmerSize);
 170
 171                                 kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
 172                                 ifstream kmerFileTest(kmerDBName.c_str());
 173                                 if(kmerFileTest){
 174                                         bool GoodFile = m->checkReleaseVersion(kmerFileTest, m->getVersion());
 175                                         if (GoodFile) {  needToGenerate = false;        }
 176                                 }
 177                         }
 178                         else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
 179                         else if(method == "blast")              {       database = new BlastDB(tempFile.substr(0,tempFile.find_last_of(".")+1), gapOpen, gapExtend, match, misMatch, "", threadID);     }
 180                         else if(method == "distance")   {       database = new DistanceDB();    }
 181                         else {
 182                                 m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
 183                                 m->mothurOutEndLine();
 184                                 database = new KmerDB(tempFile, 8);
 185                         }
 186
 187                         if (needToGenerate) {
 188                                 ifstream fastaFile;
 189                                 m->openInputFile(tempFile, fastaFile);
 190
 191                                 while (!fastaFile.eof()) {
 192                                         Sequence temp(fastaFile);
 193                                         m->gobble(fastaFile);
 194
 195                                         if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
 196
 197                                         names.push_back(temp.getName());
 198
 199                                         database->addSequence(temp);
 200                                 }
 201                                 fastaFile.close();
 202
 203                 if ((method == "kmer") && (!shortcuts)) {;} //don't print
 204                 else {database->generateDB(); }
 205
 206                         }else if ((method == "kmer") && (!needToGenerate)) {
 207                                 ifstream kmerFileTest(kmerDBName.c_str());
 208                                 database->readKmerDB(kmerFileTest);
 209
 210                                 ifstream fastaFile;
 211                                 m->openInputFile(tempFile, fastaFile);
 212
 213                                 while (!fastaFile.eof()) {
 214                                         Sequence temp(fastaFile);
 215                                         m->gobble(fastaFile);
 216
 217                                         if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
 218                                         names.push_back(temp.getName());
 219                                 }
 220                                 fastaFile.close();
 221                         }
 222         #endif
 223
 224                         database->setNumSeqs(names.size());
 225
 226                         //sanity check
 227                         bool okay = phyloTree->ErrorCheck(names);
 228
 229                         if (!okay) { m->control_pressed = true; }
 230
 231                         m->mothurOut("DONE."); m->mothurOutEndLine();
 232                         m->mothurOut("It took " + toString(time(NULL) - start) + " seconds generate search database. "); m->mothurOutEndLine();
 233                 }
 234
 235         }
 236         catch(exception& e) {
 237                 m->errorOut(e, "Classify", "generateDatabaseAndNames");
 238                 exit(1);
 239         }
 240 }
 241 /**************************************************************************************************/
 242 Classify::Classify() {          m = MothurOut::getInstance();   database = NULL;        flipped=false; }
 243 /**************************************************************************************************/
 244
 245 int Classify::readTaxonomy(string file) {
 246         try {
 247
 248                 phyloTree = new PhyloTree();
 249                 string name, taxInfo;
 250
 251                 m->mothurOutEndLine();
 252                 m->mothurOut("Reading in the " + file + " taxonomy...\t");      cout.flush();
 253         if (m->debug) { m->mothurOut("[DEBUG]: Taxonomies read in...\n"); }
 254
 255 #ifdef USE_MPI
 256                 int pid, num, processors;
 257                 vector<unsigned long long> positions;
 258                 int tag = 2001;
 259
 260                 MPI_Status status;
 261                 MPI_File inMPI;
 262                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 263                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
 264
 265                 char inFileName[1024];
 266                 strcpy(inFileName, file.c_str());
 267
 268                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 269                 //delete inFileName;
 270
 271                 if (pid == 0) {
 272                         positions = m->setFilePosEachLine(file, num);
 273
 274                         //send file positions to all processes
 275                         for(int i = 1; i < processors; i++) {
 276                                 MPI_Send(&num, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 277                                 MPI_Send(&positions[0], (num+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
 278                         }
 279                 }else{
 280                         MPI_Recv(&num, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 281                         positions.resize(num+1);
 282                         MPI_Recv(&positions[0], (num+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
 283                 }
 284
 285                 //read file
 286                 for(int i=0;i<num;i++){
 287                         //read next sequence
 288                         int length = positions[i+1] - positions[i];
 289                         char* buf4 = new char[length];
 290
 291                         MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
 292
 293                         string tempBuf = buf4;
 294                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
 295                         delete buf4;
 296
 297                         istringstream iss (tempBuf,istringstream::in);
 298                         iss >> name; m->gobble(iss);
 299             iss >> taxInfo;
 300             if (m->debug) { m->mothurOut("[DEBUG]: name = " + name + " tax = " + taxInfo + "\n"); }
 301                         taxonomy[name] = taxInfo;
 302                         phyloTree->addSeqToTree(name, taxInfo);
 303                 }
 304
 305                 MPI_File_close(&inMPI);
 306                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 307 #else
 308
 309         taxonomy.clear();
 310         m->readTax(file, taxonomy);
 311         for (map<string, string>::iterator itTax = taxonomy.begin(); itTax != taxonomy.end(); itTax++) {  phyloTree->addSeqToTree(itTax->first, itTax->second);  }
 312 #endif
 313                 phyloTree->assignHeirarchyIDs(0);
 314
 315                 phyloTree->setUp(file);
 316
 317                 m->mothurOut("DONE.");
 318                 m->mothurOutEndLine();  cout.flush();
 319
 320                 return phyloTree->getNumSeqs();
 321
 322         }
 323         catch(exception& e) {
 324                 m->errorOut(e, "Classify", "readTaxonomy");
 325                 exit(1);
 326         }
 327 }
 328 /**************************************************************************************************/
 329
 330 vector<string> Classify::parseTax(string tax) {
 331         try {
 332                 vector<string> taxons;
 333
 334                 tax = tax.substr(0, tax.length()-1);  //get rid of last ';'
 335
 336                 //parse taxonomy
 337                 string individual;
 338                 while (tax.find_first_of(';') != -1) {
 339                         individual = tax.substr(0,tax.find_first_of(';'));
 340                         tax = tax.substr(tax.find_first_of(';')+1, tax.length());
 341                         taxons.push_back(individual);
 342
 343                 }
 344                 //get last one
 345                 taxons.push_back(tax);
 346
 347                 return taxons;
 348         }
 349         catch(exception& e) {
 350                 m->errorOut(e, "Classify", "parseTax");
 351                 exit(1);
 352         }
 353 }
 354 /**************************************************************************************************/
 355