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
  17 /**************************************************************************************************/
  18 void Classify::generateDatabaseAndNames(string tfile, string tempFile, string method, int kmerSize, float gapOpen, float gapExtend, float match, float misMatch)  {
  19         try {
  20                 taxFile = tfile;
  21                 readTaxonomy(taxFile);
  22
  23                 templateFile = tempFile;
  24
  25                 int start = time(NULL);
  26                 int numSeqs = 0;
  27
  28                 m->mothurOut("Generating search database...    "); cout.flush();
  29 #ifdef USE_MPI
  30                         int pid, processors;
  31                         vector<unsigned long int> positions;
  32                         int tag = 2001;
  33
  34                         MPI_Status status;
  35                         MPI_File inMPI;
  36                         MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
  37                         MPI_Comm_size(MPI_COMM_WORLD, &processors);
  38
  39                         //char* inFileName = new char[tempFile.length()];
  40                         //memcpy(inFileName, tempFile.c_str(), tempFile.length());
  41
  42                         char inFileName[1024];
  43                         strcpy(inFileName, tempFile.c_str());
  44
  45                         MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
  46                         //delete inFileName;
  47
  48                         if (pid == 0) { //only one process needs to scan file
  49                                 positions = setFilePosFasta(tempFile, numSeqs); //fills MPIPos, returns numSeqs
  50
  51                                 //send file positions to all processes
  52                                 for(int i = 1; i < processors; i++) {
  53                                         MPI_Send(&numSeqs, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
  54                                         MPI_Send(&positions[0], (numSeqs+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
  55                                 }
  56                         }else{
  57                                 MPI_Recv(&numSeqs, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
  58                                 positions.resize(numSeqs+1);
  59                                 MPI_Recv(&positions[0], (numSeqs+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
  60                         }
  61
  62                         //create database
  63                         if(method == "kmer")                    {       database = new KmerDB(tempFile, kmerSize);                      }
  64                         else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
  65                         else if(method == "blast")              {       database = new BlastDB(gapOpen, gapExtend, match, misMatch);    }
  66                         else if(method == "distance")   {       database = new DistanceDB();    }
  67                         else {
  68                                 m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8."); m->mothurOutEndLine();
  69                                 database = new KmerDB(tempFile, 8);
  70                         }
  71
  72                         //read file
  73                         for(int i=0;i<numSeqs;i++){
  74                                 //read next sequence
  75                                 int length = positions[i+1] - positions[i];
  76                                 char* buf4 = new char[length];
  77                                 MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
  78
  79                                 string tempBuf = buf4;
  80                                 if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
  81                                 delete buf4;
  82                                 istringstream iss (tempBuf,istringstream::in);
  83
  84                                 Sequence temp(iss);
  85                                 if (temp.getName() != "") {
  86                                         names.push_back(temp.getName());
  87                                         database->addSequence(temp);
  88                                 }
  89                         }
  90
  91                         database->generateDB();
  92                         MPI_File_close(&inMPI);
  93                         MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
  94         #else
  95
  96                 //need to know number of template seqs for suffixdb
  97                 if (method == "suffix") {
  98                         ifstream inFASTA;
  99                         openInputFile(tempFile, inFASTA);
 100                         getNumSeqs(inFASTA, numSeqs);
 101                         inFASTA.close();
 102                 }
 103
 104                 bool needToGenerate = true;
 105                 string kmerDBName;
 106                 if(method == "kmer")                    {
 107                         database = new KmerDB(tempFile, kmerSize);
 108
 109                         kmerDBName = tempFile.substr(0,tempFile.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
 110                         ifstream kmerFileTest(kmerDBName.c_str());
 111                         if(kmerFileTest){       needToGenerate = false;         }
 112                 }
 113                 else if(method == "suffix")             {       database = new SuffixDB(numSeqs);                                                               }
 114                 else if(method == "blast")              {       database = new BlastDB(gapOpen, gapExtend, match, misMatch);    }
 115                 else if(method == "distance")   {       database = new DistanceDB();    }
 116                 else {
 117                         m->mothurOut(method + " is not a valid search option. I will run the command using kmer, ksize=8.");
 118                         m->mothurOutEndLine();
 119                         database = new KmerDB(tempFile, 8);
 120                 }
 121
 122                 if (needToGenerate) {
 123                         ifstream fastaFile;
 124                         openInputFile(tempFile, fastaFile);
 125
 126                         while (!fastaFile.eof()) {
 127                                 Sequence temp(fastaFile);
 128                                 gobble(fastaFile);
 129
 130                                 names.push_back(temp.getName());
 131
 132                                 database->addSequence(temp);
 133                         }
 134                         fastaFile.close();
 135
 136                         database->generateDB();
 137
 138                 }else if ((method == "kmer") && (!needToGenerate)) {
 139                         ifstream kmerFileTest(kmerDBName.c_str());
 140                         database->readKmerDB(kmerFileTest);
 141
 142                         ifstream fastaFile;
 143                         openInputFile(tempFile, fastaFile);
 144
 145                         while (!fastaFile.eof()) {
 146                                 Sequence temp(fastaFile);
 147                                 gobble(fastaFile);
 148
 149                                 names.push_back(temp.getName());
 150                         }
 151                         fastaFile.close();
 152                 }
 153 #endif
 154
 155                 database->setNumSeqs(names.size());
 156
 157                 //sanity check
 158                 bool okay = phyloTree->ErrorCheck(names);
 159
 160                 if (!okay) { m->control_pressed = true; }
 161
 162                 m->mothurOut("DONE."); m->mothurOutEndLine();
 163                 m->mothurOut("It took " + toString(time(NULL) - start) + " seconds generate search database. "); m->mothurOutEndLine();
 164
 165         }
 166         catch(exception& e) {
 167                 m->errorOut(e, "Classify", "generateDatabaseAndNames");
 168                 exit(1);
 169         }
 170 }
 171 /**************************************************************************************************/
 172 Classify::Classify() {          m = MothurOut::getInstance();   database = NULL;        }
 173 /**************************************************************************************************/
 174
 175 int Classify::readTaxonomy(string file) {
 176         try {
 177
 178                 phyloTree = new PhyloTree();
 179                 string name, taxInfo;
 180
 181                 m->mothurOutEndLine();
 182                 m->mothurOut("Reading in the " + file + " taxonomy...\t");      cout.flush();
 183
 184 #ifdef USE_MPI
 185                 int pid, num, processors;
 186                 vector<unsigned long int> positions;
 187                 int tag = 2001;
 188
 189                 MPI_Status status;
 190                 MPI_File inMPI;
 191                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 192                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
 193
 194                 //char* inFileName = new char[file.length()];
 195                 //memcpy(inFileName, file.c_str(), file.length());
 196
 197                 char inFileName[1024];
 198                 strcpy(inFileName, file.c_str());
 199
 200                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 201                 //delete inFileName;
 202
 203                 if (pid == 0) {
 204                         positions = setFilePosEachLine(file, num);
 205
 206                         //send file positions to all processes
 207                         for(int i = 1; i < processors; i++) {
 208                                 MPI_Send(&num, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 209                                 MPI_Send(&positions[0], (num+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
 210                         }
 211                 }else{
 212                         MPI_Recv(&num, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 213                         positions.resize(num+1);
 214                         MPI_Recv(&positions[0], (num+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
 215                 }
 216
 217                 //read file
 218                 for(int i=0;i<num;i++){
 219                         //read next sequence
 220                         int length = positions[i+1] - positions[i];
 221                         char* buf4 = new char[length];
 222
 223                         MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
 224
 225                         string tempBuf = buf4;
 226                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
 227                         delete buf4;
 228
 229                         istringstream iss (tempBuf,istringstream::in);
 230                         iss >> name >> taxInfo;
 231                         taxonomy[name] = taxInfo;
 232                         phyloTree->addSeqToTree(name, taxInfo);
 233                 }
 234
 235                 MPI_File_close(&inMPI);
 236                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 237 #else
 238                 ifstream inTax;
 239                 openInputFile(file, inTax);
 240
 241                 //read template seqs and save
 242                 while (!inTax.eof()) {
 243                         inTax >> name >> taxInfo;
 244
 245                         taxonomy[name] = taxInfo;
 246
 247                         phyloTree->addSeqToTree(name, taxInfo);
 248
 249                         gobble(inTax);
 250                 }
 251                 inTax.close();
 252 #endif
 253
 254                 phyloTree->assignHeirarchyIDs(0);
 255
 256                 phyloTree->setUp(file);
 257
 258                 m->mothurOut("DONE.");
 259                 m->mothurOutEndLine();  cout.flush();
 260
 261                 return phyloTree->getNumSeqs();
 262
 263         }
 264         catch(exception& e) {
 265                 m->errorOut(e, "Classify", "readTaxonomy");
 266                 exit(1);
 267         }
 268 }
 269 /**************************************************************************************************/
 270
 271 vector<string> Classify::parseTax(string tax) {
 272         try {
 273                 vector<string> taxons;
 274
 275                 tax = tax.substr(0, tax.length()-1);  //get rid of last ';'
 276
 277                 //parse taxonomy
 278                 string individual;
 279                 while (tax.find_first_of(';') != -1) {
 280                         individual = tax.substr(0,tax.find_first_of(';'));
 281                         tax = tax.substr(tax.find_first_of(';')+1, tax.length());
 282                         taxons.push_back(individual);
 283
 284                 }
 285                 //get last one
 286                 taxons.push_back(tax);
 287
 288                 return taxons;
 289         }
 290         catch(exception& e) {
 291                 m->errorOut(e, "Classify", "parseTax");
 292                 exit(1);
 293         }
 294 }
 295 /**************************************************************************************************/
 296