classifyseqscommand.cpp

   1 /*
   2  *  classifyseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 11/2/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "classifyseqscommand.h"
  11 #include "sequence.hpp"
  12 #include "bayesian.h"
  13 #include "phylotree.h"
  14 #include "phylosummary.h"
  15 #include "knn.h"
  16
  17 //**********************************************************************************************************************
  18
  19 ClassifySeqsCommand::ClassifySeqsCommand(string option)  {
  20         try {
  21                 abort = false;
  22
  23                 //allow user to run help
  24                 if(option == "help") { help(); abort = true; }
  25
  26                 else {
  27
  28                         //valid paramters for this command
  29                         string AlignArray[] =  {"template","fasta","name","group","search","ksize","method","processors","taxonomy","match","mismatch","gapopen","gapextend","numwanted","cutoff","probs","iters", "outputdir","inputdir"};
  30                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
  31
  32                         OptionParser parser(option);
  33                         map<string, string> parameters = parser.getParameters();
  34
  35                         ValidParameters validParameter;
  36                         map<string, string>::iterator it;
  37
  38                         //check to make sure all parameters are valid for command
  39                         for (it = parameters.begin(); it != parameters.end(); it++) {
  40                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  41                         }
  42
  43                         //if the user changes the output directory command factory will send this info to us in the output parameter
  44                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = "";         }
  45
  46                         //if the user changes the input directory command factory will send this info to us in the output parameter
  47                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  48                         if (inputDir == "not found"){   inputDir = "";          }
  49                         else {
  50                                 string path;
  51                                 it = parameters.find("template");
  52                                 //user has given a template file
  53                                 if(it != parameters.end()){
  54                                         path = hasPath(it->second);
  55                                         //if the user has not given a path then, add inputdir. else leave path alone.
  56                                         if (path == "") {       parameters["template"] = inputDir + it->second;         }
  57                                 }
  58
  59                                 it = parameters.find("taxonomy");
  60                                 //user has given a template file
  61                                 if(it != parameters.end()){
  62                                         path = hasPath(it->second);
  63                                         //if the user has not given a path then, add inputdir. else leave path alone.
  64                                         if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
  65                                 }
  66
  67                                 it = parameters.find("group");
  68                                 //user has given a template file
  69                                 if(it != parameters.end()){
  70                                         path = hasPath(it->second);
  71                                         //if the user has not given a path then, add inputdir. else leave path alone.
  72                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
  73                                 }
  74                         }
  75
  76                         //check for required parameters
  77                         templateFileName = validParameter.validFile(parameters, "template", true);
  78                         if (templateFileName == "not found") {
  79                                 m->mothurOut("template is a required parameter for the classify.seqs command.");
  80                                 m->mothurOutEndLine();
  81                                 abort = true;
  82                         }
  83                         else if (templateFileName == "not open") { abort = true; }
  84
  85                         groupfile = validParameter.validFile(parameters, "group", false);
  86                         if (groupfile == "not found") { groupfile = "";  }
  87                         else {
  88                                 splitAtDash(groupfile, groupfileNames);
  89
  90                                 //go through files and make sure they are good, if not, then disregard them
  91                                 for (int i = 0; i < groupfileNames.size(); i++) {
  92                                         if (inputDir != "") {
  93                                                 string path = hasPath(groupfileNames[i]);
  94                                                 //if the user has not given a path then, add inputdir. else leave path alone.
  95                                                 if (path == "") {       groupfileNames[i] = inputDir + groupfileNames[i];               }
  96                                         }
  97                                         int ableToOpen;
  98
  99                                         #ifdef USE_MPI
 100                                                 int pid;
 101                                                 MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
 102                                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 103
 104                                                 if (pid == 0) {
 105                                         #endif
 106
 107                                         ifstream in;
 108                                         ableToOpen = openInputFile(groupfileNames[i], in);
 109                                         in.close();
 110
 111                                         #ifdef USE_MPI
 112                                                         for (int j = 1; j < processors; j++) {
 113                                                                 MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
 114                                                         }
 115                                                 }else{
 116                                                         MPI_Status status;
 117                                                         MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
 118                                                 }
 119
 120                                         #endif
 121                                         if (ableToOpen == 1) {  m->mothurOut("Unable to match group file with fasta file."); m->mothurOutEndLine(); abort = true;       }
 122
 123                                 }
 124                         }
 125
 126                         fastaFileName = validParameter.validFile(parameters, "fasta", false);
 127                         if (fastaFileName == "not found") { m->mothurOut("fasta is a required parameter for the classify.seqs command."); m->mothurOutEndLine(); abort = true;  }
 128                         else {
 129                                 splitAtDash(fastaFileName, fastaFileNames);
 130
 131                                 //go through files and make sure they are good, if not, then disregard them
 132                                 for (int i = 0; i < fastaFileNames.size(); i++) {
 133                                         if (inputDir != "") {
 134                                                 string path = hasPath(fastaFileNames[i]);
 135                                                 //if the user has not given a path then, add inputdir. else leave path alone.
 136                                                 if (path == "") {       fastaFileNames[i] = inputDir + fastaFileNames[i];               }
 137                                         }
 138
 139                                         int ableToOpen;
 140
 141                                         #ifdef USE_MPI
 142                                                 int pid;
 143                                                 MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
 144                                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 145
 146                                                 if (pid == 0) {
 147                                         #endif
 148
 149                                         ifstream in;
 150                                         ableToOpen = openInputFile(fastaFileNames[i], in);
 151                                         in.close();
 152
 153                                         #ifdef USE_MPI
 154                                                         for (int j = 1; j < processors; j++) {
 155                                                                 MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
 156                                                         }
 157                                                 }else{
 158                                                         MPI_Status status;
 159                                                         MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
 160                                                 }
 161
 162                                         #endif
 163
 164                                         if (ableToOpen == 1) {
 165                                                 m->mothurOut(fastaFileNames[i] + " will be disregarded."); m->mothurOutEndLine();
 166                                                 //erase from file list
 167                                                 fastaFileNames.erase(fastaFileNames.begin()+i);
 168                                                 i--;
 169                                         }
 170
 171                                 }
 172
 173                                 //make sure there is at least one valid file left
 174                                 if (fastaFileNames.size() == 0) { m->mothurOut("no valid files."); m->mothurOutEndLine(); abort = true; }
 175                         }
 176
 177
 178                         taxonomyFileName = validParameter.validFile(parameters, "taxonomy", true);
 179                         if (taxonomyFileName == "not found") {
 180                                 m->mothurOut("taxonomy is a required parameter for the classify.seqs command.");
 181                                 m->mothurOutEndLine();
 182                                 abort = true;
 183                         }
 184                         else if (taxonomyFileName == "not open") { abort = true; }
 185
 186
 187                         namefile = validParameter.validFile(parameters, "name", false);
 188                         if (namefile == "not found") { namefile = "";  }
 189
 190                         else {
 191                                 splitAtDash(namefile, namefileNames);
 192
 193                                 //go through files and make sure they are good, if not, then disregard them
 194                                 for (int i = 0; i < namefileNames.size(); i++) {
 195                                         if (inputDir != "") {
 196                                                 string path = hasPath(namefileNames[i]);
 197                                                 //if the user has not given a path then, add inputdir. else leave path alone.
 198                                                 if (path == "") {       namefileNames[i] = inputDir + namefileNames[i];         }
 199                                         }
 200                                         int ableToOpen;
 201
 202                                         #ifdef USE_MPI
 203                                                 int pid;
 204                                                 MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
 205                                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 206
 207                                                 if (pid == 0) {
 208                                         #endif
 209
 210                                         ifstream in;
 211                                         ableToOpen = openInputFile(namefileNames[i], in);
 212                                         in.close();
 213
 214                                         #ifdef USE_MPI
 215                                                         for (int j = 1; j < processors; j++) {
 216                                                                 MPI_Send(&ableToOpen, 1, MPI_INT, j, 2001, MPI_COMM_WORLD);
 217                                                         }
 218                                                 }else{
 219                                                         MPI_Status status;
 220                                                         MPI_Recv(&ableToOpen, 1, MPI_INT, 0, 2001, MPI_COMM_WORLD, &status);
 221                                                 }
 222
 223                                         #endif
 224                                         if (ableToOpen == 1) {  m->mothurOut("Unable to match name file with fasta file."); m->mothurOutEndLine(); abort = true;        }
 225
 226                                 }
 227                         }
 228
 229                         if (namefile != "") {
 230                                 if (namefileNames.size() != fastaFileNames.size()) { abort = true; m->mothurOut("If you provide a name file, you must have one for each fasta file."); m->mothurOutEndLine(); }
 231                         }
 232
 233                         if (groupfile != "") {
 234                                 if (groupfileNames.size() != fastaFileNames.size()) { abort = true; m->mothurOut("If you provide a group file, you must have one for each fasta file."); m->mothurOutEndLine(); }
 235                         }
 236
 237                         //check for optional parameter and set defaults
 238                         // ...at some point should added some additional type checking...
 239                         string temp;
 240                         temp = validParameter.validFile(parameters, "ksize", false);            if (temp == "not found"){       temp = "8";                             }
 241                         convert(temp, kmerSize);
 242
 243                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = "1";                             }
 244                         convert(temp, processors);
 245
 246                         search = validParameter.validFile(parameters, "search", false);         if (search == "not found"){     search = "kmer";                }
 247
 248                         method = validParameter.validFile(parameters, "method", false);         if (method == "not found"){     method = "bayesian";    }
 249
 250                         temp = validParameter.validFile(parameters, "match", false);            if (temp == "not found"){       temp = "1.0";                   }
 251                         convert(temp, match);
 252
 253                         temp = validParameter.validFile(parameters, "mismatch", false);         if (temp == "not found"){       temp = "-1.0";                  }
 254                         convert(temp, misMatch);
 255
 256                         temp = validParameter.validFile(parameters, "gapopen", false);          if (temp == "not found"){       temp = "-2.0";                  }
 257                         convert(temp, gapOpen);
 258
 259                         temp = validParameter.validFile(parameters, "gapextend", false);        if (temp == "not found"){       temp = "-1.0";                  }
 260                         convert(temp, gapExtend);
 261
 262                         temp = validParameter.validFile(parameters, "numwanted", false);        if (temp == "not found"){       temp = "10";                    }
 263                         convert(temp, numWanted);
 264
 265                         temp = validParameter.validFile(parameters, "cutoff", false);           if (temp == "not found"){       temp = "0";                             }
 266                         convert(temp, cutoff);
 267
 268                         temp = validParameter.validFile(parameters, "probs", false);            if (temp == "not found"){       temp = "true";                  }
 269                         probs = isTrue(temp);
 270
 271                         temp = validParameter.validFile(parameters, "iters", false);            if (temp == "not found") { temp = "100";                        }
 272                         convert(temp, iters);
 273
 274
 275
 276                         if ((method == "bayesian") && (search != "kmer"))  {
 277                                 m->mothurOut("The bayesian method requires the kmer search." + search + "will be disregarded." ); m->mothurOutEndLine();
 278                                 search = "kmer";
 279                         }
 280                 }
 281
 282         }
 283         catch(exception& e) {
 284                 m->errorOut(e, "ClassifySeqsCommand", "ClassifySeqsCommand");
 285                 exit(1);
 286         }
 287 }
 288
 289 //**********************************************************************************************************************
 290
 291 ClassifySeqsCommand::~ClassifySeqsCommand(){
 292
 293         if (abort == false) {
 294                 for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 295         }
 296 }
 297
 298 //**********************************************************************************************************************
 299
 300 void ClassifySeqsCommand::help(){
 301         try {
 302                 m->mothurOut("The classify.seqs command reads a fasta file containing sequences and creates a .taxonomy file and a .tax.summary file.\n");
 303                 m->mothurOut("The classify.seqs command parameters are template, fasta, name, search, ksize, method, taxonomy, processors, match, mismatch, gapopen, gapextend, numwanted and probs.\n");
 304                 m->mothurOut("The template, fasta and taxonomy parameters are required. You may enter multiple fasta files by separating their names with dashes. ie. fasta=abrecovery.fasta-amzon.fasta \n");
 305                 m->mothurOut("The search parameter allows you to specify the method to find most similar template.  Your options are: suffix, kmer, blast and distance. The default is kmer.\n");
 306                 m->mothurOut("The name parameter allows you add a names file with your fasta file, if you enter multiple fasta files, you must enter matching names files for them.\n");
 307                 m->mothurOut("The group parameter allows you add a group file so you can have the summary totals broken up by group.\n");
 308                 m->mothurOut("The method parameter allows you to specify classification method to use.  Your options are: bayesian and knn. The default is bayesian.\n");
 309                 m->mothurOut("The ksize parameter allows you to specify the kmer size for finding most similar template to candidate.  The default is 8.\n");
 310                 m->mothurOut("The processors parameter allows you to specify the number of processors to use. The default is 1.\n");
 311                 #ifdef USE_MPI
 312                 m->mothurOut("When using MPI, the processors parameter is set to the number of MPI processes running. \n");
 313                 #endif
 314                 m->mothurOut("The match parameter allows you to specify the bonus for having the same base. The default is 1.0.\n");
 315                 m->mothurOut("The mistmatch parameter allows you to specify the penalty for having different bases.  The default is -1.0.\n");
 316                 m->mothurOut("The gapopen parameter allows you to specify the penalty for opening a gap in an alignment. The default is -2.0.\n");
 317                 m->mothurOut("The gapextend parameter allows you to specify the penalty for extending a gap in an alignment.  The default is -1.0.\n");
 318                 m->mothurOut("The numwanted parameter allows you to specify the number of sequence matches you want with the knn method.  The default is 10.\n");
 319                 m->mothurOut("The cutoff parameter allows you to specify a bootstrap confidence threshold for your taxonomy.  The default is 0.\n");
 320                 m->mothurOut("The probs parameter shut off the bootstrapping results for the bayesian method. The default is true, meaning you want the bootstrapping to be run.\n");
 321                 m->mothurOut("The iters parameter allows you to specify how many iterations to do when calculating the bootstrap confidence score for your taxonomy with the bayesian method.  The default is 100.\n");
 322                 m->mothurOut("The classify.seqs command should be in the following format: \n");
 323                 m->mothurOut("classify.seqs(template=yourTemplateFile, fasta=yourFastaFile, method=yourClassificationMethod, search=yourSearchmethod, ksize=yourKmerSize, taxonomy=yourTaxonomyFile, processors=yourProcessors) \n");
 324                 m->mothurOut("Example classify.seqs(fasta=amazon.fasta, template=core.filtered, method=knn, search=gotoh, ksize=8, processors=2)\n");
 325                 m->mothurOut("The .taxonomy file consists of 2 columns: 1 = your sequence name, 2 = the taxonomy for your sequence. \n");
 326                 m->mothurOut("The .tax.summary is a summary of the different taxonomies represented in your fasta file. \n");
 327                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFastaFile).\n\n");
 328         }
 329         catch(exception& e) {
 330                 m->errorOut(e, "ClassifySeqsCommand", "help");
 331                 exit(1);
 332         }
 333 }
 334
 335
 336 //**********************************************************************************************************************
 337
 338 int ClassifySeqsCommand::execute(){
 339         try {
 340                 if (abort == true) {    return 0;       }
 341
 342                 if(method == "bayesian"){       classify = new Bayesian(taxonomyFileName, templateFileName, search, kmerSize, cutoff, iters);           }
 343                 else if(method == "knn"){       classify = new Knn(taxonomyFileName, templateFileName, search, kmerSize, gapOpen, gapExtend, match, misMatch, numWanted);                               }
 344                 else {
 345                         m->mothurOut(search + " is not a valid method option. I will run the command using bayesian.");
 346                         m->mothurOutEndLine();
 347                         classify = new Bayesian(taxonomyFileName, templateFileName, search, kmerSize, cutoff, iters);
 348                 }
 349
 350                 if (m->control_pressed) { delete classify; return 0; }
 351
 352                 vector<string> outputNames;
 353
 354                 for (int s = 0; s < fastaFileNames.size(); s++) {
 355
 356                         m->mothurOut("Classifying sequences from " + fastaFileNames[s] + " ..." ); m->mothurOutEndLine();
 357
 358                         if (outputDir == "") { outputDir += hasPath(fastaFileNames[s]); }
 359                         string newTaxonomyFile = outputDir + getRootName(getSimpleName(fastaFileNames[s])) + getRootName(getSimpleName(taxonomyFileName)) + "taxonomy";
 360                         string tempTaxonomyFile = outputDir + getRootName(getSimpleName(fastaFileNames[s])) + "taxonomy.temp";
 361                         string taxSummary = outputDir + getRootName(getSimpleName(fastaFileNames[s])) + getRootName(getSimpleName(taxonomyFileName)) + "tax.summary";
 362
 363                         outputNames.push_back(newTaxonomyFile);
 364                         outputNames.push_back(taxSummary);
 365
 366                         int start = time(NULL);
 367                         int numFastaSeqs = 0;
 368                         for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 369
 370 #ifdef USE_MPI
 371                                 int pid, end, numSeqsPerProcessor;
 372                                 int tag = 2001;
 373                                 vector<long> MPIPos;
 374
 375                                 MPI_Status status;
 376                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 377                                 MPI_Comm_size(MPI_COMM_WORLD, &processors);
 378
 379                                 MPI_File inMPI;
 380                                 MPI_File outMPINewTax;
 381                                 MPI_File outMPITempTax;
 382
 383                                 int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
 384                                 int inMode=MPI_MODE_RDONLY;
 385
 386                                 //char* outNewTax = new char[newTaxonomyFile.length()];
 387                                 //memcpy(outNewTax, newTaxonomyFile.c_str(), newTaxonomyFile.length());
 388
 389                                 char outNewTax[1024];
 390                                 strcpy(outNewTax, newTaxonomyFile.c_str());
 391
 392                                 //char* outTempTax = new char[tempTaxonomyFile.length()];
 393                                 //memcpy(outTempTax, tempTaxonomyFile.c_str(), tempTaxonomyFile.length());
 394
 395                                 char outTempTax[1024];
 396                                 strcpy(outTempTax, tempTaxonomyFile.c_str());
 397
 398                                 //char* inFileName = new char[fastaFileNames[s].length()];
 399                                 //memcpy(inFileName, fastaFileNames[s].c_str(), fastaFileNames[s].length());
 400
 401                                 char inFileName[1024];
 402                                 strcpy(inFileName, fastaFileNames[s].c_str());
 403
 404                                 MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 405                                 MPI_File_open(MPI_COMM_WORLD, outNewTax, outMode, MPI_INFO_NULL, &outMPINewTax);
 406                                 MPI_File_open(MPI_COMM_WORLD, outTempTax, outMode, MPI_INFO_NULL, &outMPITempTax);
 407
 408                                 //delete outNewTax;
 409                                 //delete outTempTax;
 410                                 //delete inFileName;
 411
 412                                 if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPINewTax);   MPI_File_close(&outMPITempTax);  delete classify; return 0;  }
 413
 414                                 if(namefile != "") {  MPIReadNamesFile(namefileNames[s]);  }
 415
 416                                 if (pid == 0) { //you are the root process
 417
 418                                         MPIPos = setFilePosFasta(fastaFileNames[s], numFastaSeqs); //fills MPIPos, returns numSeqs
 419
 420                                         //send file positions to all processes
 421                                         MPI_Bcast(&numFastaSeqs, 1, MPI_INT, 0, MPI_COMM_WORLD);  //send numSeqs
 422                                         MPI_Bcast(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, MPI_COMM_WORLD); //send file pos
 423
 424                                         //figure out how many sequences you have to align
 425                                         numSeqsPerProcessor = numFastaSeqs / processors;
 426                                         int startIndex =  pid * numSeqsPerProcessor;
 427                                         if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
 428
 429
 430                                         //align your part
 431                                         driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPINewTax, outMPITempTax, MPIPos);
 432
 433                                         if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPINewTax);   MPI_File_close(&outMPITempTax);  for (int i = 0; i < outputNames.size(); i++) {    remove(outputNames[i].c_str()); } delete classify; return 0;  }
 434
 435                                         for (int i = 1; i < processors; i++) {
 436                                                 int done;
 437                                                 MPI_Recv(&done, 1, MPI_INT, i, tag, MPI_COMM_WORLD, &status);
 438                                         }
 439                                 }else{ //you are a child process
 440                                         MPI_Bcast(&numFastaSeqs, 1, MPI_INT, 0, MPI_COMM_WORLD); //get numSeqs
 441                                         MPIPos.resize(numFastaSeqs+1);
 442                                         MPI_Bcast(&MPIPos[0], (numFastaSeqs+1), MPI_LONG, 0, MPI_COMM_WORLD); //get file positions
 443
 444                                         //figure out how many sequences you have to align
 445                                         numSeqsPerProcessor = numFastaSeqs / processors;
 446                                         int startIndex =  pid * numSeqsPerProcessor;
 447                                         if(pid == (processors - 1)){    numSeqsPerProcessor = numFastaSeqs - pid * numSeqsPerProcessor;         }
 448
 449
 450                                         //align your part
 451                                         driverMPI(startIndex, numSeqsPerProcessor, inMPI, outMPINewTax, outMPITempTax, MPIPos);
 452
 453                                         if (m->control_pressed) {  MPI_File_close(&inMPI);  MPI_File_close(&outMPINewTax);   MPI_File_close(&outMPITempTax);  delete classify; return 0;  }
 454
 455                                         int done = 0;
 456                                         MPI_Send(&done, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
 457                                 }
 458
 459                                 //close files
 460                                 MPI_File_close(&inMPI);
 461                                 MPI_File_close(&outMPINewTax);
 462                                 MPI_File_close(&outMPITempTax);
 463
 464 #else
 465                         //read namefile
 466                         if(namefile != "") {
 467                                 nameMap.clear(); //remove old names
 468
 469                                 ifstream inNames;
 470                                 openInputFile(namefileNames[s], inNames);
 471
 472                                 string firstCol, secondCol;
 473                                 while(!inNames.eof()) {
 474                                         inNames >> firstCol >> secondCol; gobble(inNames);
 475                                         nameMap[firstCol] = getNumNames(secondCol);  //ex. seq1 seq1,seq3,seq5 -> seq1 = 3.
 476                                 }
 477                                 inNames.close();
 478                         }
 479
 480         #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 481                         if(processors == 1){
 482                                 ifstream inFASTA;
 483                                 openInputFile(fastaFileNames[s], inFASTA);
 484                                 numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 485                                 inFASTA.close();
 486
 487                                 lines.push_back(new linePair(0, numFastaSeqs));
 488
 489                                 driver(lines[0], newTaxonomyFile, tempTaxonomyFile, fastaFileNames[s]);
 490                         }
 491                         else{
 492                                 vector<int> positions;
 493                                 processIDS.resize(0);
 494
 495                                 ifstream inFASTA;
 496                                 openInputFile(fastaFileNames[s], inFASTA);
 497
 498                                 string input;
 499                                 while(!inFASTA.eof()){
 500                                         input = getline(inFASTA);
 501                                         if (input.length() != 0) {
 502                                                 if(input[0] == '>'){    int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1);       }
 503                                         }
 504                                 }
 505                                 inFASTA.close();
 506
 507                                 numFastaSeqs = positions.size();
 508
 509                                 int numSeqsPerProcessor = numFastaSeqs / processors;
 510
 511                                 for (int i = 0; i < processors; i++) {
 512                                         int startPos = positions[ i * numSeqsPerProcessor ];
 513                                         if(i == processors - 1){
 514                                                 numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
 515                                         }
 516                                         lines.push_back(new linePair(startPos, numSeqsPerProcessor));
 517                                 }
 518                                 createProcesses(newTaxonomyFile, tempTaxonomyFile, fastaFileNames[s]);
 519
 520                                 rename((newTaxonomyFile + toString(processIDS[0]) + ".temp").c_str(), newTaxonomyFile.c_str());
 521                                 rename((tempTaxonomyFile + toString(processIDS[0]) + ".temp").c_str(), tempTaxonomyFile.c_str());
 522
 523                                 for(int i=1;i<processors;i++){
 524                                         appendTaxFiles((newTaxonomyFile + toString(processIDS[i]) + ".temp"), newTaxonomyFile);
 525                                         appendTaxFiles((tempTaxonomyFile + toString(processIDS[i]) + ".temp"), tempTaxonomyFile);
 526                                         remove((newTaxonomyFile + toString(processIDS[i]) + ".temp").c_str());
 527                                         remove((tempTaxonomyFile + toString(processIDS[i]) + ".temp").c_str());
 528                                 }
 529
 530                         }
 531         #else
 532                         ifstream inFASTA;
 533                         openInputFile(fastaFileNames[s], inFASTA);
 534                         numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 535                         inFASTA.close();
 536
 537                         lines.push_back(new linePair(0, numFastaSeqs));
 538
 539                         driver(lines[0], newTaxonomyFile, tempTaxonomyFile, fastaFileNames[s]);
 540         #endif
 541 #endif
 542
 543                 #ifdef USE_MPI
 544                         if (pid == 0) {  //this part does not need to be paralellized
 545                 #endif
 546
 547                         m->mothurOutEndLine();
 548                         m->mothurOut("It took " + toString(time(NULL) - start) + " secs to classify " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
 549                         start = time(NULL);
 550
 551                         PhyloSummary taxaSum(taxonomyFileName, groupfileNames[s]);
 552
 553                         if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } delete classify; return 0; }
 554
 555                         if (namefile == "") {  taxaSum.summarize(tempTaxonomyFile);  }
 556                         else {
 557                                 ifstream in;
 558                                 openInputFile(tempTaxonomyFile, in);
 559
 560                                 //read in users taxonomy file and add sequences to tree
 561                                 string name, taxon;
 562                                 while(!in.eof()){
 563                                         in >> name >> taxon; gobble(in);
 564
 565                                         itNames = nameMap.find(name);
 566
 567                                         if (itNames == nameMap.end()) {
 568                                                 m->mothurOut(name + " is not in your name file please correct."); m->mothurOutEndLine(); exit(1);
 569                                         }else{
 570                                                 for (int i = 0; i < itNames->second; i++) {
 571                                                         taxaSum.addSeqToTree(name, taxon);  //add it as many times as there are identical seqs
 572                                                 }
 573                                         }
 574                                 }
 575                                 in.close();
 576                         }
 577                         remove(tempTaxonomyFile.c_str());
 578
 579                         if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } delete classify; return 0; }
 580
 581                         //print summary file
 582                         ofstream outTaxTree;
 583                         openOutputFile(taxSummary, outTaxTree);
 584                         taxaSum.print(outTaxTree);
 585                         outTaxTree.close();
 586
 587                         //output taxonomy with the unclassified bins added
 588                         ifstream inTax;
 589                         openInputFile(newTaxonomyFile, inTax);
 590
 591                         ofstream outTax;
 592                         string unclass = newTaxonomyFile + ".unclass.temp";
 593                         openOutputFile(unclass, outTax);
 594
 595                         //get maxLevel from phylotree so you know how many 'unclassified's to add
 596                         int maxLevel = taxaSum.getMaxLevel();
 597
 598                         //read taxfile - this reading and rewriting is done to preserve the confidence scores.
 599                         string name, taxon;
 600                         while (!inTax.eof()) {
 601                                 if (m->control_pressed) {  for (int i = 0; i < outputNames.size(); i++) {       remove(outputNames[i].c_str()); } remove(unclass.c_str()); delete classify; return 0; }
 602
 603                                 inTax >> name >> taxon; gobble(inTax);
 604
 605                                 string newTax = addUnclassifieds(taxon, maxLevel);
 606
 607                                 outTax << name << '\t' << newTax << endl;
 608                         }
 609                         inTax.close();
 610                         outTax.close();
 611
 612                         remove(newTaxonomyFile.c_str());
 613                         rename(unclass.c_str(), newTaxonomyFile.c_str());
 614
 615                         m->mothurOutEndLine();
 616                         m->mothurOut("It took " + toString(time(NULL) - start) + " secs to create the summary file for  " + toString(numFastaSeqs) + " sequences."); m->mothurOutEndLine(); m->mothurOutEndLine();
 617
 618                         #ifdef USE_MPI
 619                                 }
 620                         #endif
 621
 622                         m->mothurOutEndLine();
 623                         m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 624                         for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 625                         m->mothurOutEndLine();
 626                 }
 627
 628                 delete classify;
 629                 return 0;
 630         }
 631         catch(exception& e) {
 632                 m->errorOut(e, "ClassifySeqsCommand", "execute");
 633                 exit(1);
 634         }
 635 }
 636
 637 /**************************************************************************************************/
 638 string ClassifySeqsCommand::addUnclassifieds(string tax, int maxlevel) {
 639         try{
 640                 string newTax, taxon;
 641                 int level = 0;
 642
 643                 //keep what you have counting the levels
 644                 while (tax.find_first_of(';') != -1) {
 645                         //get taxon
 646                         taxon = tax.substr(0,tax.find_first_of(';'))+';';
 647                         tax = tax.substr(tax.find_first_of(';')+1, tax.length());
 648                         newTax += taxon;
 649                         level++;
 650                 }
 651
 652                 //add "unclassified" until you reach maxLevel
 653                 while (level < maxlevel) {
 654                         newTax += "unclassified;";
 655                         level++;
 656                 }
 657
 658                 return newTax;
 659         }
 660         catch(exception& e) {
 661                 m->errorOut(e, "ClassifySeqsCommand", "addUnclassifieds");
 662                 exit(1);
 663         }
 664 }
 665
 666 /**************************************************************************************************/
 667
 668 void ClassifySeqsCommand::createProcesses(string taxFileName, string tempTaxFile, string filename) {
 669         try {
 670 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 671                 int process = 0;
 672                 //              processIDS.resize(0);
 673
 674                 //loop through and create all the processes you want
 675                 while (process != processors) {
 676                         int pid = fork();
 677
 678                         if (pid > 0) {
 679                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 680                                 process++;
 681                         }else if (pid == 0){
 682                                 driver(lines[process], taxFileName + toString(getpid()) + ".temp", tempTaxFile + toString(getpid()) + ".temp", filename);
 683                                 exit(0);
 684                         }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
 685                 }
 686
 687                 //force parent to wait until all the processes are done
 688                 for (int i=0;i<processors;i++) {
 689                         int temp = processIDS[i];
 690                         wait(&temp);
 691                 }
 692 #endif
 693         }
 694         catch(exception& e) {
 695                 m->errorOut(e, "ClassifySeqsCommand", "createProcesses");
 696                 exit(1);
 697         }
 698 }
 699 /**************************************************************************************************/
 700
 701 void ClassifySeqsCommand::appendTaxFiles(string temp, string filename) {
 702         try{
 703
 704                 ofstream output;
 705                 ifstream input;
 706                 openOutputFileAppend(filename, output);
 707                 openInputFile(temp, input);
 708
 709                 while(char c = input.get()){
 710                         if(input.eof())         {       break;                  }
 711                         else                            {       output << c;    }
 712                 }
 713
 714                 input.close();
 715                 output.close();
 716         }
 717         catch(exception& e) {
 718                 m->errorOut(e, "ClassifySeqsCommand", "appendTaxFiles");
 719                 exit(1);
 720         }
 721 }
 722
 723 //**********************************************************************************************************************
 724
 725 int ClassifySeqsCommand::driver(linePair* line, string taxFName, string tempTFName, string filename){
 726         try {
 727                 ofstream outTax;
 728                 openOutputFile(taxFName, outTax);
 729
 730                 ofstream outTaxSimple;
 731                 openOutputFile(tempTFName, outTaxSimple);
 732
 733                 ifstream inFASTA;
 734                 openInputFile(filename, inFASTA);
 735
 736                 inFASTA.seekg(line->start);
 737
 738                 string taxonomy;
 739
 740                 for(int i=0;i<line->numSeqs;i++){
 741                         if (m->control_pressed) { return 0; }
 742
 743                         Sequence* candidateSeq = new Sequence(inFASTA);
 744
 745                         if (candidateSeq->getName() != "") {
 746                                 taxonomy = classify->getTaxonomy(candidateSeq);
 747
 748                                 if (m->control_pressed) { delete candidateSeq; return 0; }
 749
 750                                 if (taxonomy != "bad seq") {
 751                                         //output confidence scores or not
 752                                         if (probs) {
 753                                                 outTax << candidateSeq->getName() << '\t' << taxonomy << endl;
 754                                         }else{
 755                                                 outTax << candidateSeq->getName() << '\t' << classify->getSimpleTax() << endl;
 756                                         }
 757
 758                                         outTaxSimple << candidateSeq->getName() << '\t' << classify->getSimpleTax() << endl;
 759                                 }
 760                         }
 761                         delete candidateSeq;
 762
 763                         if((i+1) % 100 == 0){
 764                                 m->mothurOut("Classifying sequence " + toString(i+1)); m->mothurOutEndLine();
 765                         }
 766                 }
 767
 768                 inFASTA.close();
 769                 outTax.close();
 770                 outTaxSimple.close();
 771
 772                 return 1;
 773         }
 774         catch(exception& e) {
 775                 m->errorOut(e, "ClassifySeqsCommand", "driver");
 776                 exit(1);
 777         }
 778 }
 779 //**********************************************************************************************************************
 780 #ifdef USE_MPI
 781 int ClassifySeqsCommand::driverMPI(int start, int num, MPI_File& inMPI, MPI_File& newFile, MPI_File& tempFile, vector<long>& MPIPos){
 782         try {
 783                 MPI_Status statusNew;
 784                 MPI_Status statusTemp;
 785                 MPI_Status status;
 786
 787                 int pid;
 788                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 789
 790                 string taxonomy;
 791                 string outputString;
 792
 793                 for(int i=0;i<num;i++){
 794
 795                         if (m->control_pressed) { return 0; }
 796
 797                         //read next sequence
 798                         int length = MPIPos[start+i+1] - MPIPos[start+i];
 799                         char* buf4 = new char[length];
 800                         MPI_File_read_at(inMPI, MPIPos[start+i], buf4, length, MPI_CHAR, &status);
 801
 802                         string tempBuf = buf4;
 803                         if (tempBuf.length() > length) { tempBuf = tempBuf.substr(0, length);  }
 804                         istringstream iss (tempBuf,istringstream::in);
 805                         delete buf4;
 806
 807                         Sequence* candidateSeq = new Sequence(iss);
 808
 809                         if (candidateSeq->getName() != "") {
 810                                 taxonomy = classify->getTaxonomy(candidateSeq);
 811
 812                                 if (taxonomy != "bad seq") {
 813                                         //output confidence scores or not
 814                                         if (probs) {
 815                                                 outputString =  candidateSeq->getName() + "\t" + taxonomy + "\n";
 816                                         }else{
 817                                                 outputString =  candidateSeq->getName() + "\t" + classify->getSimpleTax() + "\n";
 818                                         }
 819
 820                                         int length = outputString.length();
 821                                         char* buf2 = new char[length];
 822                                         memcpy(buf2, outputString.c_str(), length);
 823
 824                                         MPI_File_write_shared(newFile, buf2, length, MPI_CHAR, &statusNew);
 825                                         delete buf2;
 826
 827                                         outputString =  candidateSeq->getName() + "\t" + classify->getSimpleTax() + "\n";
 828                                         length = outputString.length();
 829                                         char* buf = new char[length];
 830                                         memcpy(buf, outputString.c_str(), length);
 831
 832                                         MPI_File_write_shared(tempFile, buf, length, MPI_CHAR, &statusTemp);
 833                                         delete buf;
 834                                 }
 835                         }
 836                         delete candidateSeq;
 837
 838                         if((i+1) % 100 == 0){   cout << "Classifying sequence " << (i+1) << endl;       }
 839                 }
 840
 841                 if(num % 100 != 0){     cout << "Classifying sequence " << (num) << endl;       }
 842
 843
 844                 return 1;
 845         }
 846         catch(exception& e) {
 847                 m->errorOut(e, "ClassifySeqsCommand", "driverMPI");
 848                 exit(1);
 849         }
 850 }
 851
 852 //**********************************************************************************************************************
 853 int ClassifySeqsCommand::MPIReadNamesFile(string nameFilename){
 854         try {
 855
 856                 nameMap.clear(); //remove old names
 857
 858                 MPI_File inMPI;
 859                 MPI_Offset size;
 860                 MPI_Status status;
 861
 862                 //char* inFileName = new char[nameFilename.length()];
 863                 //memcpy(inFileName, nameFilename.c_str(), nameFilename.length());
 864
 865                 char inFileName[1024];
 866                 strcpy(inFileName, nameFilename.c_str());
 867
 868                 MPI_File_open(MPI_COMM_WORLD, inFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &inMPI);
 869                 MPI_File_get_size(inMPI, &size);
 870                 //delete inFileName;
 871
 872                 char* buffer = new char[size];
 873                 MPI_File_read(inMPI, buffer, size, MPI_CHAR, &status);
 874
 875                 string tempBuf = buffer;
 876                 if (tempBuf.length() > size) { tempBuf = tempBuf.substr(0, size);  }
 877                 istringstream iss (tempBuf,istringstream::in);
 878                 delete buffer;
 879
 880                 string firstCol, secondCol;
 881                 while(!iss.eof()) {
 882                         iss >> firstCol >> secondCol; gobble(iss);
 883                         nameMap[firstCol] = getNumNames(secondCol);  //ex. seq1 seq1,seq3,seq5 -> seq1 = 3.
 884                 }
 885
 886                 MPI_File_close(&inMPI);
 887
 888                 return 1;
 889         }
 890         catch(exception& e) {
 891                 m->errorOut(e, "ClassifySeqsCommand", "MPIReadNamesFile");
 892                 exit(1);
 893         }
 894 }
 895 #endif
 896 /**************************************************************************************************/