primerdesigncommand.cpp

   1 //
   2 //  primerdesigncommand.cpp
   3 //  Mothur
   4 //
   5 //  Created by Sarah Westcott on 1/18/13.
   6 //  Copyright (c) 2013 Schloss Lab. All rights reserved.
   7 //
   8
   9 #include "primerdesigncommand.h"
  10
  11 //**********************************************************************************************************************
  12 vector<string> PrimerDesignCommand::setParameters(){
  13         try {
  14                 CommandParameter plabel("label", "String", "", "", "", "", "","",false,false); parameters.push_back(plabel);
  15         CommandParameter plist("list", "InputTypes", "", "", "none", "none", "none","summary-list",false,true,true); parameters.push_back(plist);
  16                 CommandParameter pfasta("fasta", "InputTypes", "", "", "none", "none", "none","",false,true, true); parameters.push_back(pfasta);
  17                 CommandParameter pname("name", "InputTypes", "", "", "NameCount", "none", "none","",false,false,true); parameters.push_back(pname);
  18         CommandParameter pcount("count", "InputTypes", "", "", "NameCount", "none", "none","",false,false,true); parameters.push_back(pcount);
  19         CommandParameter plength("length", "Number", "", "18", "", "", "","",false,false); parameters.push_back(plength);
  20         CommandParameter pmintm("mintm", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pmintm);
  21         CommandParameter pmaxtm("maxtm", "Number", "", "-1", "", "", "","",false,false); parameters.push_back(pmaxtm);
  22         CommandParameter pprocessors("processors", "Number", "", "1", "", "", "","",false,false); parameters.push_back(pprocessors);
  23         CommandParameter potunumber("otulabel", "String", "", "", "", "", "","",false,true,true); parameters.push_back(potunumber);
  24         CommandParameter ppdiffs("pdiffs", "Number", "", "0", "", "", "","",false,false,true); parameters.push_back(ppdiffs);
  25         CommandParameter pcutoff("cutoff", "Number", "", "100", "", "", "","",false,false); parameters.push_back(pcutoff);
  26                 CommandParameter pinputdir("inputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(pinputdir);
  27                 CommandParameter poutputdir("outputdir", "String", "", "", "", "", "","",false,false); parameters.push_back(poutputdir);
  28
  29                 vector<string> myArray;
  30                 for (int i = 0; i < parameters.size(); i++) {   myArray.push_back(parameters[i].name);          }
  31                 return myArray;
  32         }
  33         catch(exception& e) {
  34                 m->errorOut(e, "PrimerDesignCommand", "setParameters");
  35                 exit(1);
  36         }
  37 }
  38 //**********************************************************************************************************************
  39 string PrimerDesignCommand::getHelpString(){
  40         try {
  41                 string helpString = "";
  42                 helpString += "The primer.design allows you to identify sequence fragments that are specific to particular OTUs.\n";
  43                 helpString += "The primer.design command parameters are: list, fasta, name, count, otulabel, cutoff, length, pdiffs, mintm, maxtm, processors and label.\n";
  44                 helpString += "The list parameter allows you to provide a list file and is required.\n";
  45         helpString += "The fasta parameter allows you to provide a fasta file and is required.\n";
  46         helpString += "The name parameter allows you to provide a name file associated with your fasta file.\n";
  47         helpString += "The count parameter allows you to provide a count file associated with your fasta file.\n";
  48         helpString += "The label parameter is used to indicate the label you want to use from your list file.\n";
  49         helpString += "The otulabel parameter is used to indicate the otu you want to use from your list file. It is required.\n";
  50         helpString += "The pdiffs parameter is used to specify the number of differences allowed in the primer. The default is 0.\n";
  51         helpString += "The length parameter is used to indicate the length of the primer. The default is 18.\n";
  52         helpString += "The mintm parameter is used to indicate minimum melting temperature.\n";
  53         helpString += "The maxtm parameter is used to indicate maximum melting temperature.\n";
  54         helpString += "The processors parameter allows you to indicate the number of processors you want to use. Default=1.\n";
  55         helpString += "The cutoff parameter allows you set a percentage of sequences that support the base. For example: cutoff=97 would only return a sequence that only showed ambiguities for bases that were not supported by at least 97% of sequences.\n";
  56                 helpString += "The primer.desing command should be in the following format: primer.design(list=yourListFile, fasta=yourFastaFile, name=yourNameFile)\n";
  57                 helpString += "primer.design(list=final.an.list, fasta=final.fasta, name=final.names, label=0.03)\n";
  58                 return helpString;
  59         }
  60         catch(exception& e) {
  61                 m->errorOut(e, "PrimerDesignCommand", "getHelpString");
  62                 exit(1);
  63         }
  64 }
  65 //**********************************************************************************************************************
  66 string PrimerDesignCommand::getOutputPattern(string type) {
  67     try {
  68         string pattern = "";
  69
  70         if (type == "fasta") {  pattern = "[filename],[distance],otu.cons.fasta"; }
  71         else if (type == "summary") {  pattern = "[filename],[distance],primer.summary"; }
  72         else if (type == "list") {  pattern = "[filename],pick,[extension]"; }
  73         else { m->mothurOut("[ERROR]: No definition for type " + type + " output pattern.\n"); m->control_pressed = true;  }
  74
  75         return pattern;
  76     }
  77     catch(exception& e) {
  78         m->errorOut(e, "PrimerDesignCommand", "getOutputPattern");
  79         exit(1);
  80     }
  81 }
  82 //**********************************************************************************************************************
  83 PrimerDesignCommand::PrimerDesignCommand(){
  84         try {
  85                 abort = true; calledHelp = true;
  86                 setParameters();
  87         vector<string> tempOutNames;
  88                 outputTypes["summary"] = tempOutNames;
  89         outputTypes["fasta"] = tempOutNames;
  90         outputTypes["list"] = tempOutNames;
  91
  92         }
  93         catch(exception& e) {
  94                 m->errorOut(e, "PrimerDesignCommand", "PrimerDesignCommand");
  95                 exit(1);
  96         }
  97 }
  98 //**********************************************************************************************************************
  99 PrimerDesignCommand::PrimerDesignCommand(string option)  {
 100         try {
 101                 abort = false; calledHelp = false;
 102
 103                 //allow user to run help
 104                 if(option == "help") { help(); abort = true; calledHelp = true; }
 105                 else if(option == "citation") { citation(); abort = true; calledHelp = true;}
 106
 107                 else {
 108                         //valid paramters for this command
 109                         vector<string> myArray = setParameters();
 110
 111                         OptionParser parser(option);
 112                         map<string,string> parameters = parser.getParameters();
 113
 114                         ValidParameters validParameter;
 115                         map<string,string>::iterator it;
 116                         //check to make sure all parameters are valid for command
 117                         for (it = parameters.begin(); it != parameters.end(); it++) {
 118                                 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) {  abort = true;  }
 119                         }
 120
 121             vector<string> tempOutNames;
 122             outputTypes["summary"] = tempOutNames;
 123             outputTypes["fasta"] = tempOutNames;
 124             outputTypes["list"] = tempOutNames;
 125
 126                         //if the user changes the input directory command factory will send this info to us in the output parameter
 127                         string inputDir = validParameter.validFile(parameters, "inputdir", false);
 128                         if (inputDir == "not found"){   inputDir = "";          }
 129                         else {
 130                                 string path;
 131                                 it = parameters.find("count");
 132                                 //user has given a template file
 133                                 if(it != parameters.end()){
 134                                         path = m->hasPath(it->second);
 135                                         //if the user has not given a path then, add inputdir. else leave path alone.
 136                                         if (path == "") {       parameters["count"] = inputDir + it->second;            }
 137                                 }
 138
 139                                 it = parameters.find("fasta");
 140                                 //user has given a template file
 141                                 if(it != parameters.end()){
 142                                         path = m->hasPath(it->second);
 143                                         //if the user has not given a path then, add inputdir. else leave path alone.
 144                                         if (path == "") {       parameters["fasta"] = inputDir + it->second;            }
 145                                 }
 146
 147                                 it = parameters.find("name");
 148                                 //user has given a template file
 149                                 if(it != parameters.end()){
 150                                         path = m->hasPath(it->second);
 151                                         //if the user has not given a path then, add inputdir. else leave path alone.
 152                                         if (path == "") {       parameters["name"] = inputDir + it->second;             }
 153                                 }
 154
 155                 it = parameters.find("list");
 156                                 //user has given a template file
 157                                 if(it != parameters.end()){
 158                                         path = m->hasPath(it->second);
 159                                         //if the user has not given a path then, add inputdir. else leave path alone.
 160                                         if (path == "") {       parameters["list"] = inputDir + it->second;             }
 161                                 }
 162             }
 163
 164                         //check for parameters
 165                         namefile = validParameter.validFile(parameters, "name", true);
 166                         if (namefile == "not open") { abort = true; }
 167                         else if (namefile == "not found") { namefile = ""; }
 168                         else { m->setNameFile(namefile); }
 169
 170             countfile = validParameter.validFile(parameters, "count", true);
 171                         if (countfile == "not open") { countfile = ""; abort = true; }
 172                         else if (countfile == "not found") { countfile = "";  }
 173                         else { m->setCountTableFile(countfile); }
 174
 175             //get fastafile - it is required
 176             fastafile = validParameter.validFile(parameters, "fasta", true);
 177                         if (fastafile == "not open") { fastafile = ""; abort=true;  }
 178                         else if (fastafile == "not found") {
 179                 fastafile = m->getFastaFile();
 180                                 if (fastafile != "") { m->mothurOut("Using " + fastafile + " as input file for the fasta parameter."); m->mothurOutEndLine(); }
 181                                 else {  m->mothurOut("You have no current fastafile and the fasta parameter is required."); m->mothurOutEndLine(); abort = true; }
 182             }else  { m->setFastaFile(fastafile); }
 183
 184             //get listfile - it is required
 185             listfile = validParameter.validFile(parameters, "list", true);
 186                         if (listfile == "not open") { listfile = ""; abort=true;  }
 187                         else if (listfile == "not found") {
 188                 listfile = m->getListFile();
 189                                 if (listfile != "") { m->mothurOut("Using " + listfile + " as input file for the list parameter."); m->mothurOutEndLine(); }
 190                                 else {  m->mothurOut("You have no current listfile and the list parameter is required."); m->mothurOutEndLine(); abort = true; }
 191             }else  { m->setListFile(listfile); }
 192
 193
 194                         if ((namefile != "") && (countfile != "")) {
 195                 m->mothurOut("[ERROR]: you may only use one of the following: name or count."); m->mothurOutEndLine(); abort = true;
 196             }
 197
 198
 199             //if the user changes the output directory command factory will send this info to us in the output parameter
 200                         outputDir = validParameter.validFile(parameters, "outputdir", false);           if (outputDir == "not found"){
 201                                 outputDir = m->hasPath(listfile); //if user entered a file with a path then preserve it
 202                         }
 203
 204             string temp = validParameter.validFile(parameters, "cutoff", false);  if (temp == "not found") { temp = "100"; }
 205                         m->mothurConvert(temp, cutoff);
 206
 207             temp = validParameter.validFile(parameters, "pdiffs", false);  if (temp == "not found") { temp = "0"; }
 208                         m->mothurConvert(temp, pdiffs);
 209
 210             temp = validParameter.validFile(parameters, "length", false);  if (temp == "not found") { temp = "18"; }
 211                         m->mothurConvert(temp, length);
 212
 213             temp = validParameter.validFile(parameters, "mintm", false);  if (temp == "not found") { temp = "-1"; }
 214                         m->mothurConvert(temp, minTM);
 215
 216             temp = validParameter.validFile(parameters, "maxtm", false);  if (temp == "not found") { temp = "-1"; }
 217                         m->mothurConvert(temp, maxTM);
 218
 219             otulabel = validParameter.validFile(parameters, "otulabel", false);  if (otulabel == "not found") { temp = ""; }
 220             if (otulabel == "") {  m->mothurOut("[ERROR]: You must provide an OTU label, aborting.\n"); abort = true; }
 221
 222             temp = validParameter.validFile(parameters, "processors", false);   if (temp == "not found"){       temp = m->getProcessors();      }
 223                         m->setProcessors(temp);
 224                         m->mothurConvert(temp, processors);
 225
 226             label = validParameter.validFile(parameters, "label", false);
 227                         if (label == "not found") { label = ""; m->mothurOut("You did not provide a label, I will use the first label in your inputfile."); m->mothurOutEndLine(); label=""; }
 228
 229             if (countfile == "") {
 230                 if (namefile == "") {
 231                     vector<string> files; files.push_back(fastafile);
 232                     parser.getNameFile(files);
 233                 }
 234             }
 235                 }
 236         }
 237         catch(exception& e) {
 238                 m->errorOut(e, "PrimerDesignCommand", "PrimerDesignCommand");
 239                 exit(1);
 240         }
 241 }
 242 //**********************************************************************************************************************
 243 int PrimerDesignCommand::execute(){
 244         try {
 245
 246                 if (abort == true) { if (calledHelp) { return 0; }  return 2;   }
 247
 248         int start = time(NULL);
 249         //////////////////////////////////////////////////////////////////////////////
 250         //              get file inputs                                             //
 251         //////////////////////////////////////////////////////////////////////////////
 252
 253         //reads list file and selects the label the users specified or the first label
 254         getListVector();
 255         vector<string> binLabels = list->getLabels();
 256         int binIndex = findIndex(otulabel, binLabels);
 257         if (binIndex == -1) { m->mothurOut("[ERROR]: You selected an OTU label that is not in your in your list file, quitting.\n"); return 0; }
 258
 259         map<string, int> nameMap;
 260         unsigned long int numSeqs;  //used to sanity check the files. numSeqs = total seqs for namefile and uniques for count.
 261                                     //list file should have all seqs if namefile was used to create it and only uniques in count file was used.
 262
 263         if (namefile != "")         {  nameMap = m->readNames(namefile, numSeqs);       }
 264         else if (countfile != "")   {  nameMap = readCount(numSeqs);                    }
 265         else { numSeqs = list->getNumSeqs();  }
 266
 267         //sanity check
 268         if (numSeqs != list->getNumSeqs()) {
 269             if (namefile != "")         {  m->mothurOut("[ERROR]: Your list file contains " + toString(list->getNumSeqs()) + " sequences, and your name file contains " + toString(numSeqs) + " sequences, aborting. Do you have the correct files? Perhaps you forgot to include the name file when you clustered? \n");   }
 270             else if (countfile != "") {
 271                 m->mothurOut("[ERROR]: Your list file contains " + toString(list->getNumSeqs()) + " sequences, and your count file contains " + toString(numSeqs) + " unique sequences, aborting. Do you have the correct files? Perhaps you forgot to include the count file when you clustered? \n");
 272             }
 273             m->control_pressed = true;
 274         }
 275
 276         if (m->control_pressed) { delete list; return 0; }
 277
 278         //////////////////////////////////////////////////////////////////////////////
 279         //              process data                                                //
 280         //////////////////////////////////////////////////////////////////////////////
 281         m->mothurOut("\nFinding consensus sequences for each otu..."); cout.flush();
 282
 283         vector<Sequence> conSeqs = createProcessesConSeqs(nameMap, numSeqs);
 284
 285         map<string, string> variables;
 286         variables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(listfile));
 287         variables["[distance]"] = list->getLabel();
 288         string consFastaFile = getOutputFileName("fasta", variables);
 289         outputNames.push_back(consFastaFile); outputTypes["fasta"].push_back(consFastaFile);
 290         ofstream out;
 291         m->openOutputFile(consFastaFile, out);
 292         for (int i = 0; i < conSeqs.size(); i++) {  conSeqs[i].printSequence(out);  }
 293         out.close();
 294
 295         m->mothurOut("Done.\n\n");
 296
 297         set<string> primers = getPrimer(conSeqs[binIndex]);
 298
 299         if (m->control_pressed) { delete list; for (int i = 0; i < outputNames.size(); i++) {   m->mothurRemove(outputNames[i]); } return 0; }
 300
 301         string consSummaryFile = getOutputFileName("summary", variables);
 302         outputNames.push_back(consSummaryFile); outputTypes["summary"].push_back(consSummaryFile);
 303         ofstream outSum;
 304         m->openOutputFile(consSummaryFile, outSum);
 305
 306         outSum << "PrimerOtu: " << otulabel << " Members: " << list->get(binIndex) << endl << "Primers\tminTm\tmaxTm" << endl;
 307
 308         //find min and max melting points
 309         vector<double> minTms;
 310         vector<double> maxTms;
 311         string primerString = "";
 312         for (set<string>::iterator it = primers.begin(); it != primers.end();) {
 313
 314             double minTm, maxTm;
 315             findMeltingPoint(*it, minTm, maxTm);
 316             if ((minTM == -1) && (maxTM == -1)) { //user did not set min or max Tm so save this primer
 317                 minTms.push_back(minTm);
 318                 maxTms.push_back(maxTm);
 319                 outSum << *it << '\t' << minTm << '\t' << maxTm << endl;
 320                 it++;
 321             }else if ((minTM == -1) && (maxTm <= maxTM)){ //user set max and no min, keep if below max
 322                 minTms.push_back(minTm);
 323                 maxTms.push_back(maxTm);
 324                 outSum << *it << '\t' << minTm << '\t' << maxTm << endl;
 325                 it++;
 326             }else if ((maxTM == -1) && (minTm >= minTM)){ //user set min and no max, keep if above min
 327                 minTms.push_back(minTm);
 328                 maxTms.push_back(maxTm);
 329                 outSum << *it << '\t' << minTm << '\t' << maxTm << endl;
 330                 it++;
 331             }else if ((maxTm <= maxTM) && (minTm >= minTM)) { //keep if above min and below max
 332                 minTms.push_back(minTm);
 333                 maxTms.push_back(maxTm);
 334                 outSum << *it << '\t' << minTm << '\t' << maxTm << endl;
 335                 it++;
 336             }else { primers.erase(it++);  } //erase because it didn't qualify
 337         }
 338
 339         outSum << "\nOTUNumber\tPrimer\tStart\tEnd\tLength\tMismatches\tminTm\tmaxTm\n";
 340         outSum.close();
 341
 342         //check each otu's conseq for each primer in otunumber
 343         set<int> otuToRemove = createProcesses(consSummaryFile, minTms, maxTms, primers, conSeqs, binIndex);
 344
 345         if (m->control_pressed) { delete list; for (int i = 0; i < outputNames.size(); i++) {   m->mothurRemove(outputNames[i]); } return 0; }
 346
 347         //print new list file
 348         map<string, string> mvariables;
 349         mvariables["[filename]"] = outputDir + m->getRootName(m->getSimpleName(listfile));
 350         mvariables["[extension]"] = m->getExtension(listfile);
 351         string newListFile = getOutputFileName("list", mvariables);
 352         ofstream outListTemp;
 353         m->openOutputFile(newListFile+".temp", outListTemp);
 354
 355         outListTemp << list->getLabel() << '\t' << (list->getNumBins()-otuToRemove.size()) << '\t';
 356         string headers = "label\tnumOtus\t";
 357         for (int j = 0; j < list->getNumBins(); j++) {
 358             if (m->control_pressed) { break; }
 359             //good otus
 360             if (otuToRemove.count(j) == 0) {
 361                 string bin = list->get(j);
 362                 if (bin != "") {  outListTemp << bin << '\t';  headers += binLabels[j] + '\t'; }
 363             }
 364         }
 365         outListTemp << endl;
 366         outListTemp.close();
 367
 368         ofstream outList;
 369         m->openOutputFile(newListFile, outList);
 370         outList << headers << endl;
 371         outList.close();
 372         m->appendFiles(newListFile+".temp", newListFile);
 373         m->mothurRemove(newListFile+".temp");
 374         outputNames.push_back(newListFile); outputTypes["list"].push_back(newListFile);
 375
 376         if (m->control_pressed) { delete list; for (int i = 0; i < outputNames.size(); i++) {   m->mothurRemove(outputNames[i]); } return 0; }
 377
 378         delete list;
 379
 380         m->mothurOut("It took " + toString(time(NULL) - start) + " secs to process " + toString(list->getNumBins()) + " OTUs.\n");
 381
 382
 383         //output files created by command
 384                 m->mothurOutEndLine();
 385                 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
 386                 for (int i = 0; i < outputNames.size(); i++) {  m->mothurOut(outputNames[i]); m->mothurOutEndLine();    }
 387                 m->mothurOutEndLine();
 388
 389         return 0;
 390
 391     }
 392         catch(exception& e) {
 393                 m->errorOut(e, "PrimerDesignCommand", "execute");
 394                 exit(1);
 395         }
 396 }
 397 //********************************************************************/
 398 //used http://www.biophp.org/minitools/melting_temperature/ as a reference to substitute degenerate bases
 399 // in order to find the min and max Tm values.
 400 //Tm =  64.9°C + 41°C x (number of G’s and C’s in the primer – 16.4)/N
 401
 402 /* A = adenine
 403  * C = cytosine
 404  * G = guanine
 405  * T = thymine
 406  * R = G A (purine)
 407  * Y = T C (pyrimidine)
 408  * K = G T (keto)
 409  * M = A C (amino)
 410  * S = G C (strong bonds)
 411  * W = A T (weak bonds)
 412  * B = G T C (all but A)
 413  * D = G A T (all but C)
 414  * H = A C T (all but G)
 415  * V = G C A (all but T)
 416  * N = A G C T (any) */
 417
 418 int PrimerDesignCommand::findMeltingPoint(string primer, double& minTm, double& maxTm){
 419     try {
 420         string minTmprimer = primer;
 421         string maxTmprimer = primer;
 422
 423         //find minimum Tm string substituting for degenerate bases
 424         for (int i = 0; i < minTmprimer.length(); i++) {
 425             minTmprimer[i] = toupper(minTmprimer[i]);
 426
 427             if (minTmprimer[i] == 'Y') { minTmprimer[i] = 'A'; }
 428             else if (minTmprimer[i] == 'R') { minTmprimer[i] = 'A'; }
 429             else if (minTmprimer[i] == 'W') { minTmprimer[i] = 'A'; }
 430             else if (minTmprimer[i] == 'K') { minTmprimer[i] = 'A'; }
 431             else if (minTmprimer[i] == 'M') { minTmprimer[i] = 'A'; }
 432             else if (minTmprimer[i] == 'D') { minTmprimer[i] = 'A'; }
 433             else if (minTmprimer[i] == 'V') { minTmprimer[i] = 'A'; }
 434             else if (minTmprimer[i] == 'H') { minTmprimer[i] = 'A'; }
 435             else if (minTmprimer[i] == 'B') { minTmprimer[i] = 'A'; }
 436             else if (minTmprimer[i] == 'N') { minTmprimer[i] = 'A'; }
 437             else if (minTmprimer[i] == 'S') { minTmprimer[i] = 'G'; }
 438         }
 439
 440         //find maximum Tm string substituting for degenerate bases
 441         for (int i = 0; i < maxTmprimer.length(); i++) {
 442             maxTmprimer[i] = toupper(maxTmprimer[i]);
 443
 444             if (maxTmprimer[i] == 'Y') { maxTmprimer[i] = 'G'; }
 445             else if (maxTmprimer[i] == 'R') { maxTmprimer[i] = 'G'; }
 446             else if (maxTmprimer[i] == 'W') { maxTmprimer[i] = 'A'; }
 447             else if (maxTmprimer[i] == 'K') { maxTmprimer[i] = 'G'; }
 448             else if (maxTmprimer[i] == 'M') { maxTmprimer[i] = 'G'; }
 449             else if (maxTmprimer[i] == 'D') { maxTmprimer[i] = 'G'; }
 450             else if (maxTmprimer[i] == 'V') { maxTmprimer[i] = 'G'; }
 451             else if (maxTmprimer[i] == 'H') { maxTmprimer[i] = 'G'; }
 452             else if (maxTmprimer[i] == 'B') { maxTmprimer[i] = 'G'; }
 453             else if (maxTmprimer[i] == 'N') { maxTmprimer[i] = 'G'; }
 454             else if (maxTmprimer[i] == 'S') { maxTmprimer[i] = 'G'; }
 455         }
 456
 457         int numGC = 0;
 458         for (int i = 0; i < minTmprimer.length(); i++) {
 459             if (minTmprimer[i] == 'G')       { numGC++; }
 460             else if (minTmprimer[i] == 'C')  { numGC++; }
 461         }
 462
 463         minTm = 64.9 + 41 * (numGC - 16.4) / (double) minTmprimer.length();
 464
 465         numGC = 0;
 466         for (int i = 0; i < maxTmprimer.length(); i++) {
 467             if (maxTmprimer[i] == 'G')       { numGC++; }
 468             else if (maxTmprimer[i] == 'C')  { numGC++; }
 469         }
 470
 471         maxTm = 64.9 + 41 * (numGC - 16.4) / (double) maxTmprimer.length();
 472
 473         return 0;
 474     }
 475         catch(exception& e) {
 476                 m->errorOut(e, "PrimerDesignCommand", "findMeltingPoint");
 477                 exit(1);
 478         }
 479 }
 480 //********************************************************************/
 481 //search for a primer over the sequence string
 482 bool PrimerDesignCommand::findPrimer(string rawSequence, string primer, vector<int>& primerStart, vector<int>& primerEnd, vector<int>& mismatches){
 483         try {
 484         bool foundAtLeastOne = false;  //innocent til proven guilty
 485
 486         //look for exact match
 487         if(rawSequence.length() < primer.length()) {  return false;  }
 488
 489         //search for primer
 490         for (int j = 0; j < rawSequence.length()-length; j++){
 491
 492             if (m->control_pressed) {  return foundAtLeastOne; }
 493
 494             string rawChunk = rawSequence.substr(j, length);
 495
 496             int numDiff = countDiffs(primer, rawChunk);
 497
 498             if(numDiff <= pdiffs){
 499                 primerStart.push_back(j);
 500                 primerEnd.push_back(j+length);
 501                 mismatches.push_back(numDiff);
 502                 foundAtLeastOne = true;
 503             }
 504         }
 505
 506                 return foundAtLeastOne;
 507
 508         }
 509         catch(exception& e) {
 510                 m->errorOut(e, "PrimerDesignCommand", "findPrimer");
 511                 exit(1);
 512         }
 513 }
 514 //********************************************************************/
 515 //find all primers for the given sequence
 516 set<string> PrimerDesignCommand::getPrimer(Sequence primerSeq){
 517         try {
 518         set<string> primers;
 519
 520         string rawSequence = primerSeq.getUnaligned();
 521
 522         for (int j = 0; j < rawSequence.length()-length; j++){
 523             if (m->control_pressed) { break; }
 524
 525             string primer = rawSequence.substr(j, length);
 526             primers.insert(primer);
 527         }
 528
 529         return primers;
 530         }
 531         catch(exception& e) {
 532                 m->errorOut(e, "PrimerDesignCommand", "getPrimer");
 533                 exit(1);
 534         }
 535 }
 536 /**************************************************************************************************/
 537 set<int> PrimerDesignCommand::createProcesses(string newSummaryFile, vector<double>& minTms, vector<double>& maxTms, set<string>& primers, vector<Sequence>& conSeqs, int binIndex) {
 538         try {
 539
 540                 vector<int> processIDS;
 541                 int process = 1;
 542         set<int> otusToRemove;
 543         int numBinsProcessed = 0;
 544
 545                 //sanity check
 546         int numBins = conSeqs.size();
 547                 if (numBins < processors) { processors = numBins; }
 548
 549                 //divide the otus between the processors
 550                 vector<linePair> lines;
 551                 int numOtusPerProcessor = numBins / processors;
 552                 for (int i = 0; i < processors; i++) {
 553                         int startIndex =  i * numOtusPerProcessor;
 554                         int endIndex = (i+1) * numOtusPerProcessor;
 555                         if(i == (processors - 1)){      endIndex = numBins;     }
 556                         lines.push_back(linePair(startIndex, endIndex));
 557                 }
 558
 559 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 560
 561                 //loop through and create all the processes you want
 562                 while (process != processors) {
 563                         int pid = fork();
 564
 565                         if (pid > 0) {
 566                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 567                                 process++;
 568                         }else if (pid == 0){
 569                 //clear old file because we append in driver
 570                 m->mothurRemove(newSummaryFile + toString(getpid()) + ".temp");
 571
 572                                 otusToRemove = driver(newSummaryFile + toString(getpid()) + ".temp", minTms, maxTms, primers, conSeqs, lines[process].start, lines[process].end, numBinsProcessed, binIndex);
 573
 574                 string tempFile = toString(getpid()) + ".otus2Remove.temp";
 575                 ofstream outTemp;
 576                 m->openOutputFile(tempFile, outTemp);
 577
 578                 outTemp << numBinsProcessed << endl;
 579                 outTemp << otusToRemove.size() << endl;
 580                 for (set<int>::iterator it = otusToRemove.begin(); it != otusToRemove.end(); it++) { outTemp << *it << endl; }
 581                 outTemp.close();
 582
 583                                 exit(0);
 584                         }else {
 585                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 586                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 587                                 exit(0);
 588                         }
 589                 }
 590
 591                 //do my part
 592                 otusToRemove = driver(newSummaryFile, minTms, maxTms, primers, conSeqs, lines[0].start, lines[0].end, numBinsProcessed, binIndex);
 593
 594                 //force parent to wait until all the processes are done
 595                 for (int i=0;i<processIDS.size();i++) {
 596                         int temp = processIDS[i];
 597                         wait(&temp);
 598                 }
 599
 600         for (int i = 0; i < processIDS.size(); i++) {
 601             string tempFile = toString(processIDS[i]) +  ".otus2Remove.temp";
 602             ifstream intemp;
 603             m->openInputFile(tempFile, intemp);
 604
 605             int num;
 606             intemp >> num; m->gobble(intemp);
 607             if (num != (lines[i+1].end - lines[i+1].start)) { m->mothurOut("[ERROR]: process " + toString(processIDS[i]) + " did not complete processing all OTUs assigned to it, quitting.\n"); m->control_pressed = true; }
 608             intemp >> num; m->gobble(intemp);
 609             for (int k = 0; k < num; k++) {
 610                 int otu;
 611                 intemp >> otu; m->gobble(intemp);
 612                 otusToRemove.insert(otu);
 613             }
 614             intemp.close();
 615             m->mothurRemove(tempFile);
 616         }
 617
 618
 619     #else
 620
 621                 //////////////////////////////////////////////////////////////////////////////////////////////////////
 622                 //Windows version shared memory, so be careful when passing variables through the primerDesignData struct.
 623                 //Above fork() will clone, so memory is separate, but that's not the case with windows,
 624                 //////////////////////////////////////////////////////////////////////////////////////////////////////
 625
 626                 vector<primerDesignData*> pDataArray;
 627                 DWORD   dwThreadIdArray[processors-1];
 628                 HANDLE  hThreadArray[processors-1];
 629
 630                 //Create processor worker threads.
 631                 for( int i=1; i<processors; i++ ){
 632                         // Allocate memory for thread data.
 633                         string extension = toString(i) + ".temp";
 634                         m->mothurRemove(newSummaryFile+extension);
 635
 636                         primerDesignData* tempPrimer = new primerDesignData((newSummaryFile+extension), m, lines[i].start, lines[i].end, minTms, maxTms, primers, conSeqs, pdiffs, binIndex, length, i);
 637                         pDataArray.push_back(tempPrimer);
 638                         processIDS.push_back(i);
 639
 640                         //MySeqSumThreadFunction is in header. It must be global or static to work with the threads.
 641                         //default security attributes, thread function name, argument to thread function, use default creation flags, returns the thread identifier
 642                         hThreadArray[i-1] = CreateThread(NULL, 0, MyPrimerThreadFunction, pDataArray[i-1], 0, &dwThreadIdArray[i-1]);
 643                 }
 644
 645
 646                 //using the main process as a worker saves time and memory
 647                 otusToRemove = driver(newSummaryFile, minTms, maxTms, primers, conSeqs, lines[0].start, lines[0].end, numBinsProcessed, binIndex);
 648
 649                 //Wait until all threads have terminated.
 650                 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
 651
 652                 //Close all thread handles and free memory allocations.
 653                 for(int i=0; i < pDataArray.size(); i++){
 654             for (set<int>::iterator it = pDataArray[i]->otusToRemove.begin(); it != pDataArray[i]->otusToRemove.end(); it++) {
 655                                 otusToRemove.insert(*it);
 656                         }
 657             int num = pDataArray[i]->numBinsProcessed;
 658             if (num != (lines[processIDS[i]].end - lines[processIDS[i]].start)) { m->mothurOut("[ERROR]: process " + toString(processIDS[i]) + " did not complete processing all OTUs assigned to it, quitting.\n"); m->control_pressed = true; }
 659                         CloseHandle(hThreadArray[i]);
 660                         delete pDataArray[i];
 661                 }
 662
 663 #endif
 664
 665                 //append output files
 666                 for(int i=0;i<processIDS.size();i++){
 667                         m->appendFiles((newSummaryFile + toString(processIDS[i]) + ".temp"), newSummaryFile);
 668                         m->mothurRemove((newSummaryFile + toString(processIDS[i]) + ".temp"));
 669                 }
 670
 671                 return otusToRemove;
 672
 673         }
 674         catch(exception& e) {
 675                 m->errorOut(e, "PrimerDesignCommand", "createProcesses");
 676                 exit(1);
 677         }
 678 }
 679 //**********************************************************************************************************************
 680 set<int> PrimerDesignCommand::driver(string summaryFileName, vector<double>& minTms, vector<double>& maxTms, set<string>& primers, vector<Sequence>& conSeqs, int start, int end, int& numBinsProcessed, int binIndex){
 681         try {
 682         set<int> otuToRemove;
 683
 684         ofstream outSum;
 685         m->openOutputFileAppend(summaryFileName, outSum);
 686
 687         for (int i = start; i < end; i++) {
 688
 689             if (m->control_pressed) { break; }
 690
 691             if (i != (binIndex)) {
 692                 int primerIndex = 0;
 693                 for (set<string>::iterator it = primers.begin(); it != primers.end(); it++) {
 694                     vector<int> primerStarts;
 695                     vector<int> primerEnds;
 696                     vector<int> mismatches;
 697
 698                     bool found = findPrimer(conSeqs[i].getUnaligned(), (*it), primerStarts, primerEnds, mismatches);
 699
 700                     //if we found it report to the table
 701                     if (found) {
 702                         for (int j = 0; j < primerStarts.size(); j++) {
 703                             outSum << (i+1) << '\t' << *it << '\t' << primerStarts[j] << '\t' << primerEnds[j] << '\t' << length << '\t' << mismatches[j] << '\t' << minTms[primerIndex] << '\t' << maxTms[primerIndex] << endl;
 704                         }
 705                         otuToRemove.insert(i);
 706                     }
 707                     primerIndex++;
 708                 }
 709             }
 710             numBinsProcessed++;
 711         }
 712         outSum.close();
 713
 714
 715         return otuToRemove;
 716     }
 717         catch(exception& e) {
 718                 m->errorOut(e, "PrimerDesignCommand", "driver");
 719                 exit(1);
 720         }
 721 }
 722 /**************************************************************************************************/
 723 vector< vector< vector<unsigned int> > > PrimerDesignCommand::driverGetCounts(map<string, int>& nameMap, unsigned long int& fastaCount, vector<unsigned int>& otuCounts, unsigned long long& start, unsigned long long& end){
 724     try {
 725         vector< vector< vector<unsigned int> > > counts;
 726         map<string, int> seq2Bin;
 727         alignedLength = 0;
 728
 729         ifstream in;
 730                 m->openInputFile(fastafile, in);
 731
 732                 in.seekg(start);
 733
 734                 bool done = false;
 735                 fastaCount = 0;
 736
 737                 while (!done) {
 738             if (m->control_pressed) { in.close(); return counts; }
 739
 740                         Sequence seq(in); m->gobble(in);
 741
 742                         if (seq.getName() != "") {
 743                 if (fastaCount == 0) { alignedLength = seq.getAligned().length(); initializeCounts(counts, alignedLength, seq2Bin, nameMap, otuCounts); }
 744                 else if (alignedLength != seq.getAligned().length()) {
 745                     m->mothurOut("[ERROR]: your sequences are not all the same length. primer.design requires sequences to be aligned."); m->mothurOutEndLine(); m->control_pressed = true; break;
 746                 }
 747
 748                 int num = 1;
 749                 map<string, int>::iterator itCount;
 750                 if (namefile != "") {
 751                     itCount = nameMap.find(seq.getName());
 752                     if (itCount == nameMap.end()) {  m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file and not in your name file, aborting."); m->mothurOutEndLine(); m->control_pressed = true; break; }
 753                     else { num = itCount->second; }
 754                     fastaCount+=num;
 755                 }else if (countfile != "") {
 756                     itCount = nameMap.find(seq.getName());
 757                     if (itCount == nameMap.end()) {  m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file and not in your count file, aborting."); m->mothurOutEndLine(); m->control_pressed = true; break; }
 758                     else { num = itCount->second; }
 759                     fastaCount++;
 760                 }else {
 761                     fastaCount++;
 762                 }
 763
 764                 //increment counts
 765                 itCount = seq2Bin.find(seq.getName());
 766                 if (itCount == seq2Bin.end()) {
 767                     if ((namefile != "") || (countfile != "")) {
 768                         m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file and not in your list file, aborting. Perhaps you forgot to include your name or count file while clustering.\n"); m->mothurOutEndLine(); m->control_pressed = true; break;
 769                     }else{
 770                         m->mothurOut("[ERROR]: " + seq.getName() + " is in your fasta file and not in your list file, aborting."); m->mothurOutEndLine(); m->control_pressed = true; break;
 771                     }
 772                 }else {
 773                     otuCounts[itCount->second] += num;
 774                     string aligned = seq.getAligned();
 775                     for (int i = 0; i < alignedLength; i++) {
 776                         char base = toupper(aligned[i]);
 777                         if (base == 'A') { counts[itCount->second][i][0]+=num; }
 778                         else if (base == 'T') { counts[itCount->second][i][1]+=num; }
 779                         else if (base == 'G') { counts[itCount->second][i][2]+=num; }
 780                         else if (base == 'C') { counts[itCount->second][i][3]+=num; }
 781                         else { counts[itCount->second][i][4]+=num; }
 782                     }
 783                 }
 784
 785             }
 786
 787 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 788             unsigned long long pos = in.tellg();
 789             if ((pos == -1) || (pos >= end)) { break; }
 790 #else
 791             if (in.eof()) { break; }
 792 #endif
 793                 }
 794
 795                 in.close();
 796
 797         return counts;
 798     }
 799         catch(exception& e) {
 800                 m->errorOut(e, "PrimerDesignCommand", "driverGetCounts");
 801                 exit(1);
 802         }
 803 }
 804 /**************************************************************************************************/
 805 vector<Sequence> PrimerDesignCommand::createProcessesConSeqs(map<string, int>& nameMap, unsigned long int& numSeqs) {
 806         try {
 807         vector< vector< vector<unsigned int> > > counts;
 808         vector<unsigned int> otuCounts;
 809                 vector<int> processIDS;
 810                 int process = 1;
 811         unsigned long int fastaCount = 0;
 812
 813 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
 814
 815         vector<unsigned long long> positions;
 816         vector<fastaLinePair> lines;
 817         positions = m->divideFile(fastafile, processors);
 818         for (int i = 0; i < (positions.size()-1); i++) {        lines.push_back(fastaLinePair(positions[i], positions[(i+1)])); }
 819
 820                 //loop through and create all the processes you want
 821                 while (process != processors) {
 822                         int pid = fork();
 823
 824                         if (pid > 0) {
 825                                 processIDS.push_back(pid);  //create map from line number to pid so you can append files in correct order later
 826                                 process++;
 827                         }else if (pid == 0){
 828                                 counts = driverGetCounts(nameMap, fastaCount, otuCounts, lines[process].start, lines[process].end);
 829
 830                 string tempFile = toString(getpid()) + ".cons_counts.temp";
 831                 ofstream outTemp;
 832                 m->openOutputFile(tempFile, outTemp);
 833
 834                 outTemp << fastaCount << endl;
 835                 //pass counts
 836                 outTemp << counts.size() << endl;
 837                 for (int i = 0; i < counts.size(); i++) {
 838                     outTemp << counts[i].size() << endl;
 839                     for (int j = 0; j < counts[i].size(); j++) {
 840                         for (int k = 0; k < 5; k++) {  outTemp << counts[i][j][k] << '\t'; }
 841                         outTemp << endl;
 842                     }
 843                 }
 844                 //pass otuCounts
 845                 outTemp << otuCounts.size() << endl;
 846                 for (int i = 0; i < otuCounts.size(); i++) { outTemp << otuCounts[i] << '\t'; }
 847                 outTemp << endl;
 848                 outTemp.close();
 849
 850                                 exit(0);
 851                         }else {
 852                                 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
 853                                 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
 854                                 exit(0);
 855                         }
 856                 }
 857
 858                 //do my part
 859                 counts = driverGetCounts(nameMap, fastaCount, otuCounts, lines[0].start, lines[0].end);
 860
 861                 //force parent to wait until all the processes are done
 862                 for (int i=0;i<processIDS.size();i++) {
 863                         int temp = processIDS[i];
 864                         wait(&temp);
 865                 }
 866
 867         for (int i = 0; i < processIDS.size(); i++) {
 868             string tempFile = toString(processIDS[i]) +  ".cons_counts.temp";
 869             ifstream intemp;
 870             m->openInputFile(tempFile, intemp);
 871
 872             unsigned long int num;
 873             intemp >> num; m->gobble(intemp); fastaCount += num;
 874             intemp >> num; m->gobble(intemp);
 875             if (num != counts.size()) { m->mothurOut("[ERROR]: " + tempFile + " was not built correctly by the child process, quitting.\n"); m->control_pressed = true; }
 876             else {
 877                 //read counts
 878                 for (int k = 0; k < num; k++) {
 879                     int alength;
 880                     intemp >> alength; m->gobble(intemp);
 881                     if (alength != alignedLength) {  m->mothurOut("[ERROR]: your sequences are not all the same length. primer.design requires sequences to be aligned."); m->mothurOutEndLine(); m->control_pressed = true; }
 882                     else {
 883                         for (int j = 0; j < alength; j++) {
 884                             for (int l = 0; l < 5; l++) {  unsigned int numTemp; intemp >> numTemp; m->gobble(intemp); counts[k][j][l] += numTemp;  }
 885                         }
 886                     }
 887                 }
 888                 //read otuCounts
 889                 intemp >> num; m->gobble(intemp);
 890                 for (int k = 0; k < num; k++) {
 891                     unsigned int numTemp; intemp >> numTemp; m->gobble(intemp);
 892                     otuCounts[k] += numTemp;
 893                 }
 894             }
 895             intemp.close();
 896             m->mothurRemove(tempFile);
 897         }
 898
 899
 900 #else
 901         unsigned long long start = 0;
 902         unsigned long long end = 1000;
 903                 counts = driverGetCounts(nameMap, fastaCount, otuCounts, start, end);
 904 #endif
 905
 906         //you will have a nameMap error if there is a namefile or countfile, but if those aren't given we want to make sure the fasta and list file match.
 907         if (fastaCount != numSeqs) {
 908             if ((namefile == "") && (countfile == ""))   {  m->mothurOut("[ERROR]: Your list file contains " + toString(list->getNumSeqs()) + " sequences, and your fasta file contains " + toString(fastaCount) + " sequences, aborting. Do you have the correct files? Perhaps you forgot to include the name or count file? \n");   }
 909             m->control_pressed = true;
 910         }
 911
 912                 vector<Sequence> conSeqs;
 913
 914         if (m->control_pressed) { return conSeqs; }
 915
 916                 //build consensus seqs
 917         string snumBins = toString(counts.size());
 918         for (int i = 0; i < counts.size(); i++) {
 919             if (m->control_pressed) { break; }
 920
 921             string otuLabel = "Otu";
 922             string sbinNumber = toString(i+1);
 923             if (sbinNumber.length() < snumBins.length()) {
 924                 int diff = snumBins.length() - sbinNumber.length();
 925                 for (int h = 0; h < diff; h++) { otuLabel += "0"; }
 926             }
 927             otuLabel += sbinNumber;
 928
 929             string cons = "";
 930             for (int j = 0; j < counts[i].size(); j++) {
 931                 cons += getBase(counts[i][j], otuCounts[i]);
 932             }
 933             Sequence consSeq(otuLabel, cons);
 934             conSeqs.push_back(consSeq);
 935         }
 936
 937         if (m->control_pressed) { conSeqs.clear(); return conSeqs; }
 938
 939         return conSeqs;
 940
 941
 942         }
 943         catch(exception& e) {
 944                 m->errorOut(e, "PrimerDesignCommand", "createProcessesConSeqs");
 945                 exit(1);
 946         }
 947 }
 948 //***************************************************************************************************************
 949
 950 char PrimerDesignCommand::getBase(vector<unsigned int> counts, int size){  //A,T,G,C,Gap
 951         try{
 952                 /* A = adenine
 953          * C = cytosine
 954          * G = guanine
 955          * T = thymine
 956          * R = G A (purine)
 957          * Y = T C (pyrimidine)
 958          * K = G T (keto)
 959          * M = A C (amino)
 960          * S = G C (strong bonds)
 961          * W = A T (weak bonds)
 962          * B = G T C (all but A)
 963          * D = G A T (all but C)
 964          * H = A C T (all but G)
 965          * V = G C A (all but T)
 966          * N = A G C T (any) */
 967
 968                 char conBase = 'N';
 969
 970                 //zero out counts that don't make the cutoff
 971                 float percentage = (100.0 - cutoff) / 100.0;
 972
 973                 for (int i = 0; i < counts.size(); i++) {
 974             float countPercentage = counts[i] / (float) size;
 975                         if (countPercentage < percentage) { counts[i] = 0; }
 976                 }
 977
 978                 //any
 979                 if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] != 0)) {  conBase = 'n'; }
 980                 //any no gap
 981                 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] == 0)) {  conBase = 'N'; }
 982                 //all but T
 983                 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] != 0)) {  conBase = 'v'; }
 984                 //all but T no gap
 985                 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] == 0)) {  conBase = 'V'; }
 986                 //all but G
 987                 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] != 0)) {  conBase = 'h'; }
 988                 //all but G no gap
 989                 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] == 0)) {  conBase = 'H'; }
 990                 //all but C
 991                 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] != 0)) {  conBase = 'd'; }
 992                 //all but C no gap
 993                 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] == 0)) {  conBase = 'D'; }
 994                 //all but A
 995                 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] != 0)) {  conBase = 'b'; }
 996                 //all but A no gap
 997                 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] == 0)) {  conBase = 'B'; }
 998                 //W = A T (weak bonds)
 999                 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] != 0)) {  conBase = 'w'; }
1000                 //W = A T (weak bonds) no gap
1001                 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] == 0)) {  conBase = 'W'; }
1002                 //S = G C (strong bonds)
1003                 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] != 0)) {  conBase = 's'; }
1004                 //S = G C (strong bonds) no gap
1005                 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] == 0)) {  conBase = 'S'; }
1006                 //M = A C (amino)
1007                 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] != 0)) {  conBase = 'm'; }
1008                 //M = A C (amino) no gap
1009                 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] == 0)) {  conBase = 'M'; }
1010                 //K = G T (keto)
1011                 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] != 0)) {  conBase = 'k'; }
1012                 //K = G T (keto) no gap
1013                 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] == 0)) {  conBase = 'K'; }
1014                 //Y = T C (pyrimidine)
1015                 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] != 0)) {  conBase = 'y'; }
1016                 //Y = T C (pyrimidine) no gap
1017                 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] == 0)) {  conBase = 'Y'; }
1018                 //R = G A (purine)
1019                 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] != 0)) {  conBase = 'r'; }
1020                 //R = G A (purine) no gap
1021                 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] == 0)) {  conBase = 'R'; }
1022                 //only A
1023                 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] != 0)) {  conBase = 'a'; }
1024                 //only A no gap
1025                 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] == 0)) {  conBase = 'A'; }
1026                 //only T
1027                 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] != 0)) {  conBase = 't'; }
1028                 //only T no gap
1029                 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] == 0)) {  conBase = 'T'; }
1030                 //only G
1031                 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] != 0)) {  conBase = 'g'; }
1032                 //only G no gap
1033                 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] == 0)) {  conBase = 'G'; }
1034                 //only C
1035                 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] != 0)) {  conBase = 'c'; }
1036                 //only C no gap
1037                 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] == 0)) {  conBase = 'C'; }
1038                 //only gap
1039                 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] != 0)) {  conBase = '-'; }
1040                 //cutoff removed all counts
1041                 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] == 0)) {  conBase = 'N'; }
1042                 else{ m->mothurOut("[ERROR]: cannot find consensus base."); m->mothurOutEndLine(); }
1043
1044                 return conBase;
1045
1046         }
1047         catch(exception& e) {
1048                 m->errorOut(e, "PrimerDesignCommand", "getBase");
1049                 exit(1);
1050         }
1051 }
1052
1053 //**********************************************************************************************************************
1054 int PrimerDesignCommand::initializeCounts(vector< vector< vector<unsigned int> > >& counts, int length, map<string, int>& seq2Bin, map<string, int>& nameMap, vector<unsigned int>& otuCounts){
1055         try {
1056         counts.clear();
1057         otuCounts.clear();
1058         seq2Bin.clear();
1059
1060         //vector< vector< vector<unsigned int> > > counts - otu < spot_in_alignment < counts_for_A,T,G,C,Gap > > >
1061         for (int i = 0; i < list->getNumBins(); i++) {
1062             string binNames = list->get(i);
1063             vector<string> names;
1064             m->splitAtComma(binNames, names);
1065             otuCounts.push_back(0);
1066
1067             //lets be smart and only map the unique names if a name or count file was given to save search time and memory
1068             if ((namefile != "") || (countfile != "")) {
1069                 for (int j = 0; j < names.size(); j++) {
1070                     map<string, int>::iterator itNames = nameMap.find(names[j]);
1071                     if (itNames != nameMap.end()) { //add name because its a unique one
1072                         seq2Bin[names[j]] = i;
1073                     }
1074                 }
1075             }else { //map everyone
1076                 for (int j = 0; j < names.size(); j++) { seq2Bin[names[j]] = i;  }
1077             }
1078
1079             vector<unsigned int> temp; temp.resize(5, 0); //A,T,G,C,Gap
1080             vector< vector<unsigned int> > temp2;
1081             for (int j = 0; j < length; j++) {
1082                 temp2.push_back(temp);
1083             }
1084             counts.push_back(temp2);
1085         }
1086
1087         return 0;
1088     }
1089         catch(exception& e) {
1090                 m->errorOut(e, "PrimerDesignCommand", "initializeCounts");
1091                 exit(1);
1092         }
1093 }
1094 //**********************************************************************************************************************
1095 map<string, int> PrimerDesignCommand::readCount(unsigned long int& numSeqs){
1096         try {
1097         map<string, int> nameMap;
1098
1099         CountTable ct;
1100         ct.readTable(countfile, false, false);
1101         vector<string> namesOfSeqs = ct.getNamesOfSeqs();
1102         numSeqs = ct.getNumUniqueSeqs();
1103
1104         for (int i = 0; i < namesOfSeqs.size(); i++) {
1105             if (m->control_pressed) { break; }
1106
1107             nameMap[namesOfSeqs[i]] = ct.getNumSeqs(namesOfSeqs[i]);
1108         }
1109
1110         return nameMap;
1111     }
1112         catch(exception& e) {
1113                 m->errorOut(e, "PrimerDesignCommand", "readCount");
1114                 exit(1);
1115         }
1116 }
1117 //**********************************************************************************************************************
1118 int PrimerDesignCommand::getListVector(){
1119         try {
1120                 InputData input(listfile, "list");
1121                 list = input.getListVector();
1122                 string lastLabel = list->getLabel();
1123
1124                 if (label == "") { label = lastLabel;  return 0; }
1125
1126                 //if the users enters label "0.06" and there is no "0.06" in their file use the next lowest label.
1127                 set<string> labels; labels.insert(label);
1128                 set<string> processedLabels;
1129                 set<string> userLabels = labels;
1130
1131                 //as long as you are not at the end of the file or done wih the lines you want
1132                 while((list != NULL) && (userLabels.size() != 0)) {
1133                         if (m->control_pressed) {  return 0;  }
1134
1135                         if(labels.count(list->getLabel()) == 1){
1136                                 processedLabels.insert(list->getLabel());
1137                                 userLabels.erase(list->getLabel());
1138                                 break;
1139                         }
1140
1141                         if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1142                                 string saveLabel = list->getLabel();
1143
1144                                 delete list;
1145                                 list = input.getListVector(lastLabel);
1146
1147                                 processedLabels.insert(list->getLabel());
1148                                 userLabels.erase(list->getLabel());
1149
1150                                 //restore real lastlabel to save below
1151                                 list->setLabel(saveLabel);
1152                                 break;
1153                         }
1154
1155                         lastLabel = list->getLabel();
1156
1157                         //get next line to process
1158                         //prevent memory leak
1159                         delete list;
1160                         list = input.getListVector();
1161                 }
1162
1163
1164                 if (m->control_pressed) {  return 0;  }
1165
1166                 //output error messages about any remaining user labels
1167                 set<string>::iterator it;
1168                 bool needToRun = false;
1169                 for (it = userLabels.begin(); it != userLabels.end(); it++) {
1170                         m->mothurOut("Your file does not include the label " + *it);
1171                         if (processedLabels.count(lastLabel) != 1) {
1172                                 m->mothurOut(". I will use " + lastLabel + "."); m->mothurOutEndLine();
1173                                 needToRun = true;
1174                         }else {
1175                                 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1176                         }
1177                 }
1178
1179                 //run last label if you need to
1180                 if (needToRun == true)  {
1181                         delete list;
1182                         list = input.getListVector(lastLabel);
1183                 }
1184
1185                 return 0;
1186         }
1187         catch(exception& e) {
1188                 m->errorOut(e, "PrimerDesignCommand", "getListVector");
1189                 exit(1);
1190         }
1191 }
1192 //********************************************************************/
1193 /* A = adenine
1194  * C = cytosine
1195  * G = guanine
1196  * T = thymine
1197  * R = G A (purine)
1198  * Y = T C (pyrimidine)
1199  * K = G T (keto)
1200  * M = A C (amino)
1201  * S = G C (strong bonds)
1202  * W = A T (weak bonds)
1203  * B = G T C (all but A)
1204  * D = G A T (all but C)
1205  * H = A C T (all but G)
1206  * V = G C A (all but T)
1207  * N = A G C T (any) */
1208 int PrimerDesignCommand::countDiffs(string oligo, string seq){
1209         try {
1210
1211                 int length = oligo.length();
1212                 int countDiffs = 0;
1213
1214                 for(int i=0;i<length;i++){
1215
1216                         oligo[i] = toupper(oligo[i]);
1217             seq[i] = toupper(seq[i]);
1218
1219                         if(oligo[i] != seq[i]){
1220                 if(oligo[i] == 'A' && (seq[i] != 'A' && seq[i] != 'M' && seq[i] != 'R' && seq[i] != 'W' && seq[i] != 'D' && seq[i] != 'H' && seq[i] != 'V'))       {    countDiffs++;   }
1221                 else if(oligo[i] == 'C' && (seq[i] != 'C' && seq[i] != 'Y' && seq[i] != 'M' && seq[i] != 'S' && seq[i] != 'B' && seq[i] != 'H' && seq[i] != 'V'))       {       countDiffs++;   }
1222                 else if(oligo[i] == 'G' && (seq[i] != 'G' && seq[i] != 'R' && seq[i] != 'K' && seq[i] != 'S' && seq[i] != 'B' && seq[i] != 'D' && seq[i] != 'V'))       {       countDiffs++;   }
1223                 else if(oligo[i] == 'T' && (seq[i] != 'T' && seq[i] != 'Y' && seq[i] != 'K' && seq[i] != 'W' && seq[i] != 'B' && seq[i] != 'D' && seq[i] != 'H'))       {       countDiffs++;   }
1224                 else if((oligo[i] == '.' || oligo[i] == '-'))           {       countDiffs++;   }
1225                                 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N'))                         {      countDiffs++;   }
1226                                 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G'))                        {   countDiffs++;   }
1227                                 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T'))                        {   countDiffs++;   }
1228                                 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A'))                        {   countDiffs++;   }
1229                                 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G'))                        {   countDiffs++;   }
1230                                 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A'))                        {   countDiffs++;   }
1231                                 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G'))                        {   countDiffs++;   }
1232                                 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G'))       {   countDiffs++;   }
1233                                 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G'))       {   countDiffs++;   }
1234                                 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C'))       {   countDiffs++;   }
1235                                 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G'))       {   countDiffs++;   }
1236             }
1237
1238                 }
1239
1240                 return countDiffs;
1241         }
1242         catch(exception& e) {
1243                 m->errorOut(e, "PrimerDesignCommand", "countDiffs");
1244                 exit(1);
1245         }
1246 }
1247 //**********************************************************************************************************************
1248 int PrimerDesignCommand::findIndex(string binLabel, vector<string> binLabels){
1249         try {
1250         int index = -1;
1251         for (int i = 0; i < binLabels.size(); i++){
1252             if (m->control_pressed) { return index; }
1253             if (m->isLabelEquivalent(binLabel, binLabels[i])) { index = i; break; }
1254         }
1255         return index;
1256     }
1257         catch(exception& e) {
1258                 m->errorOut(e, "PrimerDesignCommand", "findIndex");
1259                 exit(1);
1260         }
1261 }
1262 //**********************************************************************************************************************
1263
1264
1265