prcseqscommand.cpp

   1 //
   2 //  prcseqscommand.cpp
   3 //  Mothur
   4 //
   5 //  Created by Sarah Westcott on 3/14/12.
   6 //  Copyright (c) 2012 Schloss Lab. All rights reserved.
   7 //
   8
   9 #include "pcrseqscommand.h"
  10
  11 //**********************************************************************************************************************
  12 vector<string> PcrSeqsCommand::setParameters(){
  13         try {
  14                 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none","fasta",false,true,true); parameters.push_back(pfasta);
  15                 CommandParameter poligos("oligos", "InputTypes", "", "", "ecolioligos", "none", "none","",false,false,true); parameters.push_back(poligos);
  16         CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","name",false,false,true); parameters.push_back(pname);
  17         CommandParameter pcount("count", "InputTypes", "", "", "NameCount-CountGroup", "none", "none","count",false,false,true); parameters.push_back(pcount);
  18                 CommandParameter pgroup("group", "InputTypes", "", "", "CountGroup", "none", "none","group",false,false,true); parameters.push_back(pgroup);
  19         CommandParameter ptax("taxonomy", "InputTypes", "", "", "none", "none", "none","taxonomy",false,false,true); parameters.push_back(ptax);
  20         CommandParameter pecoli("ecoli", "InputTypes", "", "", "ecolioligos", "none", "none","",false,false); parameters.push_back(pecoli);
  21                 CommandParameter pstart("start", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pstart);
  22                 CommandParameter pend("end", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pend);
  23                 CommandParameter pnomatch("nomatch", "Multiple", "reject-keep", "reject", "", "", "","",false,false); parameters.push_back(pnomatch);
  24         CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(ppdiffs);
  25
  26                 CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false,true); parameters.push_back(pprocessors);
  27                 CommandParameter pkeepprimer("keepprimer", "Boolean", "", "F", "", "", "","",false,false); parameters.push_back(pkeepprimer);
  28         CommandParameter pkeepdots("keepdots", "Boolean", "", "T", "", "", "","",false,false); parameters.push_back(pkeepdots);
  29         CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  30                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  31
  32                 vector<string> myArray;
  33                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  34                 return myArray;
  35         }
  36         catch(exception& e) {
  37                 m->errorOut(e, "PcrSeqsCommand", "setParameters");
  38                 exit(1);
  39         }
  40 }
  41 //**********************************************************************************************************************
  42 string PcrSeqsCommand::getHelpString(){
  43         try {
  44                 string helpString = "";
  45                 helpString += "The pcr.seqs command reads a fasta file.\n";
  46         helpString += "The pcr.seqs command parameters are fasta, oligos, name, group, count, taxonomy, ecoli, start, end, nomatch, pdiffs, processors, keepprimer and keepdots.\n";
  47                 helpString += "The ecoli parameter is used to provide a fasta file containing a single reference sequence (e.g. for e. coli) this must be aligned. Mothur will trim to the start and end positions of the reference sequence.\n";
  48         helpString += "The start parameter allows you to provide a starting position to trim to.\n";
  49         helpString += "The end parameter allows you to provide a ending position to trim from.\n";
  50         helpString += "The nomatch parameter allows you to decide what to do with sequences where the primer is not found. Default=reject, meaning remove from fasta file.  if nomatch=true, then do nothing to sequence.\n";
  51         helpString += "The processors parameter allows you to use multiple processors.\n";
  52         helpString += "The keepprimer parameter allows you to keep the primer, default=false.\n";
  53         helpString += "The keepdots parameter allows you to keep the leading and trailing .'s, default=true.\n";
  54         helpString += "The pdiffs parameter is used to specify the number of differences allowed in the primer. The default is 0.\n";
  55                 helpString += "Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n";
  56                 helpString += "For more details please check out the wiki http://www.mothur.org/wiki/Pcr.seqs .\n";
  57                 return helpString;
  58         }
  59         catch(exception& e) {
  60                 m->errorOut(e, "PcrSeqsCommand", "getHelpString");
  61                 exit(1);
  62         }
  63 }
  64 //**********************************************************************************************************************
  65 string PcrSeqsCommand::getOutputPattern(string type) {
  66     try {
  67         string pattern = "";
  68
  69         if (type == "fasta")            {   pattern = "[filename],pcr,[extension]-[filename],[tag],pcr,[extension]";    }
  70         else if (type == "taxonomy")    {   pattern = "[filename],pcr,[extension]";    }
  71         else if (type == "name")        {   pattern = "[filename],pcr,[extension]";    }
  72         else if (type == "group")       {   pattern = "[filename],pcr,[extension]";    }
  73         else if (type == "count")       {   pattern = "[filename],pcr,[extension]";    }
  74         else if (type == "accnos")      {   pattern = "[filename],bad.accnos";    }
  75         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  76
  77         return pattern;
  78     }
  79     catch(exception& e) {
  80         m->errorOut(e, "PcrSeqsCommand", "getOutputPattern");
  81         exit(1);
  82     }
  83 }
  84 //**********************************************************************************************************************
  85
  86 PcrSeqsCommand::PcrSeqsCommand(){
  87         try {
  88                 abort = true; calledHelp = true;
  89                 setParameters();
  90                 vector<string> tempOutNames;
  91                 outputTypes["fasta"] = tempOutNames;
  92                 outputTypes["taxonomy"] = tempOutNames;
  93                 outputTypes["group"] = tempOutNames;
  94                 outputTypes["name"] = tempOutNames;
  95         outputTypes["count"] = tempOutNames;
  96         outputTypes["accnos"] = tempOutNames;
  97         }
  98         catch(exception& e) {
  99                 m->errorOut(e, "PcrSeqsCommand", "PcrSeqsCommand");
 100                 exit(1);
 101         }
 102 }
 103 //***************************************************************************************************************
 104
 105 PcrSeqsCommand::PcrSeqsCommand(string option)  {
 106         try {
 107
 108                 abort = false; calledHelp = false;
 109
 110                 //allow user to run help
 111                 if(option == "help") { help(); abort = true; calledHelp = true; }
 112                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 113
 114                 else {
 115                         vector<string> myArray = setParameters();
 116
 117                         OptionParser parser(option);
 118                         map<string,string> parameters = parser.getParameters();
 119
 120                         ValidParameters validParameter;
 121                         map<string,string>::iterator it;
 122
 123                         //check to make sure all parameters are valid for command
 124                         for (it = parameters.begin(); it != parameters.end(); it++) {
 125                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 126                         }
 127
 128                         //initialize outputTypes
 129                         vector<string> tempOutNames;
 130                         outputTypes["fasta"] = tempOutNames;
 131                         outputTypes["taxonomy"] = tempOutNames;
 132                         outputTypes["group"] = tempOutNames;
 133                         outputTypes["name"] = tempOutNames;
 134             outputTypes["accnos"] = tempOutNames;
 135             outputTypes["count"] = tempOutNames;
 136
 137                         //if the user changes the input directory command factory will send this info to us in the output parameter
 138                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 139                         if (inputDir == "not found"){   inputDir = "";          }
 140                         else {
 141                                 string path;
 142                                 it = parameters.find("fasta");
 143                                 //user has given a template file
 144                                 if(it != parameters.end()){
 145                                         path = m->hasPath(it->second);
 146                                         //if the user has not given a path then, add inputdir. else leave path alone.
 147                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
 148                                 }
 149
 150                                 it = parameters.find("oligos");
 151                                 //user has given a template file
 152                                 if(it != parameters.end()){
 153                                         path = m->hasPath(it->second);
 154                                         //if the user has not given a path then, add inputdir. else leave path alone.
 155                                         if (path == "") {       parameters["oligos"] = inputDir + it->second;           }
 156                                 }
 157
 158                 it = parameters.find("ecoli");
 159                                 //user has given a template file
 160                                 if(it != parameters.end()){
 161                                         path = m->hasPath(it->second);
 162                                         //if the user has not given a path then, add inputdir. else leave path alone.
 163                                         if (path == "") {       parameters["ecoli"] = inputDir + it->second;            }
 164                                 }
 165
 166                                 it = parameters.find("taxonomy");
 167                                 //user has given a template file
 168                                 if(it != parameters.end()){
 169                                         path = m->hasPath(it->second);
 170                                         //if the user has not given a path then, add inputdir. else leave path alone.
 171                                         if (path == "") {       parameters["taxonomy"] = inputDir + it->second;         }
 172                                 }
 173
 174                                 it = parameters.find("name");
 175                                 //user has given a template file
 176                                 if(it != parameters.end()){
 177                                         path = m->hasPath(it->second);
 178                                         //if the user has not given a path then, add inputdir. else leave path alone.
 179                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 180                                 }
 181
 182                 it = parameters.find("group");
 183                                 //user has given a template file
 184                                 if(it != parameters.end()){
 185                                         path = m->hasPath(it->second);
 186                                         //if the user has not given a path then, add inputdir. else leave path alone.
 187                                         if (path == "") {       parameters["group"] = inputDir + it->second;            }
 188                                 }
 189
 190                 it = parameters.find("count");
 191                                 //user has given a template file
 192                                 if(it != parameters.end()){
 193                                         path = m->hasPath(it->second);
 194                                         //if the user has not given a path then, add inputdir. else leave path alone.
 195                                         if (path == "") {       parameters["count"] = inputDir + it->second;            }
 196                                 }
 197
 198                         }
 199
 200
 201                         //check for required parameters
 202                         fastafile = validParameter.validFile(parameters, "fasta", true);
 203                         if (fastafile == "not found") {
 204                                 fastafile = m->getFastaFile();
 205                                 if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
 206                                 else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
 207                         }else if (fastafile == "not open") { fastafile = ""; abort = true; }
 208                         else { m->setFastaFile(fastafile); }
 209
 210             //if the user changes the output directory command factory will send this info to us in the output parameter
 211                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){  outputDir = m->hasPath(fastafile);      }
 212
 213                         //check for optional parameter and set defaults
 214                         // ...at some point should added some additional type checking...
 215                         string temp;
 216                         temp = validParameter.validFile(parameters, "keepprimer", false);  if (temp == "not found")    {        temp = "f";     }
 217                         keepprimer = m->isTrue(temp);
 218
 219             temp = validParameter.validFile(parameters, "keepdots", false);  if (temp == "not found")    {      temp = "t";     }
 220                         keepdots = m->isTrue(temp);
 221
 222                         temp = validParameter.validFile(parameters, "oligos", true);
 223                         if (temp == "not found"){       oligosfile = "";                }
 224                         else if(temp == "not open"){    oligosfile = ""; abort = true;  }
 225                         else                                    {       oligosfile = temp; m->setOligosFile(oligosfile);                }
 226
 227             ecolifile = validParameter.validFile(parameters, "ecoli", true);
 228                         if (ecolifile == "not found"){  ecolifile = "";         }
 229                         else if(ecolifile == "not open"){       ecolifile = ""; abort = true;   }
 230
 231             namefile = validParameter.validFile(parameters, "name", true);
 232                         if (namefile == "not found"){   namefile = "";          }
 233                         else if(namefile == "not open"){        namefile = ""; abort = true;    }
 234             else { m->setNameFile(namefile); }
 235
 236             groupfile = validParameter.validFile(parameters, "group", true);
 237                         if (groupfile == "not found"){  groupfile = "";         }
 238                         else if(groupfile == "not open"){       groupfile = ""; abort = true;   }
 239             else { m->setGroupFile(groupfile); }
 240
 241             countfile = validParameter.validFile(parameters, "count", true);
 242                         if (countfile == "not open") { countfile = ""; abort = true; }
 243                         else if (countfile == "not found") { countfile = "";  }
 244                         else { m->setCountTableFile(countfile); }
 245
 246             if ((namefile != "") && (countfile != "")) {
 247                 m->mothurOut("[ERROR]: you may only use one of the following: name or count."); m->mothurOutEndLine(); abort = true;
 248             }
 249
 250             if ((groupfile != "") && (countfile != "")) {
 251                 m->mothurOut("[ERROR]: you may only use one of the following: group or count."); m->mothurOutEndLine(); abort=true;
 252             }
 253
 254             taxfile = validParameter.validFile(parameters, "taxonomy", true);
 255                         if (taxfile == "not found"){    taxfile = "";           }
 256                         else if(taxfile == "not open"){ taxfile = ""; abort = true;     }
 257             else { m->setTaxonomyFile(taxfile); }
 258
 259                         temp = validParameter.validFile(parameters, "start", false);    if (temp == "not found") { temp = "-1"; }
 260                         m->mothurConvert(temp, start);
 261
 262             temp = validParameter.validFile(parameters, "end", false);  if (temp == "not found") { temp = "-1"; }
 263                         m->mothurConvert(temp, end);
 264
 265                         temp = validParameter.validFile(parameters, "processors", false);       if (temp == "not found"){       temp = m->getProcessors();      }
 266                         m->setProcessors(temp);
 267                         m->mothurConvert(temp, processors);
 268
 269             temp = validParameter.validFile(parameters, "pdiffs", false);               if (temp == "not found") { temp = "0"; }
 270                         m->mothurConvert(temp, pdiffs);
 271
 272             nomatch = validParameter.validFile(parameters, "nomatch", false);   if (nomatch == "not found") { nomatch = "reject"; }
 273
 274             if ((nomatch != "reject") && (nomatch != "keep")) { m->mothurOut("[ERROR]: " + nomatch + " is not a valid entry for nomatch. Choices are reject and keep.\n");  abort = true; }
 275
 276             //didnt set anything
 277                         if ((oligosfile == "") && (ecolifile == "") && (start == -1) && (end == -1)) {
 278                 m->mothurOut("[ERROR]: You did not set any options. Please provide an oligos or ecoli file, or set start or end.\n"); abort = true;
 279             }
 280
 281             if ((oligosfile == "") && (ecolifile == "") && (start < 0) && (end == -1)) { m->mothurOut("[ERROR]: Invalid start value.\n"); abort = true; }
 282
 283             if ((ecolifile != "") && (start != -1) && (end != -1)) {
 284                 m->mothurOut("[ERROR]: You provided an ecoli file , but set the start or end parameters. Unsure what you intend.  When you provide the ecoli file, mothur thinks you want to use the start and end of the sequence in the ecoli file.\n"); abort = true;
 285             }
 286
 287
 288             if ((oligosfile != "") && (ecolifile != "")) {
 289                  m->mothurOut("[ERROR]: You can not use an ecoli file at the same time as an oligos file.\n"); abort = true;
 290             }
 291
 292                         //check to make sure you didn't forget the name file by mistake
 293                         if (countfile == "") {
 294                 if (namefile == "") {
 295                     vector<string> files; files.push_back(fastafile);
 296                     parser.getNameFile(files);
 297                 }
 298             }
 299                 }
 300
 301         }
 302         catch(exception& e) {
 303                 m->errorOut(e, "PcrSeqsCommand", "PcrSeqsCommand");
 304                 exit(1);
 305         }
 306 }
 307 //***************************************************************************************************************
 308
 309 int PcrSeqsCommand::execute(){
 310         try{
 311
 312                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 313
 314         int start = time(NULL);
 315
 316         string thisOutputDir = outputDir;
 317                 if (outputDir == "") {  thisOutputDir += m->hasPath(fastafile);  }
 318         map<string, string> variables;
 319         variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(fastafile));
 320         variables["[extension]"] = m->getExtension(fastafile);
 321                 string trimSeqFile = getOutputFileName("fasta",variables);
 322                 outputNames.push_back(trimSeqFile); outputTypes["fasta"].push_back(trimSeqFile);
 323         variables["[tag]"] = "scrap";
 324         string badSeqFile = getOutputFileName("fasta",variables);
 325
 326
 327         length = 0;
 328                 if(oligosfile != ""){    readOligos();     if (m->debug) { m->mothurOut("[DEBUG]: read oligos file. numprimers = " + toString(primers.size()) + ", revprimers = " + toString(revPrimer.size()) + ".\n"); } }  if (m->control_pressed) {  return 0; }
 329         if(ecolifile != "") {    readEcoli();      }  if (m->control_pressed) {  return 0; }
 330
 331         vector<unsigned long long> positions;
 332         int numFastaSeqs = 0;
 333 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 334         positions = m->divideFile(fastafile, processors);
 335         for (int i = 0; i < (positions.size()-1); i++) {        lines.push_back(linePair(positions[i], positions[(i+1)]));      }
 336 #else
 337         if (processors == 1) {
 338             lines.push_back(linePair(0, 1000));
 339         }else {
 340             positions = m->setFilePosFasta(fastafile, numFastaSeqs);
 341             if (positions.size() < processors) { processors = positions.size(); }
 342
 343             //figure out how many sequences you have to process
 344             int numSeqsPerProcessor = numFastaSeqs / processors;
 345             for (int i = 0; i < processors; i++) {
 346                 int startIndex =  i * numSeqsPerProcessor;
 347                 if(i == (processors - 1)){      numSeqsPerProcessor = numFastaSeqs - i * numSeqsPerProcessor;   }
 348                 lines.push_back(linePair(positions[startIndex], numSeqsPerProcessor));
 349             }
 350         }
 351 #endif
 352         if (m->control_pressed) {  return 0; }
 353
 354         set<string> badNames;
 355         if(processors == 1) {    numFastaSeqs = driverPcr(fastafile, trimSeqFile, badSeqFile, badNames, lines[0]);   }
 356         else                {    numFastaSeqs = createProcesses(fastafile, trimSeqFile, badSeqFile, badNames);       }
 357
 358                 if (m->control_pressed) {  return 0; }
 359
 360         //don't write or keep if blank
 361         if (badNames.size() != 0)   { writeAccnos(badNames);        }
 362         if (m->isBlank(badSeqFile)) { m->mothurRemove(badSeqFile);  }
 363         else { outputNames.push_back(badSeqFile); outputTypes["fasta"].push_back(badSeqFile); }
 364
 365         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
 366         if (namefile != "")                     {               readName(badNames);             }
 367         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
 368         if (groupfile != "")            {               readGroup(badNames);    }
 369         if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
 370                 if (taxfile != "")                      {               readTax(badNames);              }
 371                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
 372                 if (countfile != "")                    {               readCount(badNames);            }
 373                 if (m->control_pressed) { for (int i = 0; i < outputNames.size(); i++) {        m->mothurRemove(outputNames[i]); } return 0; }
 374
 375         m->mothurOutEndLine();
 376                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 377                 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
 378                 m->mothurOutEndLine();
 379                 m->mothurOutEndLine();
 380
 381                 //set fasta file as new current fastafile
 382                 string current = "";
 383                 itTypes = outputTypes.find("fasta");
 384                 if (itTypes != outputTypes.end()) {
 385                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
 386                 }
 387
 388                 itTypes = outputTypes.find("name");
 389                 if (itTypes != outputTypes.end()) {
 390                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setNameFile(current); }
 391                 }
 392
 393                 itTypes = outputTypes.find("group");
 394                 if (itTypes != outputTypes.end()) {
 395                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setGroupFile(current); }
 396                 }
 397
 398                 itTypes = outputTypes.find("accnos");
 399                 if (itTypes != outputTypes.end()) {
 400                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setAccnosFile(current); }
 401                 }
 402
 403                 itTypes = outputTypes.find("taxonomy");
 404                 if (itTypes != outputTypes.end()) {
 405                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setTaxonomyFile(current); }
 406                 }
 407
 408         itTypes = outputTypes.find("count");
 409                 if (itTypes != outputTypes.end()) {
 410                         if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setCountTableFile(current); }
 411                 }
 412
 413                 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to screen " + toString(numFastaSeqs) + " sequences.");
 414                 m->mothurOutEndLine();
 415
 416
 417                 return 0;
 418
 419         }
 420         catch(exception& e) {
 421                 m->errorOut(e, "PcrSeqsCommand", "execute");
 422                 exit(1);
 423         }
 424 }
 425 /**************************************************************************************************/
 426 int PcrSeqsCommand::createProcesses(string filename, string goodFileName, string badFileName, set<string>& badSeqNames) {
 427         try {
 428
 429         vector<int> processIDS;
 430         int process = 1;
 431                 int num = 0;
 432
 433 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 434
 435                 //loop through and create all the processes you want
 436                 while (process != processors) {
 437                         int pid = fork();
 438
 439                         if (pid > 0) {
 440                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 441                                 process++;
 442                         }else if (pid == 0){
 443                                 num = driverPcr(filename, goodFileName + toString(getpid()) + ".temp", badFileName + toString(getpid()) + ".temp", badSeqNames, lines[process]);
 444
 445                                 //pass numSeqs to parent
 446                                 ofstream out;
 447                                 string tempFile = filename + toString(getpid()) + ".num.temp";
 448                                 m->openOutputFile(tempFile, out);
 449                                 out << num << '\t' << badSeqNames.size() << endl;
 450                 for (set<string>::iterator it = badSeqNames.begin(); it != badSeqNames.end(); it++) {
 451                     out << (*it) << endl;
 452                 }
 453                                 out.close();
 454
 455                                 exit(0);
 456                         }else {
 457                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 458                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 459                                 exit(0);
 460                         }
 461                 }
 462
 463         num = driverPcr(filename, goodFileName, badFileName, badSeqNames, lines[0]);
 464
 465                 //force parent to wait until all the processes are done
 466                 for (int i=0;i<processIDS.size();i++) {
 467                         int temp = processIDS[i];
 468                         wait(&temp);
 469                 }
 470
 471                 for (int i = 0; i < processIDS.size(); i++) {
 472                         ifstream in;
 473                         string tempFile =  filename + toString(processIDS[i]) + ".num.temp";
 474                         m->openInputFile(tempFile, in);
 475             int numBadNames = 0; string name = "";
 476                         if (!in.eof()) { int tempNum = 0; in >> tempNum >> numBadNames; num += tempNum; m->gobble(in); }
 477             for (int j = 0; j < numBadNames; j++) {
 478                 in >> name; m->gobble(in);
 479                 badSeqNames.insert(name);
 480             }
 481                         in.close(); m->mothurRemove(tempFile);
 482
 483             m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
 484             m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
 485
 486             m->appendFiles((badFileName + toString(processIDS[i]) + ".temp"), badFileName);
 487             m->mothurRemove((badFileName + toString(processIDS[i]) + ".temp"));
 488                 }
 489     #else
 490
 491         //////////////////////////////////////////////////////////////////////////////////////////////////////
 492                 //Windows version shared memory, so be careful when passing variables through the sumScreenData struct.
 493                 //Above fork() will clone, so memory is separate, but that's not the case with windows,
 494                 //Taking advantage of shared memory to allow both threads to add info to badSeqNames.
 495                 //////////////////////////////////////////////////////////////////////////////////////////////////////
 496
 497                 vector<pcrData*> pDataArray;
 498                 DWORD   dwThreadIdArray[processors-1];
 499                 HANDLE  hThreadArray[processors-1];
 500
 501                 //Create processor worker threads.
 502                 for( int i=0; i<processors-1; i++ ){
 503
 504             string extension = "";
 505             if (i!=0) {extension += toString(i) + ".temp"; processIDS.push_back(i); }
 506
 507                         // Allocate memory for thread data.
 508                         pcrData* tempPcr = new pcrData(filename, goodFileName+extension, badFileName+extension, m, oligosfile, ecolifile, primers, revPrimer, nomatch, keepprimer, keepdots, start, end, length, pdiffs, lines[i].start, lines[i].end);
 509                         pDataArray.push_back(tempPcr);
 510
 511                         //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
 512                         hThreadArray[i] = CreateThread(NULL, 0, MyPcrThreadFunction, pDataArray[i], 0, &dwThreadIdArray[i]);
 513                 }
 514
 515         //do your part
 516         num = driverPcr(filename, (goodFileName+toString(processors-1)+".temp"), (badFileName+toString(processors-1)+".temp"),badSeqNames, lines[processors-1]);
 517         processIDS.push_back(processors-1);
 518
 519                 //Wait until all threads have terminated.
 520                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 521
 522                 //Close all thread handles and free memory allocations.
 523                 for(int i=0; i < pDataArray.size(); i++){
 524                         num += pDataArray[i]->count;
 525             if (pDataArray[i]->count != pDataArray[i]->fend) {
 526                 m->mothurOut("[ERROR]: process " + toString(i) + " only processed " + toString(pDataArray[i]->count) + " of " + toString(pDataArray[i]->fend) + " sequences assigned to it, quitting. \n"); m->control_pressed = true;
 527             }
 528             for (set<string>::iterator it = pDataArray[i]->badSeqNames.begin(); it != pDataArray[i]->badSeqNames.end(); it++) { badSeqNames.insert(*it);       }
 529                         CloseHandle(hThreadArray[i]);
 530                         delete pDataArray[i];
 531                 }
 532
 533         for (int i = 0; i < processIDS.size(); i++) {
 534             m->appendFiles((goodFileName + toString(processIDS[i]) + ".temp"), goodFileName);
 535             m->mothurRemove((goodFileName + toString(processIDS[i]) + ".temp"));
 536
 537             m->appendFiles((badFileName + toString(processIDS[i]) + ".temp"), badFileName);
 538             m->mothurRemove((badFileName + toString(processIDS[i]) + ".temp"));
 539                 }
 540
 541 #endif
 542
 543         return num;
 544
 545         }
 546         catch(exception& e) {
 547                 m->errorOut(e, "PcrSeqsCommand", "createProcesses");
 548                 exit(1);
 549         }
 550 }
 551
 552 //**********************************************************************************************************************
 553 int PcrSeqsCommand::driverPcr(string filename, string goodFasta, string badFasta, set<string>& badSeqNames, linePair filePos){
 554         try {
 555                 ofstream goodFile;
 556                 m->openOutputFile(goodFasta, goodFile);
 557
 558         ofstream badFile;
 559                 m->openOutputFile(badFasta, badFile);
 560
 561                 ifstream inFASTA;
 562                 m->openInputFile(filename, inFASTA);
 563
 564                 inFASTA.seekg(filePos.start);
 565
 566                 bool done = false;
 567                 int count = 0;
 568         set<int> lengths;
 569
 570         //pdiffs, bdiffs, primers, barcodes, revPrimers
 571         map<string, int> faked;
 572         TrimOligos trim(pdiffs, 0, primers, faked, revPrimer);
 573
 574                 while (!done) {
 575
 576                         if (m->control_pressed) {  break; }
 577
 578                         Sequence currSeq(inFASTA); m->gobble(inFASTA);
 579
 580             string trashCode = "";
 581                         if (currSeq.getName() != "") {
 582
 583                 if (m->debug) { m->mothurOut("[DEBUG]: seq name = " + currSeq.getName() + ".\n"); }
 584
 585                 bool goodSeq = true;
 586                 if (oligosfile != "") {
 587                     map<int, int> mapAligned;
 588                     bool aligned = isAligned(currSeq.getAligned(), mapAligned);
 589
 590                     //process primers
 591                     if (primers.size() != 0) {
 592                         int primerStart = 0; int primerEnd = 0;
 593                         bool good = trim.findForward(currSeq, primerStart, primerEnd);
 594
 595                         if(!good){      if (nomatch == "reject") { goodSeq = false; } trashCode += "f"; }
 596                         else{
 597                             //are you aligned
 598                             if (aligned) {
 599                                 if (!keepprimer)    {
 600                                     if (keepdots)   { currSeq.filterToPos(mapAligned[primerEnd]);   }
 601                                     else            { currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerEnd]));                                              }
 602                                 }
 603                                 else                {
 604                                     if (keepdots)   { currSeq.filterToPos(mapAligned[primerStart]);  }
 605                                     else            { currSeq.setAligned(currSeq.getAligned().substr(mapAligned[primerStart]));                                             }
 606                                 }
 607                                 isAligned(currSeq.getAligned(), mapAligned);
 608                             }else {
 609                                 if (!keepprimer)    { currSeq.setAligned(currSeq.getUnaligned().substr(primerEnd)); }
 610                                 else                { currSeq.setAligned(currSeq.getUnaligned().substr(primerStart)); }
 611                             }
 612                         }
 613                     }
 614
 615                     //process reverse primers
 616                     if (revPrimer.size() != 0) {
 617                         int primerStart = 0; int primerEnd = 0;
 618                         bool good = trim.findReverse(currSeq, primerStart, primerEnd);
 619                         if(!good){      if (nomatch == "reject") { goodSeq = false; } trashCode += "r"; }
 620                         else{
 621                               //are you aligned
 622                             if (aligned) {
 623                                 if (!keepprimer)    {
 624                                     if (keepdots)   { currSeq.filterFromPos(mapAligned[primerStart]); }
 625                                     else            { currSeq.setAligned(currSeq.getAligned().substr(0, mapAligned[primerStart]));   }
 626                                 }
 627                                 else                {
 628                                     if (keepdots)   { currSeq.filterFromPos(mapAligned[primerEnd]); }
 629                                     else            { currSeq.setAligned(currSeq.getAligned().substr(0, mapAligned[primerEnd]));   }
 630                                 }
 631                             }
 632                             else {
 633                                 if (!keepprimer)    { currSeq.setAligned(currSeq.getUnaligned().substr(0, primerStart));   }
 634                                 else                { currSeq.setAligned(currSeq.getUnaligned().substr(0, primerEnd));     }
 635                             }
 636                         }
 637                     }
 638                 }else if (ecolifile != "") {
 639                     //make sure the seqs are aligned
 640                     lengths.insert(currSeq.getAligned().length());
 641                     if (lengths.size() > 1) { m->mothurOut("[ERROR]: seqs are not aligned. When using start and end your sequences must be aligned.\n"); m->control_pressed = true; break; }
 642                     else if (currSeq.getAligned().length() != length) {
 643                         m->mothurOut("[ERROR]: seqs are not the same length as ecoli seq. When using ecoli option your sequences must be aligned and the same length as the ecoli sequence.\n"); m->control_pressed = true; break;
 644                     }else {
 645                         if (keepdots)   {
 646                             currSeq.filterToPos(start);
 647                             currSeq.filterFromPos(end);
 648                         }else {
 649                             string seqString = currSeq.getAligned().substr(0, end);
 650                             seqString = seqString.substr(start);
 651                             currSeq.setAligned(seqString);
 652                         }
 653                     }
 654                 }else{ //using start and end to trim
 655                     //make sure the seqs are aligned
 656                     lengths.insert(currSeq.getAligned().length());
 657                     if (lengths.size() > 1) { m->mothurOut("[ERROR]: seqs are not aligned. When using start and end your sequences must be aligned.\n"); m->control_pressed = true; break; }
 658                     else {
 659                         if (end != -1) {
 660                             if (end > currSeq.getAligned().length()) {  m->mothurOut("[ERROR]: end is longer than your sequence length, aborting.\n"); m->control_pressed = true; break; }
 661                             else {
 662                                 if (keepdots)   { currSeq.filterFromPos(end); }
 663                                 else {
 664                                     string seqString = currSeq.getAligned().substr(0, end);
 665                                     currSeq.setAligned(seqString);
 666                                 }
 667                             }
 668                         }
 669                         if (start != -1) {
 670                             if (keepdots)   {  currSeq.filterToPos(start);  }
 671                             else {
 672                                 string seqString = currSeq.getAligned().substr(start);
 673                                 currSeq.setAligned(seqString);
 674                             }
 675                         }
 676                     }
 677                 }
 678
 679                 //trimming removed all bases
 680                 if (currSeq.getUnaligned() == "") { goodSeq = false; }
 681
 682                                 if(goodSeq == 1)    {   currSeq.printSequence(goodFile);        }
 683                                 else {
 684                     badSeqNames.insert(currSeq.getName());
 685                     currSeq.setName(currSeq.getName() + '|' + trashCode);
 686                     currSeq.printSequence(badFile);
 687                 }
 688                 count++;
 689                         }
 690
 691 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 692             unsigned long long pos = inFASTA.tellg();
 693             if ((pos == -1) || (pos >= filePos.end)) { break; }
 694 #else
 695             if (inFASTA.eof()) { break; }
 696 #endif
 697
 698                         //report progress
 699                         if((count) % 100 == 0){ m->mothurOutJustToScreen("Processing sequence: " + toString(count)+"\n");               }
 700                 }
 701                 //report progress
 702                 if((count) % 100 != 0){ m->mothurOutJustToScreen("Processing sequence: " + toString(count)+"\n");       }
 703
 704         badFile.close();
 705                 goodFile.close();
 706                 inFASTA.close();
 707
 708                 return count;
 709         }
 710         catch(exception& e) {
 711                 m->errorOut(e, "PcrSeqsCommand", "driverPcr");
 712                 exit(1);
 713         }
 714 }
 715 //********************************************************************/
 716 bool PcrSeqsCommand::isAligned(string seq, map<int, int>& aligned){
 717         try {
 718         aligned.clear();
 719         bool isAligned = false;
 720
 721         int countBases = 0;
 722         for (int i = 0; i < seq.length(); i++) {
 723             if (!isalpha(seq[i])) { isAligned = true; }
 724             else { aligned[countBases] = i; countBases++; } //maps location in unaligned -> location in aligned.
 725         }                                                   //ie. the 3rd base may be at spot 10 in the alignment
 726                                                             //later when we trim we want to trim from spot 10.
 727         return isAligned;
 728     }
 729         catch(exception& e) {
 730                 m->errorOut(e, "PcrSeqsCommand", "isAligned");
 731                 exit(1);
 732         }
 733 }
 734 //********************************************************************/
 735 string PcrSeqsCommand::reverseOligo(string oligo){
 736         try {
 737         string reverse = "";
 738
 739         for(int i=oligo.length()-1;i>=0;i--){
 740
 741             if(oligo[i] == 'A')         {       reverse += 'T'; }
 742             else if(oligo[i] == 'T'){   reverse += 'A'; }
 743             else if(oligo[i] == 'U'){   reverse += 'A'; }
 744
 745             else if(oligo[i] == 'G'){   reverse += 'C'; }
 746             else if(oligo[i] == 'C'){   reverse += 'G'; }
 747
 748             else if(oligo[i] == 'R'){   reverse += 'Y'; }
 749             else if(oligo[i] == 'Y'){   reverse += 'R'; }
 750
 751             else if(oligo[i] == 'M'){   reverse += 'K'; }
 752             else if(oligo[i] == 'K'){   reverse += 'M'; }
 753
 754             else if(oligo[i] == 'W'){   reverse += 'W'; }
 755             else if(oligo[i] == 'S'){   reverse += 'S'; }
 756
 757             else if(oligo[i] == 'B'){   reverse += 'V'; }
 758             else if(oligo[i] == 'V'){   reverse += 'B'; }
 759
 760             else if(oligo[i] == 'D'){   reverse += 'H'; }
 761             else if(oligo[i] == 'H'){   reverse += 'D'; }
 762
 763             else                                                {       reverse += 'N'; }
 764         }
 765
 766
 767         return reverse;
 768     }
 769         catch(exception& e) {
 770                 m->errorOut(e, "PcrSeqsCommand", "reverseOligo");
 771                 exit(1);
 772         }
 773 }
 774
 775 //***************************************************************************************************************
 776 bool PcrSeqsCommand::readOligos(){
 777         try {
 778                 ifstream inOligos;
 779                 m->openInputFile(oligosfile, inOligos);
 780
 781                 string type, oligo, group;
 782         int primerCount = 0;
 783
 784                 while(!inOligos.eof()){
 785
 786                         inOligos >> type;
 787
 788                         if(type[0] == '#'){ //ignore
 789                                 while (!inOligos.eof()) {       char c = inOligos.get();  if (c == 10 || c == 13){      break;  }       } // get rest of line if there's any crap there
 790                                 m->gobble(inOligos);
 791                         }else{
 792                                 m->gobble(inOligos);
 793                                 //make type case insensitive
 794                                 for(int i=0;i<type.length();i++){       type[i] = toupper(type[i]);  }
 795
 796                                 inOligos >> oligo;
 797
 798                                 for(int i=0;i<oligo.length();i++){
 799                                         oligo[i] = toupper(oligo[i]);
 800                                         if(oligo[i] == 'U')     {       oligo[i] = 'T'; }
 801                                 }
 802
 803                                 if(type == "FORWARD"){
 804                                         // get rest of line in case there is a primer name
 805                                         while (!inOligos.eof()) {
 806                         char c = inOligos.get();
 807                         if (c == 10 || c == 13 || c == -1){     break;  }
 808                         else if (c == 32 || c == 9){;} //space or tab
 809                                         }
 810                                         primers[oligo] = primerCount; primerCount++;
 811                 }else if(type == "REVERSE"){
 812                     string oligoRC = reverseOligo(oligo);
 813                     revPrimer.push_back(oligoRC);
 814                     //cout << "oligo = " << oligo << " reverse = " << oligoRC << endl;
 815                                 }else if(type == "BARCODE"){
 816                                         inOligos >> group;
 817                                 }else if((type == "LINKER")||(type == "SPACER")) {;}
 818                                 else{   m->mothurOut(type + " is not recognized as a valid type. Choices are forward, reverse, linker, spacer and barcode. Ignoring " + oligo + "."); m->mothurOutEndLine(); m->control_pressed = true; }
 819                         }
 820                         m->gobble(inOligos);
 821                 }
 822                 inOligos.close();
 823
 824                 if ((primers.size() == 0) && (revPrimer.size() == 0)) {
 825                         m->mothurOut("[ERROR]: your oligos file does not contain valid primers or reverse primers.  Please correct."); m->mothurOutEndLine();
 826             m->control_pressed = true;
 827                         return false;
 828                 }
 829
 830         return true;
 831
 832     }catch(exception& e) {
 833                 m->errorOut(e, "PcrSeqsCommand", "readOligos");
 834                 exit(1);
 835         }
 836 }
 837 //***************************************************************************************************************
 838 bool PcrSeqsCommand::readEcoli(){
 839         try {
 840                 ifstream in;
 841                 m->openInputFile(ecolifile, in);
 842
 843         //read seq
 844         if (!in.eof()){
 845             Sequence ecoli(in);
 846             length = ecoli.getAligned().length();
 847             start = ecoli.getStartPos();
 848             end = ecoli.getEndPos();
 849         }else { in.close(); m->control_pressed = true; return false; }
 850         in.close();
 851
 852         return true;
 853     }
 854         catch(exception& e) {
 855                 m->errorOut(e, "PcrSeqsCommand", "readEcoli");
 856                 exit(1);
 857         }
 858
 859 }
 860 //***************************************************************************************************************
 861 int PcrSeqsCommand::writeAccnos(set<string> badNames){
 862         try {
 863                 string thisOutputDir = outputDir;
 864                 if (outputDir == "") {  thisOutputDir += m->hasPath(fastafile);  }
 865         map<string, string> variables;
 866                 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(fastafile));
 867                 string outputFileName = getOutputFileName("accnos",variables);
 868         outputNames.push_back(outputFileName); outputTypes["accnos"].push_back(outputFileName);
 869
 870         ofstream out;
 871         m->openOutputFile(outputFileName, out);
 872
 873         for (set<string>::iterator it = badNames.begin(); it != badNames.end(); it++) {
 874             if (m->control_pressed) { break; }
 875             out << (*it) << endl;
 876         }
 877
 878         out.close();
 879         return 0;
 880     }
 881         catch(exception& e) {
 882                 m->errorOut(e, "PcrSeqsCommand", "writeAccnos");
 883                 exit(1);
 884         }
 885
 886 }
 887 //***************************************************************************************************************
 888 int PcrSeqsCommand::readName(set<string>& names){
 889         try {
 890                 string thisOutputDir = outputDir;
 891                 if (outputDir == "") {  thisOutputDir += m->hasPath(namefile);  }
 892                 map<string, string> variables;
 893                 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(namefile));
 894         variables["[extension]"] = m->getExtension(namefile);
 895                 string outputFileName = getOutputFileName("name", variables);
 896
 897                 ofstream out;
 898                 m->openOutputFile(outputFileName, out);
 899
 900                 ifstream in;
 901                 m->openInputFile(namefile, in);
 902                 string name, firstCol, secondCol;
 903
 904                 bool wroteSomething = false;
 905                 int removedCount = 0;
 906
 907                 while(!in.eof()){
 908                         if (m->control_pressed) { in.close();  out.close();  m->mothurRemove(outputFileName);  return 0; }
 909
 910                         in >> firstCol;         m->gobble(in);
 911                         in >> secondCol;
 912
 913             string savedSecond = secondCol;
 914                         vector<string> parsedNames;
 915                         m->splitAtComma(secondCol, parsedNames);
 916
 917                         vector<string> validSecond;  validSecond.clear();
 918                         for (int i = 0; i < parsedNames.size(); i++) {
 919                                 if (names.count(parsedNames[i]) == 0) {
 920                                         validSecond.push_back(parsedNames[i]);
 921                                 }
 922                         }
 923
 924                         if (validSecond.size() != parsedNames.size()) {  //we want to get rid of someone, so get rid of everyone
 925                                 for (int i = 0; i < parsedNames.size(); i++) {  names.insert(parsedNames[i]);  }
 926                                 removedCount += parsedNames.size();
 927                         }else {
 928                 out << firstCol << '\t' << savedSecond << endl;
 929                 wroteSomething = true;
 930             }
 931                         m->gobble(in);
 932                 }
 933                 in.close();
 934                 out.close();
 935
 936                 if (wroteSomething == false) {  m->mothurOut("Your file contains only sequences from the .accnos file."); m->mothurOutEndLine();  }
 937                 outputTypes["name"].push_back(outputFileName); outputNames.push_back(outputFileName);
 938
 939                 m->mothurOut("Removed " + toString(removedCount) + " sequences from your name file."); m->mothurOutEndLine();
 940
 941                 return 0;
 942         }
 943         catch(exception& e) {
 944                 m->errorOut(e, "PcrSeqsCommand", "readName");
 945                 exit(1);
 946         }
 947 }
 948 //**********************************************************************************************************************
 949 int PcrSeqsCommand::readGroup(set<string> names){
 950         try {
 951                 string thisOutputDir = outputDir;
 952                 if (outputDir == "") {  thisOutputDir += m->hasPath(groupfile);  }
 953                 map<string, string> variables;
 954                 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(groupfile));
 955         variables["[extension]"] = m->getExtension(groupfile);
 956                 string outputFileName = getOutputFileName("group", variables);
 957
 958                 ofstream out;
 959                 m->openOutputFile(outputFileName, out);
 960
 961                 ifstream in;
 962                 m->openInputFile(groupfile, in);
 963                 string name, group;
 964
 965                 bool wroteSomething = false;
 966                 int removedCount = 0;
 967
 968                 while(!in.eof()){
 969                         if (m->control_pressed) { in.close();  out.close();  m->mothurRemove(outputFileName);  return 0; }
 970
 971                         in >> name;                             //read from first column
 972                         in >> group;                    //read from second column
 973
 974                         //if this name is in the accnos file
 975                         if (names.count(name) == 0) {
 976                                 wroteSomething = true;
 977                                 out << name << '\t' << group << endl;
 978                         }else {  removedCount++;  }
 979
 980                         m->gobble(in);
 981                 }
 982                 in.close();
 983                 out.close();
 984
 985                 if (wroteSomething == false) {  m->mothurOut("Your file contains only sequences from the .accnos file."); m->mothurOutEndLine();  }
 986                 outputTypes["group"].push_back(outputFileName); outputNames.push_back(outputFileName);
 987
 988                 m->mothurOut("Removed " + toString(removedCount) + " sequences from your group file."); m->mothurOutEndLine();
 989
 990
 991                 return 0;
 992         }
 993         catch(exception& e) {
 994                 m->errorOut(e, "PcrSeqsCommand", "readGroup");
 995                 exit(1);
 996         }
 997 }
 998 //**********************************************************************************************************************
 999 int PcrSeqsCommand::readTax(set<string> names){
1000         try {
1001                 string thisOutputDir = outputDir;
1002                 if (outputDir == "") {  thisOutputDir += m->hasPath(taxfile);  }
1003                 map<string, string> variables;
1004                 variables["[filename]"] = thisOutputDir + m->getRootName(m->getSimpleName(taxfile));
1005         variables["[extension]"] = m->getExtension(taxfile);
1006                 string outputFileName = getOutputFileName("taxonomy", variables);
1007
1008                 ofstream out;
1009                 m->openOutputFile(outputFileName, out);
1010
1011                 ifstream in;
1012                 m->openInputFile(taxfile, in);
1013                 string name, tax;
1014
1015                 bool wroteSomething = false;
1016                 int removedCount = 0;
1017
1018                 while(!in.eof()){
1019                         if (m->control_pressed) { in.close();  out.close();  m->mothurRemove(outputFileName);  return 0; }
1020
1021                         in >> name;                             //read from first column
1022                         in >> tax;                      //read from second column
1023
1024                         //if this name is in the accnos file
1025                         if (names.count(name) == 0) {
1026                                 wroteSomething = true;
1027                                 out << name << '\t' << tax << endl;
1028                         }else {  removedCount++;  }
1029
1030                         m->gobble(in);
1031                 }
1032                 in.close();
1033                 out.close();
1034
1035                 if (wroteSomething == false) {  m->mothurOut("Your file contains only sequences from the .accnos file."); m->mothurOutEndLine();  }
1036                 outputTypes["taxonomy"].push_back(outputFileName); outputNames.push_back(outputFileName);
1037
1038                 m->mothurOut("Removed " + toString(removedCount) + " sequences from your taxonomy file."); m->mothurOutEndLine();
1039
1040                 return 0;
1041         }
1042         catch(exception& e) {
1043                 m->errorOut(e, "PcrSeqsCommand", "readTax");
1044                 exit(1);
1045         }
1046 }
1047 //***************************************************************************************************************
1048 int PcrSeqsCommand::readCount(set<string> badSeqNames){
1049         try {
1050                 ifstream in;
1051                 m->openInputFile(countfile, in);
1052                 set<string>::iterator it;
1053
1054                 map<string, string> variables;
1055                 variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(countfile));
1056         variables["[extension]"] = m->getExtension(countfile);
1057                 string goodCountFile = getOutputFileName("count", variables);
1058
1059         outputNames.push_back(goodCountFile);  outputTypes["count"].push_back(goodCountFile);
1060                 ofstream goodCountOut;  m->openOutputFile(goodCountFile, goodCountOut);
1061
1062         string headers = m->getline(in); m->gobble(in);
1063         goodCountOut << headers << endl;
1064
1065         string name, rest; int thisTotal, removedCount; removedCount = 0;
1066         bool wroteSomething = false;
1067         while (!in.eof()) {
1068
1069                         if (m->control_pressed) { goodCountOut.close(); in.close(); m->mothurRemove(goodCountFile); return 0; }
1070
1071                         in >> name; m->gobble(in);
1072             in >> thisTotal; m->gobble(in);
1073             rest = m->getline(in); m->gobble(in);
1074
1075                         if (badSeqNames.count(name) != 0) { removedCount+=thisTotal; }
1076                         else{
1077                 wroteSomething = true;
1078                                 goodCountOut << name << '\t' << thisTotal << '\t' << rest << endl;
1079                         }
1080                 }
1081                 in.close();
1082                 goodCountOut.close();
1083
1084         if (m->control_pressed) { m->mothurRemove(goodCountFile);   }
1085
1086         if (wroteSomething == false) {  m->mothurOut("Your count file contains only sequences from the .accnos file."); m->mothurOutEndLine(); }
1087
1088         //check for groups that have been eliminated
1089         CountTable ct;
1090         if (ct.testGroups(goodCountFile)) {
1091             ct.readTable(goodCountFile);
1092             ct.printTable(goodCountFile);
1093         }
1094
1095                 if (m->control_pressed) { m->mothurRemove(goodCountFile);   }
1096
1097         m->mothurOut("Removed " + toString(removedCount) + " sequences from your count file."); m->mothurOutEndLine();
1098
1099
1100                 return 0;
1101
1102         }
1103         catch(exception& e) {
1104                 m->errorOut(e, "PcrSeqsCommand", "readCOunt");
1105                 exit(1);
1106         }
1107 }
1108 /**************************************************************************************/
1109
1110