trimflowscommand.cpp

   1 /*
   2  *  trimflowscommand.cpp
   3  *  Mothur
   4  *
   5  *  Created by Pat Schloss on 12/22/10.
   6  *  Copyright 2010 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "trimflowscommand.h"
  11 #include "needlemanoverlap.hpp"
  12
  13 //**********************************************************************************************************************
  14
  15 vector<string> TrimFlowsCommand::getValidParameters(){
  16         try {
  17                 string Array[] =  {"flow", "maxflows", "minflows",
  18                         "fasta", "minlength", "maxlength", "maxhomop", "signal", "noise"
  19                         "oligos", "pdiffs", "bdiffs", "tdiffs",  "order",
  20                         "allfiles", "processors",
  21                         "outputdir","inputdir"
  22
  23                 };
  24                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  25                 return myArray;
  26         }
  27         catch(exception& e) {
  28                 m->errorOut(e, "TrimFlowsCommand", "getValidParameters");
  29                 exit(1);
  30         }
  31 }
  32
  33 //**********************************************************************************************************************
  34
  35 vector<string> TrimFlowsCommand::getRequiredParameters(){
  36         try {
  37                 string Array[] =  {"flow"};
  38                 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
  39                 return myArray;
  40         }
  41         catch(exception& e) {
  42                 m->errorOut(e, "TrimFlowsCommand", "getRequiredParameters");
  43                 exit(1);
  44         }
  45 }
  46
  47 //**********************************************************************************************************************
  48
  49 vector<string> TrimFlowsCommand::getRequiredFiles(){
  50         try {
  51                 vector<string> myArray;
  52                 return myArray;
  53         }
  54         catch(exception& e) {
  55                 m->errorOut(e, "TrimFlowsCommand", "getRequiredFiles");
  56                 exit(1);
  57         }
  58 }
  59
  60 //**********************************************************************************************************************
  61
  62 TrimFlowsCommand::TrimFlowsCommand(){
  63         try {
  64                 abort = true; calledHelp = true;
  65                 vector<string> tempOutNames;
  66                 outputTypes["flow"] = tempOutNames;
  67                 outputTypes["fasta"] = tempOutNames;
  68         }
  69         catch(exception& e) {
  70                 m->errorOut(e, "TrimFlowsCommand", "TrimFlowsCommand");
  71                 exit(1);
  72         }
  73 }
  74
  75 //***************************************************************************************************************
  76
  77 TrimFlowsCommand::~TrimFlowsCommand(){  /*      do nothing      */      }
  78
  79 //***************************************************************************************************************
  80
  81 void TrimFlowsCommand::help(){
  82         try {
  83                 m->mothurOut("The trim.flows command reads a flowgram file and creates .....\n");
  84                 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n");
  85                 m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/Trim.flows.\n\n");
  86
  87         }
  88         catch(exception& e) {
  89                 m->errorOut(e, "TrimFlowsCommand", "help");
  90                 exit(1);
  91         }
  92 }
  93
  94 //**********************************************************************************************************************
  95
  96 TrimFlowsCommand::TrimFlowsCommand(string option)  {
  97         try {
  98
  99                 abort = false; calledHelp = false;
 100                 comboStarts = 0;
 101
 102                 //allow user to run help
 103                 if(option == "help") { help(); abort = true; calledHelp = true; }
 104
 105                 else {
 106                         //valid paramters for this command
 107                         string AlignArray[] =  {"flow", "maxflows", "minflows",
 108                                 "fasta", "minlength", "maxlength", "maxhomop", "signal", "noise",
 109                                 "oligos", "pdiffs", "bdiffs", "tdiffs", "order",
 110                                 "allfiles", "processors",
 111
 112                                 //                      "group",
 113                                 "outputdir","inputdir"
 114
 115                         };
 116
 117                         vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
 118
 119                         OptionParser parser(option);
 120                         map<string,string> parameters = parser.getParameters();
 121
 122                         ValidParameters validParameter;
 123                         map<string,string>::iterator it;
 124
 125                         //check to make sure all parameters are valid for command
 126                         for (it = parameters.begin(); it != parameters.end(); it++) {
 127                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 128                         }
 129
 130                         //initialize outputTypes
 131                         vector<string> tempOutNames;
 132                         outputTypes["flow"] = tempOutNames;
 133                         outputTypes["fasta"] = tempOutNames;
 134
 135                         //if the user changes the input directory command factory will send this info to us in the output parameter
 136                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 137                         if (inputDir == "not found"){   inputDir = "";          }
 138                         else {
 139                                 string path;
 140                                 it = parameters.find("flow");
 141                                 //user has given a template file
 142                                 if(it != parameters.end()){
 143                                         path = m->hasPath(it->second);
 144                                         //if the user has not given a path then, add inputdir. else leave path alone.
 145                                         if (path == "") {       parameters["flow"] = inputDir + it->second;             }
 146                                 }
 147
 148                                 it = parameters.find("oligos");
 149                                 //user has given a template file
 150                                 if(it != parameters.end()){
 151                                         path = m->hasPath(it->second);
 152                                         //if the user has not given a path then, add inputdir. else leave path alone.
 153                                         if (path == "") {       parameters["oligos"] = inputDir + it->second;           }
 154                                 }
 155
 156
 157 //                              it = parameters.find("group");
 158 //                              //user has given a template file
 159 //                              if(it != parameters.end()){
 160 //                                      path = m->hasPath(it->second);
 161 //                                      //if the user has not given a path then, add inputdir. else leave path alone.
 162 //                                      if (path == "") {       parameters["group"] = inputDir + it->second;            }
 163 //                              }
 164                         }
 165
 166
 167                         //check for required parameters
 168                         flowFileName = validParameter.validFile(parameters, "flow", true);
 169                         if (flowFileName == "not found") { m->mothurOut("flow is a required parameter for the trim.flows command."); m->mothurOutEndLine(); abort = true; }
 170                         else if (flowFileName == "not open") { abort = true; }
 171
 172                         //if the user changes the output directory command factory will send this info to us in the output parameter
 173                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 174                                 outputDir = "";
 175                                 outputDir += m->hasPath(flowFileName); //if user entered a file with a path then preserve it
 176                         }
 177
 178
 179                         //check for optional parameter and set defaults
 180                         // ...at some point should added some additional type checking...
 181
 182                         string temp;
 183                         temp = validParameter.validFile(parameters, "minflows", false); if (temp == "not found") { temp = "360"; }
 184                         convert(temp, minFlows);
 185
 186                         temp = validParameter.validFile(parameters, "maxflows", false); if (temp == "not found") { temp = "720"; }
 187                         convert(temp, maxFlows);
 188
 189
 190                         temp = validParameter.validFile(parameters, "oligos", true);
 191                         if (temp == "not found")        {       oligoFileName = "";             }
 192                         else if(temp == "not open")     {       abort = true;                   }
 193                         else                                            {       oligoFileName = temp;   }
 194
 195                         temp = validParameter.validFile(parameters, "fasta", false);            if (temp == "not found"){       fasta = 0;              }
 196                         else if(m->isTrue(temp))        {       fasta = 1;      }
 197
 198                         temp = validParameter.validFile(parameters, "maxhomop", false);         if (temp == "not found"){       temp = "9";             }
 199                         convert(temp, maxHomoP);
 200
 201                         temp = validParameter.validFile(parameters, "signal", false);           if (temp == "not found"){       temp = "0.50";  }
 202                         convert(temp, signal);
 203
 204                         temp = validParameter.validFile(parameters, "noise", false);            if (temp == "not found"){       temp = "0.70";  }
 205                         convert(temp, noise);
 206
 207                         temp = validParameter.validFile(parameters, "minlength", false);        if (temp == "not found"){       temp = "0";             }
 208                         convert(temp, minLength);
 209
 210                         temp = validParameter.validFile(parameters, "maxlength", false);        if (temp == "not found"){       temp = "0";             }
 211                         convert(temp, maxLength);
 212
 213                         temp = validParameter.validFile(parameters, "bdiffs", false);           if (temp == "not found"){       temp = "0";             }
 214                         convert(temp, bdiffs);
 215
 216                         temp = validParameter.validFile(parameters, "pdiffs", false);           if (temp == "not found"){       temp = "0";             }
 217                         convert(temp, pdiffs);
 218
 219                         temp = validParameter.validFile(parameters, "tdiffs", false);
 220                         if (temp == "not found"){ int tempTotal = pdiffs + bdiffs;  temp = toString(tempTotal); }
 221                         convert(temp, tdiffs);
 222                         if(tdiffs == 0){        tdiffs = bdiffs + pdiffs;       }
 223
 224                         temp = validParameter.validFile(parameters, "allfiles", false);         if (temp == "not found"){ temp = "T";           }
 225                         allFiles = m->isTrue(temp);
 226
 227                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){ temp = "1";           }
 228                         convert(temp, processors);
 229
 230                         flowOrder = validParameter.validFile(parameters, "order", false);
 231                         if (flowOrder == "not found"){ flowOrder = "TACG";              }
 232                         else if(flowOrder.length() != 4){
 233                                 m->mothurOut("The value of the order option must be four bases long\n");
 234                         }
 235
 236                         if(oligoFileName == ""){        allFiles = 0;           }
 237
 238                         numFPrimers = 0;
 239                         numRPrimers = 0;
 240                 }
 241
 242         }
 243         catch(exception& e) {
 244                 m->errorOut(e, "TrimFlowsCommand", "TrimSeqsCommand");
 245                 exit(1);
 246         }
 247 }
 248
 249 //***************************************************************************************************************
 250
 251 int TrimFlowsCommand::execute(){
 252         try{
 253
 254                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 255
 256                 string trimFlowFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "trim.flow";
 257                 outputNames.push_back(trimFlowFileName); outputTypes["flow"].push_back(trimFlowFileName);
 258
 259                 string scrapFlowFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "scrap.flow";
 260                 outputNames.push_back(scrapFlowFileName); outputTypes["flow"].push_back(scrapFlowFileName);
 261
 262                 string fastaFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.fasta";
 263                 if(fasta){
 264                         outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
 265                 }
 266
 267                 vector<unsigned long int> flowFilePos = getFlowFileBreaks();
 268                 for (int i = 0; i < (flowFilePos.size()-1); i++) {
 269                         lines.push_back(new linePair(flowFilePos[i], flowFilePos[(i+1)]));
 270                 }
 271
 272                 vector<vector<string> > barcodePrimerComboFileNames;
 273                 if(oligoFileName != ""){
 274                         getOligos(barcodePrimerComboFileNames);
 275                 }
 276
 277 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 278                 if(processors == 1){
 279                         driverCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames, lines[0]);
 280                 }else{
 281                         createProcessesCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames);
 282                 }
 283 #else
 284                 driverCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames, lines[0]);
 285 #endif
 286
 287                 if (m->control_pressed) {  return 0; }
 288
 289                 string flowFilesFileName;
 290                 ofstream output;
 291
 292                 if(allFiles){
 293
 294                         flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files";
 295                         m->openOutputFile(flowFilesFileName, output);
 296
 297                         for(int i=0;i<barcodePrimerComboFileNames.size();i++){
 298                                 for(int j=0;j<barcodePrimerComboFileNames[0].size();j++){
 299
 300                                         FILE * pFile;
 301                                         unsigned long int size;
 302
 303                                         //get num bytes in file
 304                                         pFile = fopen (barcodePrimerComboFileNames[i][j].c_str(),"rb");
 305                                         if (pFile==NULL) perror ("Error opening file");
 306                                         else{
 307                                                 fseek (pFile, 0, SEEK_END);
 308                                                 size=ftell (pFile);
 309                                                 fclose (pFile);
 310                                         }
 311
 312                                         if(size < 10){
 313                                                 remove(barcodePrimerComboFileNames[i][j].c_str());
 314                                         }
 315                                         else{
 316                                                 output << barcodePrimerComboFileNames[i][j] << endl;
 317                                         }
 318                                 }
 319                         }
 320                         output.close();
 321                 }
 322                 else{
 323                         flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files";
 324                         m->openOutputFile(flowFilesFileName, output);
 325
 326                         output << trimFlowFileName << endl;
 327
 328                         output.close();
 329                 }
 330                 outputTypes["flow.files"].push_back(flowFilesFileName);
 331                 outputNames.push_back(flowFileName);
 332
 333                 //set fasta file as new current fastafile
 334                 string current = "";
 335                 itTypes = outputTypes.find("fasta");
 336                 if (itTypes != outputTypes.end()) {
 337                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
 338                 }
 339
 340                 m->mothurOutEndLine();
 341                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 342                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 343                 m->mothurOutEndLine();
 344
 345                 return 0;
 346         }
 347         catch(exception& e) {
 348                 m->errorOut(e, "TrimSeqsCommand", "execute");
 349                 exit(1);
 350         }
 351 }
 352
 353 //***************************************************************************************************************
 354
 355 int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileName, string scrapFlowFileName, string fastaFileName, vector<vector<string> > barcodePrimerComboFileNames, linePair* line){
 356
 357         try {
 358
 359                 ofstream trimFlowFile;
 360                 m->openOutputFile(trimFlowFileName, trimFlowFile);
 361                 trimFlowFile.setf(ios::fixed, ios::floatfield); trimFlowFile.setf(ios::showpoint);
 362
 363                 ofstream scrapFlowFile;
 364                 m->openOutputFile(scrapFlowFileName, scrapFlowFile);
 365                 scrapFlowFile.setf(ios::fixed, ios::floatfield); scrapFlowFile.setf(ios::showpoint);
 366
 367                 if(line->start == 4){
 368                         scrapFlowFile << maxFlows << endl;
 369                         trimFlowFile << maxFlows << endl;
 370                         if(allFiles){
 371                                 for(int i=0;i<barcodePrimerComboFileNames.size();i++){
 372                                         for(int j=0;j<barcodePrimerComboFileNames[0].size();j++){
 373                                                 //                              barcodePrimerComboFileNames[i][j] += toString(getpid()) + ".temp";
 374                                                 ofstream temp;
 375                                                 m->openOutputFile(barcodePrimerComboFileNames[i][j], temp);
 376                                                 temp << maxFlows << endl;
 377                                                 temp.close();
 378                                         }
 379                                 }
 380                         }
 381                 }
 382
 383                 FlowData flowData(numFlows, signal, noise, maxHomoP, flowOrder);
 384
 385                 ofstream fastaFile;
 386                 if(fasta){      m->openOutputFile(fastaFileName, fastaFile);    }
 387
 388                 ifstream flowFile;
 389                 m->openInputFile(flowFileName, flowFile);
 390
 391                 flowFile.seekg(line->start);
 392
 393                 int count = 0;
 394                 bool moreSeqs = 1;
 395
 396                 while(moreSeqs) {
 397
 398                         int success = 1;
 399                         int currentSeqDiffs = 0;
 400                         string trashCode = "";
 401
 402                         flowData.getNext(flowFile);
 403                         flowData.capFlows(maxFlows);
 404
 405                         Sequence currSeq = flowData.getSequence();
 406                         if(!flowData.hasMinFlows(minFlows)){    //screen to see if sequence is of a minimum number of flows
 407                                 success = 0;
 408                                 trashCode += 'l';
 409                         }
 410
 411                         if(minLength > 0 || maxLength > 0){     //screen to see if sequence is above and below a specific number of bases
 412                                 int seqLength = currSeq.getNumBases();
 413                                 if(seqLength < minLength || seqLength > maxLength){
 414                                         success = 0;
 415                                         trashCode += 'l';
 416                                 }
 417                         }
 418
 419                         int primerIndex = 0;
 420                         int barcodeIndex = 0;
 421
 422                         if(barcodes.size() != 0){
 423                                 success = stripBarcode(currSeq, barcodeIndex);
 424                                 if(success > bdiffs)            {       trashCode += 'b';       }
 425                                 else{ currentSeqDiffs += success;  }
 426                         }
 427
 428                         if(numFPrimers != 0){
 429                                 success = stripForward(currSeq, primerIndex);
 430                                 if(success > pdiffs)            {       trashCode += 'f';       }
 431                                 else{ currentSeqDiffs += success;  }
 432                         }
 433
 434                         if (currentSeqDiffs > tdiffs)   {       trashCode += 't';   }
 435
 436                         if(numRPrimers != 0){
 437                                 success = stripReverse(currSeq);
 438                                 if(!success)                            {       trashCode += 'r';       }
 439                         }
 440
 441                         if(trashCode.length() == 0){
 442
 443                                 flowData.printFlows(trimFlowFile);
 444
 445                                 if(fasta)       {       currSeq.printSequence(fastaFile);       }
 446
 447                                 if(allFiles){
 448                                         ofstream output;
 449                                         m->openOutputFileAppend(barcodePrimerComboFileNames[barcodeIndex][primerIndex], output);
 450                                         output.setf(ios::fixed, ios::floatfield); trimFlowFile.setf(ios::showpoint);
 451
 452                                         flowData.printFlows(output);
 453                                         output.close();
 454                                 }
 455                         }
 456                         else{
 457                                 flowData.printFlows(scrapFlowFile, trashCode);
 458                         }
 459
 460                         count++;
 461
 462                         //report progress
 463                         if((count) % 10000 == 0){       m->mothurOut(toString(count)); m->mothurOutEndLine();           }
 464
 465 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
 466                         unsigned long int pos = flowFile.tellg();
 467
 468                         if ((pos == -1) || (pos >= line->end)) { break; }
 469 #else
 470                         if (flowFile.eof()) { break; }
 471 #endif
 472
 473                 }
 474                 //report progress
 475                 if((count) % 10000 != 0){       m->mothurOut(toString(count)); m->mothurOutEndLine();           }
 476
 477                 trimFlowFile.close();
 478                 scrapFlowFile.close();
 479                 flowFile.close();
 480                 if(fasta){      fastaFile.close();      }
 481
 482                 return count;
 483         }
 484         catch(exception& e) {
 485                 m->errorOut(e, "TrimSeqsCommand", "driverCreateTrim");
 486                 exit(1);
 487         }
 488 }
 489
 490 //***************************************************************************************************************
 491
 492 void TrimFlowsCommand::getOligos(vector<vector<string> >& outFlowFileNames){
 493         try {
 494                 ifstream oligosFile;
 495                 m->openInputFile(oligoFileName, oligosFile);
 496
 497                 string type, oligo, group;
 498
 499                 int indexPrimer = 0;
 500                 int indexBarcode = 0;
 501
 502                 while(!oligosFile.eof()){
 503
 504                         oligosFile >> type; m->gobble(oligosFile);      //get the first column value of the row - is it a comment or a feature we are interested in?
 505
 506                         if(type[0] == '#'){     //igore the line because there's a comment
 507                                 while (!oligosFile.eof())       {       char c = oligosFile.get(); if (c == 10 || c == 13){     break;  }       } // get rest of line if there's any crap there
 508                         }
 509                         else{                           //there's a feature we're interested in
 510
 511                                 for(int i=0;i<type.length();i++){       type[i] = toupper(type[i]);  }                                  //make type case insensitive
 512
 513                                 oligosFile >> oligo;    //get the DNA sequence for the feature
 514
 515                                 for(int i=0;i<oligo.length();i++){      //make type case insensitive and change any U's to T's
 516                                         oligo[i] = toupper(oligo[i]);
 517                                         if(oligo[i] == 'U')     {       oligo[i] = 'T'; }
 518                                 }
 519
 520                                 if(type == "FORWARD"){  //if the feature is a forward primer...
 521                                         group = "";
 522
 523                                         while (!oligosFile.eof())       {       // get rest of line in case there is a primer name = will have the name of the primer
 524                                                 char c = oligosFile.get();
 525                                                 if (c == 10 || c == 13){        break;  }
 526                                                 else if (c == 32 || c == 9){;} //space or tab
 527                                                 else {  group += c;  }
 528                                         }
 529
 530                                         //have we seen this primer already?
 531                                         map<string, int>::iterator itPrimer = primers.find(oligo);
 532                                         if (itPrimer != primers.end()) { m->mothurOut("primer " + oligo + " is in your oligos file already."); m->mothurOutEndLine();  }
 533
 534                                         primers[oligo]=indexPrimer; indexPrimer++;
 535                                         primerNameVector.push_back(group);
 536
 537                                 }
 538                                 else if(type == "REVERSE"){
 539                                         Sequence oligoRC("reverse", oligo);
 540                                         oligoRC.reverseComplement();
 541                                         revPrimer.push_back(oligoRC.getUnaligned());
 542                                 }
 543                                 else if(type == "BARCODE"){
 544                                         oligosFile >> group;
 545
 546                                         //check for repeat barcodes
 547                                         map<string, int>::iterator itBar = barcodes.find(oligo);
 548                                         if (itBar != barcodes.end()) { m->mothurOut("barcode " + oligo + " is in your oligos file already."); m->mothurOutEndLine();  }
 549
 550                                         barcodes[oligo]=indexBarcode; indexBarcode++;
 551                                         barcodeNameVector.push_back(group);
 552                                 }
 553                                 else{
 554                                         m->mothurOut(type + " is not recognized as a valid type. Choices are forward, reverse, and barcode. Ignoring " + oligo + "."); m->mothurOutEndLine();
 555                                 }
 556                         }
 557
 558                         m->gobble(oligosFile);
 559                 }
 560                 oligosFile.close();
 561
 562                 if(barcodeNameVector.size() == 0 && primerNameVector[0] == ""){ allFiles = 0;   }
 563
 564                 //add in potential combos
 565                 if(barcodeNameVector.size() == 0){
 566                         barcodes[""] = 0;
 567                         barcodeNameVector.push_back("");
 568                 }
 569
 570                 if(primerNameVector.size() == 0){
 571                         primers[""] = 0;
 572                         primerNameVector.push_back("");
 573                 }
 574
 575
 576                 outFlowFileNames.resize(barcodeNameVector.size());
 577                 for(int i=0;i<outFlowFileNames.size();i++){
 578                         outFlowFileNames[i].assign(primerNameVector.size(), "");
 579                 }
 580
 581                 if(allFiles){
 582
 583                         for(map<string, int>::iterator itBar = barcodes.begin();itBar != barcodes.end();itBar++){
 584                                 for(map<string, int>::iterator itPrimer = primers.begin();itPrimer != primers.end(); itPrimer++){
 585
 586                                         string primerName = primerNameVector[itPrimer->second];
 587                                         string barcodeName = barcodeNameVector[itBar->second];
 588
 589                                         string comboGroupName = "";
 590                                         string fileName = "";
 591
 592                                         if(primerName == ""){
 593                                                 comboGroupName = barcodeNameVector[itBar->second];
 594                                                 fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + comboGroupName + ".flow";
 595                                         }
 596                                         else{
 597                                                 if(barcodeName == ""){
 598                                                         comboGroupName = primerNameVector[itPrimer->second];
 599                                                 }
 600                                                 else{
 601                                                         comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second];
 602                                                 }
 603                                                 fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + comboGroupName + ".flow";
 604                                         }
 605
 606                                         outputNames.push_back(fileName);
 607                                         outputTypes["flow"].push_back(fileName);
 608                                         outFlowFileNames[itBar->second][itPrimer->second] = fileName;
 609
 610                                         ofstream temp;
 611                                         m->openOutputFile(fileName, temp);
 612                                         temp.close();
 613                                 }
 614                         }
 615                 }
 616
 617                 numFPrimers = primers.size();
 618                 numRPrimers = revPrimer.size();
 619
 620         }
 621         catch(exception& e) {
 622                 m->errorOut(e, "TrimSeqsCommand", "getOligos");
 623                 exit(1);
 624         }
 625 }
 626
 627 //***************************************************************************************************************
 628
 629 int TrimFlowsCommand::stripBarcode(Sequence& seq, int& group){
 630         try {
 631
 632                 string rawSequence = seq.getUnaligned();
 633                 int success = bdiffs + 1;       //guilty until proven innocent
 634
 635                 //can you find the barcode
 636                 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
 637                         string oligo = it->first;
 638                         if(rawSequence.length() < oligo.length()){      //let's just assume that the barcodes are the same length
 639                                 success = bdiffs + 10;                                  //if the sequence is shorter than the barcode then bail out
 640                                 break;
 641                         }
 642
 643                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
 644                                 group = it->second;
 645                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
 646
 647                                 success = 0;
 648                                 break;
 649                         }
 650                 }
 651
 652                 //if you found the barcode or if you don't want to allow for diffs
 653                 if ((bdiffs == 0) || (success == 0)) { return success;  }
 654
 655                 else { //try aligning and see if you can find it
 656
 657                         int maxLength = 0;
 658
 659                         Alignment* alignment;
 660                         if (barcodes.size() > 0) {
 661                                 map<string,int>::iterator it=barcodes.begin();
 662
 663                                 for(it;it!=barcodes.end();it++){
 664                                         if(it->first.length() > maxLength){
 665                                                 maxLength = it->first.length();
 666                                         }
 667                                 }
 668                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1));
 669
 670                         }else{ alignment = NULL; }
 671
 672                         //can you find the barcode
 673                         int minDiff = 1e6;
 674                         int minCount = 1;
 675                         int minGroup = -1;
 676                         int minPos = 0;
 677
 678                         for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
 679                                 string oligo = it->first;
 680                                 //                              int length = oligo.length();
 681
 682                                 if(rawSequence.length() < maxLength){   //let's just assume that the barcodes are the same length
 683                                         success = bdiffs + 10;
 684                                         break;
 685                                 }
 686
 687                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 688                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+bdiffs));
 689                                 oligo = alignment->getSeqAAln();
 690                                 string temp = alignment->getSeqBAln();
 691
 692                                 int alnLength = oligo.length();
 693
 694                                 for(int i=oligo.length()-1;i>=0;i--){
 695                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
 696                                 }
 697                                 oligo = oligo.substr(0,alnLength);
 698                                 temp = temp.substr(0,alnLength);
 699
 700                                 int numDiff = countDiffs(oligo, temp);
 701
 702                                 if(numDiff < minDiff){
 703                                         minDiff = numDiff;
 704                                         minCount = 1;
 705                                         minGroup = it->second;
 706                                         minPos = 0;
 707                                         for(int i=0;i<alnLength;i++){
 708                                                 if(temp[i] != '-'){
 709                                                         minPos++;
 710                                                 }
 711                                         }
 712                                 }
 713                                 else if(numDiff == minDiff){
 714                                         minCount++;
 715                                 }
 716
 717                         }
 718
 719                         if(minDiff > bdiffs)    {       success = minDiff;              }       //no good matches
 720                         else if(minCount > 1)   {       success = bdiffs + 100; }       //can't tell the difference between multiple barcodes
 721                         else{                                                                                                   //use the best match
 722                                 group = minGroup;
 723                                 seq.setUnaligned(rawSequence.substr(minPos));
 724                                 success = minDiff;
 725                         }
 726
 727                         if (alignment != NULL) {  delete alignment;  }
 728
 729                 }
 730
 731                 return success;
 732
 733         }
 734         catch(exception& e) {
 735                 m->errorOut(e, "TrimFlowsCommand", "stripBarcode");
 736                 exit(1);
 737         }
 738
 739 }
 740
 741 //***************************************************************************************************************
 742
 743 int TrimFlowsCommand::stripForward(Sequence& seq, int& group){
 744         try {
 745
 746                 string rawSequence = seq.getUnaligned();
 747                 int success = pdiffs + 1;       //guilty until proven innocent
 748
 749                 //can you find the primer
 750                 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
 751                         string oligo = it->first;
 752                         if(rawSequence.length() < oligo.length()){      //let's just assume that the primers are the same length
 753                                 success = pdiffs + 10;                                  //if the sequence is shorter than the barcode then bail out
 754                                 break;
 755                         }
 756
 757                         if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
 758                                 group = it->second;
 759                                 seq.setUnaligned(rawSequence.substr(oligo.length()));
 760                                 success = 0;
 761                                 break;
 762                         }
 763                 }
 764
 765                 //if you found the barcode or if you don't want to allow for diffs
 766                 if ((pdiffs == 0) || (success == 0)) {  return success;  }
 767
 768                 else { //try aligning and see if you can find it
 769
 770                         int maxLength = 0;
 771
 772                         Alignment* alignment;
 773                         if (primers.size() > 0) {
 774                                 map<string,int>::iterator it=primers.begin();
 775
 776                                 for(it;it!=primers.end();it++){
 777                                         if(it->first.length() > maxLength){
 778                                                 maxLength = it->first.length();
 779                                         }
 780                                 }
 781                                 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+pdiffs+1));
 782
 783                         }else{ alignment = NULL; }
 784
 785                         //can you find the barcode
 786                         int minDiff = 1e6;
 787                         int minCount = 1;
 788                         int minGroup = -1;
 789                         int minPos = 0;
 790
 791                         for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
 792                                 string oligo = it->first;
 793                                 //                              int length = oligo.length();
 794
 795                                 if(rawSequence.length() < maxLength){
 796                                         success = pdiffs + 100;
 797                                         break;
 798                                 }
 799
 800                                 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
 801                                 alignment->align(oligo, rawSequence.substr(0,oligo.length()+pdiffs));
 802                                 oligo = alignment->getSeqAAln();
 803                                 string temp = alignment->getSeqBAln();
 804
 805                                 int alnLength = oligo.length();
 806
 807                                 for(int i=oligo.length()-1;i>=0;i--){
 808                                         if(oligo[i] != '-'){    alnLength = i+1;        break;  }
 809                                 }
 810                                 oligo = oligo.substr(0,alnLength);
 811                                 temp = temp.substr(0,alnLength);
 812
 813                                 int numDiff = countDiffs(oligo, temp);
 814
 815                                 if(numDiff < minDiff){
 816                                         minDiff = numDiff;
 817                                         minCount = 1;
 818                                         minGroup = it->second;
 819                                         minPos = 0;
 820                                         for(int i=0;i<alnLength;i++){
 821                                                 if(temp[i] != '-'){
 822                                                         minPos++;
 823                                                 }
 824                                         }
 825                                 }
 826                                 else if(numDiff == minDiff){
 827                                         minCount++;
 828                                 }
 829
 830                         }
 831
 832                         if(minDiff > pdiffs)    {       success = minDiff;              }       //no good matches
 833                         else if(minCount > 1)   {       success = pdiffs + 10;  }       //can't tell the difference between multiple primers
 834                         else{                                                                                                   //use the best match
 835                                 group = minGroup;
 836                                 seq.setUnaligned(rawSequence.substr(minPos));
 837                                 success = minDiff;
 838                         }
 839
 840                         if (alignment != NULL) {  delete alignment;  }
 841
 842                 }
 843
 844                 return success;
 845
 846         }
 847         catch(exception& e) {
 848                 m->errorOut(e, "TrimFlowsCommand", "stripForward");
 849                 exit(1);
 850         }
 851 }
 852
 853 //***************************************************************************************************************
 854
 855 bool TrimFlowsCommand::stripReverse(Sequence& seq){
 856         try {
 857
 858                 string rawSequence = seq.getUnaligned();
 859                 bool success = 0;       //guilty until proven innocent
 860
 861                 for(int i=0;i<numRPrimers;i++){
 862                         string oligo = revPrimer[i];
 863
 864                         if(rawSequence.length() < oligo.length()){
 865                                 success = 0;
 866                                 break;
 867                         }
 868
 869                         if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
 870                                 seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
 871                                 success = 1;
 872                                 break;
 873                         }
 874                 }
 875
 876                 return success;
 877
 878         }
 879         catch(exception& e) {
 880                 m->errorOut(e, "TrimFlowsCommand", "stripReverse");
 881                 exit(1);
 882         }
 883 }
 884
 885
 886 //***************************************************************************************************************
 887
 888 bool TrimFlowsCommand::compareDNASeq(string oligo, string seq){
 889         try {
 890                 bool success = 1;
 891                 int length = oligo.length();
 892
 893                 for(int i=0;i<length;i++){
 894
 895                         if(oligo[i] != seq[i]){
 896                                 if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C')    {       success = 0;    }
 897                                 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N'))                                {       success = 0;    }
 898                                 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G'))                                    {       success = 0;    }
 899                                 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T'))                                    {       success = 0;    }
 900                                 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A'))                                    {       success = 0;    }
 901                                 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G'))                                    {       success = 0;    }
 902                                 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A'))                                    {       success = 0;    }
 903                                 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G'))                                    {       success = 0;    }
 904                                 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G'))   {       success = 0;    }
 905                                 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G'))   {       success = 0;    }
 906                                 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C'))   {       success = 0;    }
 907                                 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G'))   {       success = 0;    }
 908
 909                                 if(success == 0)        {       break;   }
 910                         }
 911                         else{
 912                                 success = 1;
 913                         }
 914                 }
 915
 916                 return success;
 917         }
 918         catch(exception& e) {
 919                 m->errorOut(e, "TrimFlowsCommand", "compareDNASeq");
 920                 exit(1);
 921         }
 922
 923 }
 924
 925 //***************************************************************************************************************
 926
 927 int TrimFlowsCommand::countDiffs(string oligo, string seq){
 928         try {
 929
 930                 int length = oligo.length();
 931                 int countDiffs = 0;
 932
 933                 for(int i=0;i<length;i++){
 934
 935                         if(oligo[i] != seq[i]){
 936                                 if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C' || oligo[i] == '-' || oligo[i] == '.')      {       countDiffs++;   }
 937                                 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N'))                                {       countDiffs++;   }
 938                                 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G'))                                    {       countDiffs++;   }
 939                                 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T'))                                    {       countDiffs++;   }
 940                                 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A'))                                    {       countDiffs++;   }
 941                                 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G'))                                    {       countDiffs++;   }
 942                                 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A'))                                    {       countDiffs++;   }
 943                                 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G'))                                    {       countDiffs++;   }
 944                                 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G'))   {       countDiffs++;   }
 945                                 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G'))   {       countDiffs++;   }
 946                                 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C'))   {       countDiffs++;   }
 947                                 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G'))   {       countDiffs++;   }
 948                         }
 949
 950                 }
 951
 952                 return countDiffs;
 953         }
 954         catch(exception& e) {
 955                 m->errorOut(e, "TrimFlowsCommand", "countDiffs");
 956                 exit(1);
 957         }
 958
 959 }
 960
 961 /**************************************************************************************************/
 962
 963 vector<unsigned long int> TrimFlowsCommand::getFlowFileBreaks() {
 964
 965         try{
 966
 967                 vector<unsigned long int> filePos;
 968                 filePos.push_back(0);
 969
 970                 FILE * pFile;
 971                 unsigned long int size;
 972
 973                 //get num bytes in file
 974                 pFile = fopen (flowFileName.c_str(),"rb");
 975                 if (pFile==NULL) perror ("Error opening file");
 976                 else{
 977                         fseek (pFile, 0, SEEK_END);
 978                         size=ftell (pFile);
 979                         fclose (pFile);
 980                 }
 981
 982                 //estimate file breaks
 983                 unsigned long int chunkSize = 0;
 984                 chunkSize = size / processors;
 985
 986                 //file too small to divide by processors
 987                 if (chunkSize == 0)  {  processors = 1; filePos.push_back(size); return filePos;        }
 988
 989                 //for each process seekg to closest file break and search for next '>' char. make that the filebreak
 990                 for (int i = 0; i < processors; i++) {
 991                         unsigned long int spot = (i+1) * chunkSize;
 992
 993                         ifstream in;
 994                         m->openInputFile(flowFileName, in);
 995                         in.seekg(spot);
 996
 997                         string dummy = m->getline(in);
 998
 999                         //there was not another sequence before the end of the file
1000                         unsigned long int sanityPos = in.tellg();
1001
1002 //                      if (sanityPos == -1) {  break;  }
1003 //                      else {  filePos.push_back(newSpot);  }
1004                         if (sanityPos == -1) {  break;  }
1005                         else {  filePos.push_back(sanityPos);  }
1006
1007                         in.close();
1008                 }
1009
1010                 //save end pos
1011                 filePos.push_back(size);
1012
1013                 //sanity check filePos
1014                 for (int i = 0; i < (filePos.size()-1); i++) {
1015                         if (filePos[(i+1)] <= filePos[i]) {  filePos.erase(filePos.begin()+(i+1)); i--; }
1016                 }
1017
1018                 ifstream in;
1019                 m->openInputFile(flowFileName, in);
1020                 in >> numFlows;
1021                 m->gobble(in);
1022                 unsigned long int spot = in.tellg();
1023                 filePos[0] = spot;
1024                 in.close();
1025
1026                 processors = (filePos.size() - 1);
1027
1028                 return filePos;
1029         }
1030         catch(exception& e) {
1031                 m->errorOut(e, "TrimSeqsCommand", "getFlowFileBreaks");
1032                 exit(1);
1033         }
1034 }
1035
1036 /**************************************************************************************************/
1037
1038 int TrimFlowsCommand::createProcessesCreateTrim(string flowFileName, string trimFlowFileName, string scrapFlowFileName, string fastaFileName, vector<vector<string> > barcodePrimerComboFileNames){
1039
1040         try {
1041 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1042                 int process = 1;
1043                 int exitCommand = 1;
1044                 processIDS.clear();
1045
1046                 //loop through and create all the processes you want
1047                 while (process != processors) {
1048                         int pid = fork();
1049
1050                         if (pid > 0) {
1051                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
1052                                 process++;
1053                         }else if (pid == 0){
1054
1055                                 vector<vector<string> > tempBarcodePrimerComboFileNames = barcodePrimerComboFileNames;
1056                                 if(allFiles){
1057                                         for(int i=0;i<tempBarcodePrimerComboFileNames.size();i++){
1058                                                 for(int j=0;j<tempBarcodePrimerComboFileNames[0].size();j++){
1059                                                         tempBarcodePrimerComboFileNames[i][j] += toString(getpid()) + ".temp";
1060                                                         ofstream temp;
1061                                                         m->openOutputFile(tempBarcodePrimerComboFileNames[i][j], temp);
1062                                                         temp.close();
1063
1064                                                 }
1065                                         }
1066                                 }
1067                                 driverCreateTrim(flowFileName,
1068                                                                  (trimFlowFileName + toString(getpid()) + ".temp"),
1069                                                                  (scrapFlowFileName + toString(getpid()) + ".temp"),
1070                                                                  (fastaFileName + toString(getpid()) + ".temp"),
1071                                                                  tempBarcodePrimerComboFileNames, lines[process]);
1072
1073                                 exit(0);
1074                         }else {
1075                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1076                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1077                                 exit(0);
1078                         }
1079                 }
1080
1081                 //parent do my part
1082                 ofstream temp;
1083                 m->openOutputFile(trimFlowFileName, temp);
1084                 temp.close();
1085
1086                 m->openOutputFile(scrapFlowFileName, temp);
1087                 temp.close();
1088
1089                 if(fasta){
1090                         m->openOutputFile(fastaFileName, temp);
1091                         temp.close();
1092                 }
1093
1094                 driverCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames, lines[0]);
1095
1096                 //force parent to wait until all the processes are done
1097                 for (int i=0;i<processIDS.size();i++) {
1098                         int temp = processIDS[i];
1099                         wait(&temp);
1100                 }
1101
1102                 //append files
1103                 m->mothurOutEndLine();
1104                 for(int i=0;i<processIDS.size();i++){
1105
1106                         m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
1107
1108                         m->appendFiles((trimFlowFileName + toString(processIDS[i]) + ".temp"), trimFlowFileName);
1109                         remove((trimFlowFileName + toString(processIDS[i]) + ".temp").c_str());
1110 //                      m->mothurOut("\tDone with trim.flow file"); m->mothurOutEndLine();
1111
1112                         m->appendFiles((scrapFlowFileName + toString(processIDS[i]) + ".temp"), scrapFlowFileName);
1113                         remove((scrapFlowFileName + toString(processIDS[i]) + ".temp").c_str());
1114 //                      m->mothurOut("\tDone with scrap.flow file"); m->mothurOutEndLine();
1115
1116                         if(fasta){
1117                                 m->appendFiles((fastaFileName + toString(processIDS[i]) + ".temp"), fastaFileName);
1118                                 remove((fastaFileName + toString(processIDS[i]) + ".temp").c_str());
1119 //                              m->mothurOut("\tDone with flow.fasta file"); m->mothurOutEndLine();
1120                         }
1121                         if(allFiles){
1122                                 for (int j = 0; j < barcodePrimerComboFileNames.size(); j++) {
1123                                         for (int k = 0; k < barcodePrimerComboFileNames[0].size(); k++) {
1124                                                 m->appendFiles((barcodePrimerComboFileNames[j][k] + toString(processIDS[i]) + ".temp"), barcodePrimerComboFileNames[j][k]);
1125                                                 remove((barcodePrimerComboFileNames[j][k] + toString(processIDS[i]) + ".temp").c_str());
1126                                         }
1127                                 }
1128                         }
1129                 }
1130
1131                 return exitCommand;
1132 #endif
1133         }
1134         catch(exception& e) {
1135                 m->errorOut(e, "TrimFlowsCommand", "createProcessesCreateTrim");
1136                 exit(1);
1137         }
1138 }
1139
1140 //***************************************************************************************************************