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<long> 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                 database->setNumSeqs(names.size());
 155
 156                 //sanity check
 157                 bool okay = phyloTree->ErrorCheck(names);
 158
 159                 if (!okay) { m->control_pressed = true; }
 160
 161                 m->mothurOut("DONE."); m->mothurOutEndLine();
 162                 m->mothurOut("It took " + toString(time(NULL) - start) + " seconds generate search database. "); m->mothurOutEndLine();
 163
 164         }
 165         catch(exception& e) {
 166                 m->errorOut(e, "Classify", "generateDatabaseAndNames");
 167                 exit(1);
 168         }
 169 }
 170 /**************************************************************************************************/
 171 Classify::Classify() {          m = MothurOut::getInstance();   database = NULL;        }
 172 /**************************************************************************************************/
 173
 174 int Classify::readTaxonomy(string file) {
 175         try {
 176
 177                 phyloTree = new PhyloTree();
 178                 string name, taxInfo;
 179
 180                 m->mothurOutEndLine();
 181                 m->mothurOut("Reading in the " + file + " taxonomy...\t");      cout.flush();
 182
 183 #ifdef USE_MPI
 184                 int pid, num, processors;
 185                 vector<long> positions;
 186                 int tag = 2001;
 187
 188                 MPI_Status status;
 189                 MPI_File inMPI;
 190                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 191                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
 192
 193                 //char* inFileName = new char[file.length()];
 194                 //memcpy(inFileName, file.c_str(), file.length());
 195
 196                 char inFileName[1024];
 197                 strcpy(inFileName, file.c_str());
 198
 199                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 200                 //delete inFileName;
 201
 202                 if (pid == 0) {
 203                         positions = setFilePosEachLine(file, num);
 204
 205                         //send file positions to all processes
 206                         for(int i = 1; i < processors; i++) {
 207                                 MPI_Send(&num, 1, MPI_INT, i, tag, MPI_COMM_WORLD);
 208                                 MPI_Send(&positions[0], (num+1), MPI_LONG, i, tag, MPI_COMM_WORLD);
 209                         }
 210                 }else{
 211                         MPI_Recv(&num, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 212                         positions.resize(num+1);
 213                         MPI_Recv(&positions[0], (num+1), MPI_LONG, 0, tag, MPI_COMM_WORLD, &status);
 214                 }
 215
 216                 //read file
 217                 for(int i=0;i<num;i++){
 218                         //read next sequence
 219                         int length = positions[i+1] - positions[i];
 220                         char* buf4 = new char[length];
 221
 222                         MPI_File_read_at(inMPI, positions[i], buf4, length, MPI_CHAR, &status);
 223
 224                         string tempBuf = buf4;
 225                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length); }
 226                         delete buf4;
 227
 228                         istringstream iss (tempBuf,istringstream::in);
 229                         iss >> name >> taxInfo;
 230                         taxonomy[name] = taxInfo;
 231                         phyloTree->addSeqToTree(name, taxInfo);
 232                 }
 233
 234                 MPI_File_close(&inMPI);
 235                 MPI_Barrier(MPI_COMM_WORLD); //make everyone wait - just in case
 236 #else
 237                 ifstream inTax;
 238                 openInputFile(file, inTax);
 239
 240                 //read template seqs and save
 241                 while (!inTax.eof()) {
 242                         inTax >> name >> taxInfo;
 243
 244                         taxonomy[name] = taxInfo;
 245
 246                         phyloTree->addSeqToTree(name, taxInfo);
 247
 248                         gobble(inTax);
 249                 }
 250                 inTax.close();
 251 #endif
 252
 253                 phyloTree->assignHeirarchyIDs(0);
 254
 255                 phyloTree->setUp(file);
 256
 257                 m->mothurOut("DONE.");
 258                 m->mothurOutEndLine();  cout.flush();
 259
 260                 return phyloTree->getNumSeqs();
 261
 262         }
 263         catch(exception& e) {
 264                 m->errorOut(e, "Classify", "readTaxonomy");
 265                 exit(1);
 266         }
 267 }
 268 /**************************************************************************************************/
 269
 270 vector<string> Classify::parseTax(string tax) {
 271         try {
 272                 vector<string> taxons;
 273
 274                 tax = tax.substr(0, tax.length()-1);  //get rid of last ';'
 275
 276                 //parse taxonomy
 277                 string individual;
 278                 while (tax.find_first_of(';') != -1) {
 279                         individual = tax.substr(0,tax.find_first_of(';'));
 280                         tax = tax.substr(tax.find_first_of(';')+1, tax.length());
 281                         taxons.push_back(individual);
 282
 283                 }
 284                 //get last one
 285                 taxons.push_back(tax);
 286
 287                 return taxons;
 288         }
 289         catch(exception& e) {
 290                 m->errorOut(e, "Classify", "parseTax");
 291                 exit(1);
 292         }
 293 }
 294 /**************************************************************************************************/
 295