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 int> 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, "");       }
 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                                 database->generateDB();
 204
 205                         }else if ((method == "kmer") && (!needToGenerate)) {
 206                                 ifstream kmerFileTest(kmerDBName.c_str());
 207                                 database->readKmerDB(kmerFileTest);
 208
 209                                 ifstream fastaFile;
 210                                 m->openInputFile(tempFile, fastaFile);
 211
 212                                 while (!fastaFile.eof()) {
 213                                         Sequence temp(fastaFile);
 214                                         m->gobble(fastaFile);
 215
 216                                         if (rdb->save) { rdb->referenceSeqs.push_back(temp); }
 217                                         names.push_back(temp.getName());
 218                                 }
 219                                 fastaFile.close();
 220                         }
 221         #endif
 222
 223                         database->setNumSeqs(names.size());
 224
 225                         //sanity check
 226                         bool okay = phyloTree->ErrorCheck(names);
 227
 228                         if (!okay) { m->control_pressed = true; }
 229
 230                         m->mothurOut("DONE."); m->mothurOutEndLine();
 231                         m->mothurOut("It took " + toString(time(NULL) - start) + " seconds generate search database. "); m->mothurOutEndLine();
 232                 }
 233
 234         }
 235         catch(exception& e) {
 236                 m->errorOut(e, "Classify", "generateDatabaseAndNames");
 237                 exit(1);
 238         }
 239 }
 240 /**************************************************************************************************/
 241 Classify::Classify() {          m = MothurOut::getInstance();   database = NULL;        }
 242 /**************************************************************************************************/
 243
 244 int Classify::readTaxonomy(string file) {
 245         try {
 246
 247                 phyloTree = new PhyloTree();
 248                 string name, taxInfo;
 249
 250                 m->mothurOutEndLine();
 251                 m->mothurOut("Reading in the " + file + " taxonomy...\t");      cout.flush();
 252
 253 #ifdef USE_MPI
 254                 int pid, num, processors;
 255                 vector<unsigned long int> positions;
 256                 int tag = 2001;
 257
 258                 MPI_Status status;
 259                 MPI_File inMPI;
 260                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 261                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
 262
 263                 //char* inFileName = new char[file.length()];
 264                 //memcpy(inFileName, file.c_str(), file.length());
 265
 266                 char inFileName[1024];
 267                 strcpy(inFileName, file.c_str());
 268
 269                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 270                 //delete inFileName;
 271
 272                 if (pid == 0) {
 273                         positions = m->setFilePosEachLine(file, num);
 274
 275                         //send file positions to all processes
 276                         for(int i = 1; i < processors; i++) {
 277                                 MPI_Send(&num, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 278                                 MPI_Send(&positions[0], (num+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
 279                         }
 280                 }else{
 281                         MPI_Recv(&num, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 282                         positions.resize(num+1);
 283                         MPI_Recv(&positions[0], (num+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
 284                 }
 285
 286                 //read file
 287                 for(int i=0;i<num;i++){
 288                         //read next sequence
 289                         int length = positions[i+1] - positions[i];
 290                         char* buf4 = new char[length];
 291
 292                         MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
 293
 294                         string tempBuf = buf4;
 295                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
 296                         delete buf4;
 297
 298                         istringstream iss (tempBuf,istringstream::in);
 299                         iss >> name >> taxInfo;
 300                         taxonomy[name] = taxInfo;
 301                         phyloTree->addSeqToTree(name, taxInfo);
 302                 }
 303
 304                 MPI_File_close(&inMPI);
 305                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 306 #else
 307                 ifstream inTax;
 308                 m->openInputFile(file, inTax);
 309
 310                 //read template seqs and save
 311                 while (!inTax.eof()) {
 312                         inTax >> name >> taxInfo;
 313
 314                         taxonomy[name] = taxInfo;
 315
 316                         phyloTree->addSeqToTree(name, taxInfo);
 317
 318                         m->gobble(inTax);
 319                 }
 320                 inTax.close();
 321 #endif
 322
 323                 phyloTree->assignHeirarchyIDs(0);
 324
 325                 phyloTree->setUp(file);
 326
 327                 m->mothurOut("DONE.");
 328                 m->mothurOutEndLine();  cout.flush();
 329
 330                 return phyloTree->getNumSeqs();
 331
 332         }
 333         catch(exception& e) {
 334                 m->errorOut(e, "Classify", "readTaxonomy");
 335                 exit(1);
 336         }
 337 }
 338 /**************************************************************************************************/
 339
 340 vector<string> Classify::parseTax(string tax) {
 341         try {
 342                 vector<string> taxons;
 343
 344                 tax = tax.substr(0, tax.length()-1);  //get rid of last ';'
 345
 346                 //parse taxonomy
 347                 string individual;
 348                 while (tax.find_first_of(';') != -1) {
 349                         individual = tax.substr(0,tax.find_first_of(';'));
 350                         tax = tax.substr(tax.find_first_of(';')+1, tax.length());
 351                         taxons.push_back(individual);
 352
 353                 }
 354                 //get last one
 355                 taxons.push_back(tax);
 356
 357                 return taxons;
 358         }
 359         catch(exception& e) {
 360                 m->errorOut(e, "Classify", "parseTax");
 361                 exit(1);
 362         }
 363 }
 364 /**************************************************************************************************/
 365