filterseqscommand.cpp

   1 /*
   2  *  filterseqscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Thomas Ryabin on 5/4/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "filterseqscommand.h"
  11 #include "sequence.hpp"
  12
  13 /**************************************************************************************/
  14
  15 FilterSeqsCommand::FilterSeqsCommand(string option)  {
  16         try {
  17                 abort = false;
  18                 filterFileName = "";
  19
  20                 //allow user to run help
  21                 if(option == "help") { help(); abort = true; }
  22
  23                 else {
  24                         //valid paramters for this command
  25                         string Array[] =  {"fasta", "trump", "soft", "hard", "vertical", "outputdir","inputdir", "processors"};
  26                         vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  27
  28                         OptionParser parser(option);
  29                         map<string,string> parameters = parser.getParameters();
  30
  31                         ValidParameters validParameter;
  32                         map<string,string>::iterator it;
  33
  34                         //check to make sure all parameters are valid for command
  35                         for (it = parameters.begin(); it != parameters.end(); it++) {
  36                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
  37                         }
  38
  39                         //if the user changes the input directory command factory will send this info to us in the output parameter
  40                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
  41                         if (inputDir == "not found"){   inputDir = "";          }
  42                         else {
  43                                 string path;
  44                                 it = parameters.find("fasta");
  45                                 //user has given a template file
  46                                 if(it != parameters.end()){
  47                                         path = hasPath(it->second);
  48                                         //if the user has not given a path then, add inputdir. else leave path alone.
  49                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
  50                                 }
  51
  52                                 it = parameters.find("hard");
  53                                 //user has given a template file
  54                                 if(it != parameters.end()){
  55                                         path = hasPath(it->second);
  56                                         //if the user has not given a path then, add inputdir. else leave path alone.
  57                                         if (path == "") {       parameters["hard"] = inputDir + it->second;             }
  58                                 }
  59                         }
  60
  61                         //check for required parameters
  62                         fasta = validParameter.validFile(parameters, "fasta", false);
  63                         if (fasta == "not found") { m->mothurOut("fasta is a required parameter for the filter.seqs command."); m->mothurOutEndLine(); abort = true;  }
  64                         else {
  65                                 splitAtDash(fasta, fastafileNames);
  66
  67                                 //go through files and make sure they are good, if not, then disregard them
  68                                 for (int i = 0; i < fastafileNames.size(); i++) {
  69                                         if (inputDir != "") {
  70                                                 string path = hasPath(fastafileNames[i]);
  71                                                 //if the user has not given a path then, add inputdir. else leave path alone.
  72                                                 if (path == "") {       fastafileNames[i] = inputDir + fastafileNames[i];               }
  73                                         }
  74
  75                                         int ableToOpen;
  76                                         ifstream in;
  77                                         ableToOpen = openInputFile(fastafileNames[i], in);
  78                                         if (ableToOpen == 1) {
  79                                                 m->mothurOut(fastafileNames[i] + " will be disregarded."); m->mothurOutEndLine();
  80                                                 //erase from file list
  81                                                 fastafileNames.erase(fastafileNames.begin()+i);
  82                                                 i--;
  83                                         }else{
  84                                                 string simpleName = getSimpleName(fastafileNames[i]);
  85                                                 filterFileName += simpleName.substr(0, simpleName.find_first_of('.'));
  86                                         }
  87                                         in.close();
  88                                 }
  89
  90                                 //make sure there is at least one valid file left
  91                                 if (fastafileNames.size() == 0) { m->mothurOut("no valid files."); m->mothurOutEndLine(); abort = true; }
  92                         }
  93
  94                         if (!abort) {
  95                                 //if the user changes the output directory command factory will send this info to us in the output parameter
  96                                 outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
  97                                         outputDir = "";
  98                                         outputDir += hasPath(fastafileNames[0]); //if user entered a file with a path then preserve it
  99                                 }
 100                         }
 101                         //check for optional parameter and set defaults
 102                         // ...at some point should added some additional type checking...
 103
 104                         string temp;
 105                         temp = validParameter.validFile(parameters, "trump", false);                    if (temp == "not found") { temp = "*"; }
 106                         trump = temp[0];
 107
 108                         temp = validParameter.validFile(parameters, "soft", false);                             if (temp == "not found") { soft = 0; }
 109                         else {  soft = (float)atoi(temp.c_str()) / 100.0;  }
 110
 111                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = "1";                             }
 112                         convert(temp, processors);
 113
 114                         hard = validParameter.validFile(parameters, "hard", true);                              if (hard == "not found") { hard = ""; }
 115                         else if (hard == "not open") { abort = true; }
 116
 117                         vertical = validParameter.validFile(parameters, "vertical", false);             if (vertical == "not found") { vertical = "T"; }
 118
 119                         numSeqs = 0;
 120
 121                 }
 122
 123         }
 124         catch(exception& e) {
 125                 m->errorOut(e, "FilterSeqsCommand", "FilterSeqsCommand");
 126                 exit(1);
 127         }
 128 }
 129
 130 //**********************************************************************************************************************
 131
 132 void FilterSeqsCommand::help(){
 133         try {
 134
 135                 m->mothurOut("The filter.seqs command reads a file containing sequences and creates a .filter and .filter.fasta file.\n");
 136                 m->mothurOut("The filter.seqs command parameters are fasta, trump, soft, hard and vertical. \n");
 137                 m->mothurOut("The fasta parameter is required. You may enter several fasta files to build the filter from and filter, by separating their names with -'s.\n");
 138                 m->mothurOut("For example: fasta=abrecovery.fasta-amazon.fasta \n");
 139                 m->mothurOut("The trump parameter .... The default is ...\n");
 140                 m->mothurOut("The soft parameter .... The default is ....\n");
 141                 m->mothurOut("The hard parameter .... The default is ....\n");
 142                 m->mothurOut("The vertical parameter .... The default is T.\n");
 143                 m->mothurOut("The filter.seqs command should be in the following format: \n");
 144                 m->mothurOut("filter.seqs(fasta=yourFastaFile, trump=yourTrump, soft=yourSoft, hard=yourHard, vertical=yourVertical) \n");
 145                 m->mothurOut("Example filter.seqs(fasta=abrecovery.fasta, trump=..., soft=..., hard=..., vertical=T).\n");
 146                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n\n");
 147
 148         }
 149         catch(exception& e) {
 150                 m->errorOut(e, "FilterSeqsCommand", "help");
 151                 exit(1);
 152         }
 153 }
 154
 155 /**************************************************************************************/
 156
 157 int FilterSeqsCommand::execute() {
 158         try {
 159
 160                 if (abort == true) { return 0; }
 161
 162                 ifstream inFASTA;
 163                 openInputFile(fastafileNames[0], inFASTA);
 164
 165                 Sequence testSeq(inFASTA);
 166                 alignmentLength = testSeq.getAlignLength();
 167                 inFASTA.close();
 168
 169                 ////////////create filter/////////////////
 170
 171                 filter = createFilter();
 172
 173                 if (m->control_pressed) { return 0; }
 174
 175                 ofstream outFilter;
 176
 177                 string filterFile = outputDir + filterFileName + ".filter";
 178                 openOutputFile(filterFile, outFilter);
 179                 outFilter << filter << endl;
 180                 outFilter.close();
 181                 outputNames.push_back(filterFile);
 182
 183
 184                 ////////////run filter/////////////////
 185
 186                 filterSequences();
 187
 188                 int filteredLength = 0;
 189                 for(int i=0;i<alignmentLength;i++){
 190                         if(filter[i] == '1'){   filteredLength++;       }
 191                 }
 192
 193                 if (m->control_pressed) {  for(int i = 0; i < outputNames.size(); i++) { remove(outputNames[i].c_str()); }  return 0; }
 194
 195
 196                 m->mothurOutEndLine();
 197                 m->mothurOut("Length of filtered alignment: " + toString(filteredLength)); m->mothurOutEndLine();
 198                 m->mothurOut("Number of columns removed: " + toString((alignmentLength-filteredLength))); m->mothurOutEndLine();
 199                 m->mothurOut("Length of the original alignment: " + toString(alignmentLength)); m->mothurOutEndLine();
 200                 m->mothurOut("Number of sequences used to construct filter: " + toString(numSeqs)); m->mothurOutEndLine();
 201
 202
 203                 m->mothurOutEndLine();
 204                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 205                 for(int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();      }
 206                 m->mothurOutEndLine();
 207
 208                 return 0;
 209
 210         }
 211         catch(exception& e) {
 212                 m->errorOut(e, "FilterSeqsCommand", "execute");
 213                 exit(1);
 214         }
 215 }
 216 /**************************************************************************************/
 217 int FilterSeqsCommand::filterSequences() {
 218         try {
 219
 220                 for (int s = 0; s < fastafileNames.size(); s++) {
 221
 222                                 for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 223
 224                                 string filteredFasta = outputDir + getRootName(getSimpleName(fastafileNames[s])) + "filter.fasta";
 225 #ifdef USE_MPI
 226                                 int pid, start, end;
 227                                 int tag = 2001;
 228
 229                                 MPI_Status status;
 230                                 MPI_Comm_size(MPI_COMM_WORLD, &processors); //set processors to the number of mpi processes running
 231                                 MPI_Comm_rank(MPI_COMM_WORLD, &pid); //find out who we are
 232
 233                                 MPI_File outMPI;
 234                                 MPI_File inMPI;
 235                                 int outMode=MPI_MODE_CREATE|MPI_MODE_WRONLY;
 236                                 int inMode=MPI_MODE_RDONLY;
 237
 238                                 char outFilename[filteredFasta.length()];
 239                                 strcpy(outFilename, filteredFasta.c_str());
 240
 241                                 char inFileName[fastafileNames[s].length()];
 242                                 strcpy(inFileName, fastafileNames[s].c_str());
 243
 244                                 MPI_File_open(MPI_COMM_WORLD, inFileName, inMode, MPI_INFO_NULL, &inMPI);  //comm, filename, mode, info, filepointer
 245                                 MPI_File_open(MPI_COMM_WORLD, outFilename, outMode, MPI_INFO_NULL, &outMPI);
 246
 247                                 if (pid == 0) { //you are the root process
 248
 249                                         setLines(fastafileNames[s]);
 250
 251                                         for (int j = 0; j < lines.size(); j++) { //each process
 252                                                 if (j != 0) { //don't send to yourself
 253                                                         MPI_Send(&lines[j]->start, 1, MPI_INT, j, tag, MPI_COMM_WORLD); //start position in file
 254                                                         MPI_Send(&bufferSizes[j], 1, MPI_INT, j, tag, MPI_COMM_WORLD); //how bytes for the read
 255                                                 }
 256                                         }
 257
 258                                         //read your peice of file
 259                                         char buf[bufferSizes[0]];
 260                                         MPI_File_read_at(inMPI, lines[0]->start, buf, bufferSizes[0], MPI_CHAR, &status);
 261                                         istringstream iss (buf,istringstream::in);
 262
 263                                         //do your part
 264                                         driverMPIRun(iss, outMPI);
 265
 266                                         //wait on chidren
 267                                         for(int i = 1; i < processors; i++) {
 268                                                 char buf[4];
 269                                                 MPI_Recv(buf, 4, MPI_CHAR, i, tag, MPI_COMM_WORLD, &status);
 270                                         }
 271
 272                                 }else { //you are a child process
 273                                         //receive your section of file
 274                                         int startPos, numLines, bufferSize;
 275                                         MPI_Recv(&startPos, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 276                                         MPI_Recv(&bufferSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 277
 278                                         //read your peice of file
 279                                         char buf2[bufferSize];
 280                                         MPI_File_read_at(inMPI, startPos, buf2, bufferSize, MPI_CHAR, &status);
 281                                         istringstream iss (buf2,istringstream::in);
 282
 283                                         //do your part
 284                                         driverMPIRun(iss, outMPI);
 285
 286                                         char buf[4];
 287                                         strcpy(buf, "done");
 288
 289                                         //tell parent you are done.
 290                                         MPI_Send(buf, 4, MPI_CHAR, 0, tag, MPI_COMM_WORLD);
 291                                 }
 292
 293                                 MPI_File_close(&outMPI);
 294                                 MPI_File_close(&inMPI);
 295
 296 #else
 297                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 298                                 if(processors == 1){
 299                                         ifstream inFASTA;
 300                                         openInputFile(fastafileNames[s], inFASTA);
 301                                         int numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 302                                         inFASTA.close();
 303
 304                                         lines.push_back(new linePair(0, numFastaSeqs));
 305
 306                                         driverRunFilter(filter, filteredFasta, fastafileNames[s], lines[0]);
 307                                 }else{
 308                                         setLines(fastafileNames[s]);
 309                                         createProcessesRunFilter(filter, fastafileNames[s]);
 310
 311                                         rename((fastafileNames[s] + toString(processIDS[0]) + ".temp").c_str(), filteredFasta.c_str());
 312
 313                                         //append fasta files
 314                                         for(int i=1;i<processors;i++){
 315                                                 appendAlignFiles((fastafileNames[s] + toString(processIDS[i]) + ".temp"), filteredFasta);
 316                                                 remove((fastafileNames[s] + toString(processIDS[i]) + ".temp").c_str());
 317                                         }
 318                                 }
 319
 320                                 if (m->control_pressed) {  return 1; }
 321                 #else
 322                                 ifstream inFASTA;
 323                                 openInputFile(fastafileNames[s], inFASTA);
 324                                 int numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 325                                 inFASTA.close();
 326
 327                                 lines.push_back(new linePair(0, numFastaSeqs));
 328
 329                                 driverRunFilter(filter, filteredFasta, fastafileNames[s], lines[0]);
 330
 331                                 if (m->control_pressed) {  return 1; }
 332                 #endif
 333 #endif
 334                         outputNames.push_back(filteredFasta);
 335                 }
 336
 337                 return 0;
 338         }
 339         catch(exception& e) {
 340                 m->errorOut(e, "FilterSeqsCommand", "filterSequences");
 341                 exit(1);
 342         }
 343 }
 344 /**************************************************************************************/
 345 int FilterSeqsCommand::driverMPIRun(istringstream& in, MPI_File& outMPI) {
 346         try {
 347                 string outputString = "";
 348                 int count = 0;
 349                 MPI_Status status;
 350
 351                 while (!in.eof()) {
 352
 353                         Sequence seq(in); gobble(in);
 354
 355                         if (seq.getName() != "") {
 356                                 string align = seq.getAligned();
 357                                 string filterSeq = "";
 358
 359                                 for(int j=0;j<alignmentLength;j++){
 360                                         if(filter[j] == '1'){
 361                                                 filterSeq += align[j];
 362                                         }
 363                                 }
 364
 365                                 count++;
 366                                 outputString += ">" + seq.getName() + "\n" + filterSeq + "\n";
 367
 368                                 if(count % 10 == 0){ //output to file
 369                                         //send results to parent
 370                                         int length = outputString.length();
 371                                         char buf[length];
 372                                         strcpy(buf, outputString.c_str());
 373
 374                                         MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
 375                                         outputString = "";
 376                                 }
 377
 378                         }
 379                 }
 380
 381                 if(outputString != ""){ //output to file
 382                         //send results to parent
 383                         int length = outputString.length();
 384                         char buf[length];
 385                         strcpy(buf, outputString.c_str());
 386
 387                         MPI_File_write_shared(outMPI, buf, length, MPI_CHAR, &status);
 388                         outputString = "";
 389                 }
 390
 391
 392                 return 0;
 393         }
 394         catch(exception& e) {
 395                 m->errorOut(e, "FilterSeqsCommand", "driverRunFilter");
 396                 exit(1);
 397         }
 398 }
 399 /**************************************************************************************/
 400 int FilterSeqsCommand::driverRunFilter(string F, string outputFilename, string inputFilename, linePair* line) {
 401         try {
 402                 ofstream out;
 403                 openOutputFile(outputFilename, out);
 404
 405                 ifstream in;
 406                 openInputFile(inputFilename, in);
 407
 408                 in.seekg(line->start);
 409
 410                 for(int i=0;i<line->numSeqs;i++){
 411
 412                                 if (m->control_pressed) { in.close(); out.close(); return 0; }
 413
 414                                 Sequence seq(in);
 415                                 if (seq.getName() != "") {
 416                                         string align = seq.getAligned();
 417                                         string filterSeq = "";
 418
 419                                         for(int j=0;j<alignmentLength;j++){
 420                                                 if(filter[j] == '1'){
 421                                                         filterSeq += align[j];
 422                                                 }
 423                                         }
 424
 425                                         out << '>' << seq.getName() << endl << filterSeq << endl;
 426                                 }
 427                                 gobble(in);
 428                 }
 429                 out.close();
 430                 in.close();
 431
 432                 return 0;
 433         }
 434         catch(exception& e) {
 435                 m->errorOut(e, "FilterSeqsCommand", "driverRunFilter");
 436                 exit(1);
 437         }
 438 }
 439 /**************************************************************************************************/
 440
 441 int FilterSeqsCommand::createProcessesRunFilter(string F, string filename) {
 442         try {
 443 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 444                 int process = 0;
 445                 int exitCommand = 1;
 446                 processIDS.clear();
 447
 448                 //loop through and create all the processes you want
 449                 while (process != processors) {
 450                         int pid = fork();
 451
 452                         if (pid > 0) {
 453                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 454                                 process++;
 455                         }else if (pid == 0){
 456                                 string filteredFasta = filename + toString(getpid()) + ".temp";
 457                                 driverRunFilter(F, filteredFasta, filename, lines[process]);
 458                                 exit(0);
 459                         }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
 460                 }
 461
 462                 //force parent to wait until all the processes are done
 463                 for (int i=0;i<processors;i++) {
 464                         int temp = processIDS[i];
 465                         wait(&temp);
 466                 }
 467
 468                 return exitCommand;
 469 #endif
 470         }
 471         catch(exception& e) {
 472                 m->errorOut(e, "FilterSeqsCommand", "createProcessesRunFilter");
 473                 exit(1);
 474         }
 475 }
 476 /**************************************************************************************/
 477 string FilterSeqsCommand::createFilter() {
 478         try {
 479                 string filterString = "";
 480
 481                 Filters F;
 482
 483                 if (soft != 0)                  {  F.setSoft(soft);             }
 484                 if (trump != '*')               {  F.setTrump(trump);   }
 485
 486                 F.setLength(alignmentLength);
 487
 488                 if(soft != 0 || isTrue(vertical)){
 489                         F.initialize();
 490                 }
 491
 492                 if(hard.compare("") != 0)       {       F.doHard(hard);         }
 493                 else                                            {       F.setFilter(string(alignmentLength, '1'));      }
 494
 495                 numSeqs = 0;
 496
 497                 if(trump != '*' || isTrue(vertical) || soft != 0){
 498                         for (int s = 0; s < fastafileNames.size(); s++) {
 499
 500 #ifdef USE_MPI
 501                                 int pid, rc, ierr;
 502                                 int Atag = 1; int Ttag = 2; int Ctag = 3; int Gtag = 4; int Gaptag = 5;
 503                                 int tag = 2001;
 504
 505                                 MPI_Status status;
 506                                 MPI_File in;
 507                                 rc = MPI_Comm_size(MPI_COMM_WORLD, &processors);
 508                                 rc = MPI_Comm_rank(MPI_COMM_WORLD, &pid);
 509
 510                                 char tempFileName[fastafileNames[s].length()];
 511                                 strcpy(tempFileName, fastafileNames[s].c_str());
 512                 cout << pid  << " tempFileName " << tempFileName << endl;
 513                                 MPI_File_open(MPI_COMM_WORLD, tempFileName, MPI_MODE_RDONLY, MPI_INFO_NULL, &in);  //comm, filename, mode, info, filepointer
 514                 cout << pid  << " here" << endl;
 515                                 if (pid == 0) { //you are the root process
 516                                                 setLines(fastafileNames[s]);
 517                                 cout << pid  << " after setlines" << endl;
 518                                                 for (int j = 0; j < lines.size(); j++) { //each process
 519                                                         if (j != 0) { //don't send to yourself
 520                                                                 MPI_Send(&lines[j]->start, 1, MPI_INT, j, tag, MPI_COMM_WORLD); //start position in file
 521                                                                 MPI_Send(&bufferSizes[j], 1, MPI_INT, j, tag, MPI_COMM_WORLD); //how bytes for the read
 522                                                         }
 523                                                 }
 524                                 cout << pid << " done sending" << endl;
 525                                                 char buf[bufferSizes[0]];
 526                                                 MPI_File_read_at(in, 0, buf, bufferSizes[0], MPI_CHAR, &status);
 527                         cout << pid << " done reading" << endl;
 528                                                 string tempBuf = buf;
 529                                                 if (tempBuf.length() > bufferSizes[0]) { tempBuf = tempBuf.substr(0, bufferSizes[0]); }
 530
 531                                                 MPICreateFilter(F, tempBuf);
 532                                 cout << pid << "done with mpi create filter " << endl;
 533                                                 if (m->control_pressed) { MPI_File_close(&in); return filterString; }
 534
 535                                                 vector<int> temp; temp.resize(alignmentLength+1);
 536
 537                                                 //get the frequencies from the child processes
 538                                                 for(int i = 0; i < ((processors-1)*5); i++) {
 539                                                         MPI_Recv(&temp[0], (alignmentLength+1), MPI_INT, MPI_ANY_SOURCE, tag, MPI_COMM_WORLD, &status);
 540                                                         int receiveTag = temp[temp.size()-1];  //child process added a int to the end to indicate what letter count this is for
 541
 542                                                         if (receiveTag == Atag) { //you are recieveing the A frequencies
 543                                                                 for (int k = 0; k < alignmentLength; k++) {             F.a[k] += temp[k];      }
 544                                                         }else if (receiveTag == Ttag) { //you are recieveing the T frequencies
 545                                                                 for (int k = 0; k < alignmentLength; k++) {             F.t[k] += temp[k];      }
 546                                                         }else if (receiveTag == Ctag) { //you are recieveing the C frequencies
 547                                                                 for (int k = 0; k < alignmentLength; k++) {             F.c[k] += temp[k];      }
 548                                                         }else if (receiveTag == Gtag) { //you are recieveing the G frequencies
 549                                                                 for (int k = 0; k < alignmentLength; k++) {             F.g[k] += temp[k];      }
 550                                                         }else if (receiveTag == Gaptag) { //you are recieveing the gap frequencies
 551                                                                 for (int k = 0; k < alignmentLength; k++) {             F.gap[k] += temp[k];    }
 552                                                         }
 553                                                 }
 554
 555
 556                                 }else { //i am the child process
 557                         cout << pid << endl;
 558                                         int startPos, bufferSize;
 559                                         ierr = MPI_Recv(&startPos, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 560                                         ierr = MPI_Recv(&bufferSize, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &status);
 561                         cout << pid << '\t' << startPos << '\t' << bufferSize << endl;
 562                                         //send freqs
 563                                         char buf2[bufferSize];
 564                                         MPI_File_read_at( in, startPos, buf2, bufferSize, MPI_CHAR, &status);
 565
 566                                         string tempBuf = buf2;
 567                                         if (tempBuf.length() > bufferSize) { tempBuf = tempBuf.substr(0, bufferSize); }
 568
 569                                         MPICreateFilter(F, tempBuf);
 570                                 cout << pid << "done with mpi create filter " << endl;
 571                                         if (m->control_pressed) { MPI_File_close(&in); return filterString; }
 572
 573                                         //send my fequency counts
 574                                         F.a.push_back(Atag);
 575                                         int ierr = MPI_Send(&(F.a[0]), (alignmentLength+1), MPI_INT, 0, tag, MPI_COMM_WORLD);
 576                                         F.t.push_back(Ttag);
 577                                         ierr = MPI_Send (&(F.t[0]), (alignmentLength+1), MPI_INT, 0, tag, MPI_COMM_WORLD);
 578                                         F.c.push_back(Ctag);
 579                                         ierr = MPI_Send(&(F.c[0]), (alignmentLength+1), MPI_INT, 0, tag, MPI_COMM_WORLD);
 580                                         F.g.push_back(Gtag);
 581                                         ierr = MPI_Send(&(F.g[0]), (alignmentLength+1), MPI_INT, 0, tag, MPI_COMM_WORLD);
 582                                         F.gap.push_back(Gaptag);
 583                                         ierr = MPI_Send(&(F.gap[0]), (alignmentLength+1), MPI_INT, 0, tag, MPI_COMM_WORLD);
 584                                 }
 585
 586                                 MPI_Barrier(MPI_COMM_WORLD);
 587                                 MPI_File_close(&in);
 588
 589 #else
 590                 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 591                                 if(processors == 1){
 592                                         ifstream inFASTA;
 593                                         openInputFile(fastafileNames[s], inFASTA);
 594                                         int numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 595                                         inFASTA.close();
 596
 597                                         numSeqs += numFastaSeqs;
 598
 599                                         for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 600                                         lines.push_back(new linePair(0, numFastaSeqs));
 601
 602                                         driverCreateFilter(F, fastafileNames[s], lines[0]);
 603                                 }else{
 604                                         setLines(fastafileNames[s]);
 605                                         createProcessesCreateFilter(F, fastafileNames[s]);
 606                                 }
 607
 608                                 if (m->control_pressed) {  return filterString; }
 609                 #else
 610                                 ifstream inFASTA;
 611                                 openInputFile(fastafileNames[s], inFASTA);
 612                                 int numFastaSeqs=count(istreambuf_iterator<char>(inFASTA),istreambuf_iterator<char>(), '>');
 613                                 inFASTA.close();
 614
 615                                 numSeqs += numFastaSeqs;
 616
 617                                 for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 618                                 lines.push_back(new linePair(0, numFastaSeqs));
 619
 620                                 driverCreateFilter(F, fastafileNames[s], lines[0]);
 621                                 if (m->control_pressed) {  return filterString; }
 622                 #endif
 623 #endif
 624
 625                         }
 626                 }
 627
 628                 F.setNumSeqs(numSeqs);
 629
 630                 if(isTrue(vertical) == 1)       {       F.doVertical(); }
 631                 if(soft != 0)                           {       F.doSoft();             }
 632
 633                 filterString = F.getFilter();
 634
 635                 return filterString;
 636         }
 637         catch(exception& e) {
 638                 m->errorOut(e, "FilterSeqsCommand", "createFilter");
 639                 exit(1);
 640         }
 641 }
 642 /**************************************************************************************/
 643 int FilterSeqsCommand::driverCreateFilter(Filters& F, string filename, linePair* line) {
 644         try {
 645
 646                 ifstream in;
 647                 openInputFile(filename, in);
 648
 649                 in.seekg(line->start);
 650
 651                 for(int i=0;i<line->numSeqs;i++){
 652
 653                         if (m->control_pressed) { in.close(); return 1; }
 654
 655                         Sequence seq(in);
 656                         if (seq.getName() != "") {
 657                                         if (seq.getAligned().length() != alignmentLength) { m->mothurOut("Sequences are not all the same length, please correct."); m->mothurOutEndLine(); m->control_pressed = true;  }
 658
 659                                         if(trump != '*'){       F.doTrump(seq); }
 660                                         if(isTrue(vertical) || soft != 0){      F.getFreqs(seq);        }
 661                                         cout.flush();
 662                         }
 663
 664                         //report progress
 665                         if((i+1) % 100 == 0){   m->mothurOut(toString(i+1)); m->mothurOutEndLine();             }
 666                 }
 667
 668                 //report progress
 669                 if((line->numSeqs) % 100 != 0){ m->mothurOut(toString(line->numSeqs)); m->mothurOutEndLine();           }
 670
 671                 in.close();
 672
 673                 return 0;
 674         }
 675         catch(exception& e) {
 676                 m->errorOut(e, "FilterSeqsCommand", "driverCreateFilter");
 677                 exit(1);
 678         }
 679 }
 680 /**************************************************************************************/
 681 int FilterSeqsCommand::MPICreateFilter(Filters& F, string input) {
 682         try {
 683
 684                 vector<string> seqStrings;
 685                 parseBuffer(input, seqStrings);
 686
 687                 for(int i=0;i<seqStrings.size();i++){
 688
 689                         if (seqStrings[i].length() != alignmentLength) {  cout << i << '\t' << seqStrings[i].length() << "Sequences are not all the same length, please correct." << endl; m->control_pressed = true;  }
 690
 691                         if (m->control_pressed) { return 1; }
 692
 693                         Sequence seq("", seqStrings[i]);
 694
 695                         if(trump != '*'){       F.doTrump(seq); }
 696                         if(isTrue(vertical) || soft != 0){      F.getFreqs(seq);        }
 697                         cout.flush();
 698
 699                         //report progress
 700                         if((i+1) % 100 == 0){   m->mothurOut(toString(i+1)); m->mothurOutEndLine();             }
 701                 }
 702
 703                 //report progress
 704                 if((seqStrings.size()) % 100 != 0){     m->mothurOut(toString(seqStrings.size())); m->mothurOutEndLine();               }
 705
 706                 return 0;
 707         }
 708         catch(exception& e) {
 709                 m->errorOut(e, "FilterSeqsCommand", "MPICreateFilter");
 710                 exit(1);
 711         }
 712 }
 713
 714 /**************************************************************************************************/
 715
 716 int FilterSeqsCommand::createProcessesCreateFilter(Filters& F, string filename) {
 717         try {
 718 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 719                 int process = 0;
 720                 int exitCommand = 1;
 721                 processIDS.clear();
 722
 723                 //loop through and create all the processes you want
 724                 while (process != processors) {
 725                         int pid = vfork();
 726
 727                         if (pid > 0) {
 728                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 729                                 process++;
 730                         }else if (pid == 0){
 731                                 driverCreateFilter(F, filename, lines[process]);
 732                                 exit(0);
 733                         }else { m->mothurOut("unable to spawn the necessary processes."); m->mothurOutEndLine(); exit(0); }
 734                 }
 735
 736                 //force parent to wait until all the processes are done
 737                 for (int i=0;i<processors;i++) {
 738                         int temp = processIDS[i];
 739                         wait(&temp);
 740                 }
 741
 742                 return exitCommand;
 743 #endif
 744         }
 745         catch(exception& e) {
 746                 m->errorOut(e, "FilterSeqsCommand", "createProcessesCreateFilter");
 747                 exit(1);
 748         }
 749 }
 750 /**************************************************************************************************/
 751
 752 int FilterSeqsCommand::setLines(string filename) {
 753         try {
 754                 for (int i = 0; i < lines.size(); i++) {  delete lines[i];  }  lines.clear();
 755                 vector<long int> positions;
 756                 bufferSizes.clear();
 757
 758                 ifstream inFASTA;
 759                 openInputFile(filename, inFASTA);
 760
 761                 string input;
 762                 while(!inFASTA.eof()){
 763                         input = getline(inFASTA);
 764
 765                         if (input.length() != 0) {
 766                                 if(input[0] == '>'){ long int pos = inFASTA.tellg(); positions.push_back(pos - input.length() - 1);     }
 767                         }
 768                 }
 769                 inFASTA.close();
 770
 771                 int numFastaSeqs = positions.size();
 772
 773                 FILE * pFile;
 774                 long size;
 775
 776                 //get num bytes in file
 777                 pFile = fopen (filename.c_str(),"rb");
 778                 if (pFile==NULL) perror ("Error opening file");
 779                 else{
 780                         fseek (pFile, 0, SEEK_END);
 781                         size=ftell (pFile);
 782                         fclose (pFile);
 783                 }
 784
 785                 numSeqs += numFastaSeqs;
 786
 787                 int numSeqsPerProcessor = numFastaSeqs / processors;
 788
 789                 for (int i = 0; i < processors; i++) {
 790
 791                         long int startPos = positions[ i * numSeqsPerProcessor ];
 792                         if(i == processors - 1){
 793                                 numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;
 794                                 bufferSizes.push_back(size - startPos);
 795                         }else{
 796                                 long int myEnd = positions[ (i+1) * numSeqsPerProcessor ];
 797                                 bufferSizes.push_back(myEnd-startPos);
 798                         }
 799                         lines.push_back(new linePair(startPos, numSeqsPerProcessor));
 800                 }
 801
 802                 return 0;
 803         }
 804         catch(exception& e) {
 805                 m->errorOut(e, "FilterSeqsCommand", "setLines");
 806                 exit(1);
 807         }
 808 }
 809 /**************************************************************************************************/
 810 int FilterSeqsCommand::parseBuffer(string file, vector<string>& seqs) {
 811         try {
 812                 istringstream iss (file); //,istringstream::in
 813                 string name, seqstring;
 814
 815                 while (!iss.eof()) {
 816
 817                         if (m->control_pressed) { return 0; }
 818
 819                         Sequence seq(iss); gobble(iss);
 820
 821                         if (seq.getName() != "") {
 822                                 seqs.push_back(seq.getAligned());
 823                         }
 824                 }
 825
 826                 return 0;
 827         }
 828         catch(exception& e) {
 829                 m->errorOut(e, "FilterSeqsCommand", "parseBuffer");
 830                 exit(1);
 831         }
 832 }
 833 /**************************************************************************************/