2 // primerdesigncommand.cpp
5 // Created by Sarah Westcott on 1/18/13.
6 // Copyright (c) 2013 Schloss Lab. All rights reserved.
9 #include "primerdesigncommand.h"
11 //**********************************************************************************************************************
12 vector<string> PrimerDesignCommand::setParameters(){
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);
29 vector<string> myArray;
30 for (int i = 0; i < parameters.size(); i++) { myArray.push_back(parameters[i].name); }
34 m->errorOut(e, "PrimerDesignCommand", "setParameters");
38 //**********************************************************************************************************************
39 string PrimerDesignCommand::getHelpString(){
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";
61 m->errorOut(e, "PrimerDesignCommand", "getHelpString");
65 //**********************************************************************************************************************
66 string PrimerDesignCommand::getOutputPattern(string type) {
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; }
78 m->errorOut(e, "PrimerDesignCommand", "getOutputPattern");
82 //**********************************************************************************************************************
83 PrimerDesignCommand::PrimerDesignCommand(){
85 abort = true; calledHelp = true;
87 vector<string> tempOutNames;
88 outputTypes["summary"] = tempOutNames;
89 outputTypes["fasta"] = tempOutNames;
90 outputTypes["list"] = tempOutNames;
94 m->errorOut(e, "PrimerDesignCommand", "PrimerDesignCommand");
98 //**********************************************************************************************************************
99 PrimerDesignCommand::PrimerDesignCommand(string option) {
101 abort = false; calledHelp = false;
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;}
108 //valid paramters for this command
109 vector<string> myArray = setParameters();
111 OptionParser parser(option);
112 map<string,string> parameters = parser.getParameters();
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; }
121 vector<string> tempOutNames;
122 outputTypes["summary"] = tempOutNames;
123 outputTypes["fasta"] = tempOutNames;
124 outputTypes["list"] = tempOutNames;
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 = ""; }
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; }
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; }
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; }
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; }
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); }
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); }
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); }
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); }
194 if ((namefile != "") && (countfile != "")) {
195 m->mothurOut("[ERROR]: you may only use one of the following: name or count."); m->mothurOutEndLine(); abort = true;
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
204 string temp = validParameter.validFile(parameters, "cutoff", false); if (temp == "not found") { temp = "100"; }
205 m->mothurConvert(temp, cutoff);
207 temp = validParameter.validFile(parameters, "pdiffs", false); if (temp == "not found") { temp = "0"; }
208 m->mothurConvert(temp, pdiffs);
210 temp = validParameter.validFile(parameters, "length", false); if (temp == "not found") { temp = "18"; }
211 m->mothurConvert(temp, length);
213 temp = validParameter.validFile(parameters, "mintm", false); if (temp == "not found") { temp = "-1"; }
214 m->mothurConvert(temp, minTM);
216 temp = validParameter.validFile(parameters, "maxtm", false); if (temp == "not found") { temp = "-1"; }
217 m->mothurConvert(temp, maxTM);
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; }
222 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = m->getProcessors(); }
223 m->setProcessors(temp);
224 m->mothurConvert(temp, processors);
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=""; }
229 if (countfile == "") {
230 if (namefile == "") {
231 vector<string> files; files.push_back(fastafile);
232 parser.getNameFile(files);
237 catch(exception& e) {
238 m->errorOut(e, "PrimerDesignCommand", "PrimerDesignCommand");
242 //**********************************************************************************************************************
243 int PrimerDesignCommand::execute(){
246 if (abort == true) { if (calledHelp) { return 0; } return 2; }
248 int start = time(NULL);
249 //////////////////////////////////////////////////////////////////////////////
250 // get file inputs //
251 //////////////////////////////////////////////////////////////////////////////
253 //reads list file and selects the label the users specified or the first label
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; }
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.
263 if (namefile != "") { nameMap = m->readNames(namefile, numSeqs); }
264 else if (countfile != "") { nameMap = readCount(numSeqs); }
265 else { numSeqs = list->getNumSeqs(); }
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");
273 m->control_pressed = true;
276 if (m->control_pressed) { delete list; return 0; }
278 //////////////////////////////////////////////////////////////////////////////
280 //////////////////////////////////////////////////////////////////////////////
281 m->mothurOut("\nFinding consensus sequences for each otu..."); cout.flush();
283 vector<Sequence> conSeqs = createProcessesConSeqs(nameMap, numSeqs);
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);
291 m->openOutputFile(consFastaFile, out);
292 for (int i = 0; i < conSeqs.size(); i++) { conSeqs[i].printSequence(out); }
295 m->mothurOut("Done.\n\n");
297 set<string> primers = getPrimer(conSeqs[binIndex]);
299 if (m->control_pressed) { delete list; for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
301 string consSummaryFile = getOutputFileName("summary", variables);
302 outputNames.push_back(consSummaryFile); outputTypes["summary"].push_back(consSummaryFile);
304 m->openOutputFile(consSummaryFile, outSum);
306 outSum << "PrimerOtu: " << otulabel << " Members: " << list->get(binIndex) << endl << "Primers\tminTm\tmaxTm" << endl;
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();) {
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;
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;
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;
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;
336 }else { primers.erase(it++); } //erase because it didn't qualify
339 outSum << "\nOTUNumber\tPrimer\tStart\tEnd\tLength\tMismatches\tminTm\tmaxTm\n";
342 //check each otu's conseq for each primer in otunumber
343 set<int> otuToRemove = createProcesses(consSummaryFile, minTms, maxTms, primers, conSeqs, binIndex);
345 if (m->control_pressed) { delete list; for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
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);
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; }
360 if (otuToRemove.count(j) == 0) {
361 string bin = list->get(j);
362 if (bin != "") { outListTemp << bin << '\t'; headers += binLabels[j] + '\t'; }
369 m->openOutputFile(newListFile, outList);
370 outList << headers << endl;
372 m->appendFiles(newListFile+".temp", newListFile);
373 m->mothurRemove(newListFile+".temp");
374 outputNames.push_back(newListFile); outputTypes["list"].push_back(newListFile);
376 if (m->control_pressed) { delete list; for (int i = 0; i < outputNames.size(); i++) { m->mothurRemove(outputNames[i]); } return 0; }
380 m->mothurOut("It took " + toString(time(NULL) - start) + " secs to process " + toString(list->getNumBins()) + " OTUs.\n");
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();
392 catch(exception& e) {
393 m->errorOut(e, "PrimerDesignCommand", "execute");
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
407 * Y = T C (pyrimidine)
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) */
418 int PrimerDesignCommand::findMeltingPoint(string primer, double& minTm, double& maxTm){
420 string minTmprimer = primer;
421 string maxTmprimer = primer;
423 //find minimum Tm string substituting for degenerate bases
424 for (int i = 0; i < minTmprimer.length(); i++) {
425 minTmprimer[i] = toupper(minTmprimer[i]);
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'; }
440 //find maximum Tm string substituting for degenerate bases
441 for (int i = 0; i < maxTmprimer.length(); i++) {
442 maxTmprimer[i] = toupper(maxTmprimer[i]);
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'; }
458 for (int i = 0; i < minTmprimer.length(); i++) {
459 if (minTmprimer[i] == 'G') { numGC++; }
460 else if (minTmprimer[i] == 'C') { numGC++; }
463 minTm = 64.9 + 41 * (numGC - 16.4) / (double) minTmprimer.length();
466 for (int i = 0; i < maxTmprimer.length(); i++) {
467 if (maxTmprimer[i] == 'G') { numGC++; }
468 else if (maxTmprimer[i] == 'C') { numGC++; }
471 maxTm = 64.9 + 41 * (numGC - 16.4) / (double) maxTmprimer.length();
475 catch(exception& e) {
476 m->errorOut(e, "PrimerDesignCommand", "findMeltingPoint");
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){
484 bool foundAtLeastOne = false; //innocent til proven guilty
486 //look for exact match
487 if(rawSequence.length() < primer.length()) { return false; }
490 for (int j = 0; j < rawSequence.length()-length; j++){
492 if (m->control_pressed) { return foundAtLeastOne; }
494 string rawChunk = rawSequence.substr(j, length);
496 int numDiff = countDiffs(primer, rawChunk);
498 if(numDiff <= pdiffs){
499 primerStart.push_back(j);
500 primerEnd.push_back(j+length);
501 mismatches.push_back(numDiff);
502 foundAtLeastOne = true;
506 return foundAtLeastOne;
509 catch(exception& e) {
510 m->errorOut(e, "PrimerDesignCommand", "findPrimer");
514 //********************************************************************/
515 //find all primers for the given sequence
516 set<string> PrimerDesignCommand::getPrimer(Sequence primerSeq){
520 string rawSequence = primerSeq.getUnaligned();
522 for (int j = 0; j < rawSequence.length()-length; j++){
523 if (m->control_pressed) { break; }
525 string primer = rawSequence.substr(j, length);
526 primers.insert(primer);
531 catch(exception& e) {
532 m->errorOut(e, "PrimerDesignCommand", "getPrimer");
536 /**************************************************************************************************/
537 set<int> PrimerDesignCommand::createProcesses(string newSummaryFile, vector<double>& minTms, vector<double>& maxTms, set<string>& primers, vector<Sequence>& conSeqs, int binIndex) {
540 vector<int> processIDS;
542 set<int> otusToRemove;
543 int numBinsProcessed = 0;
546 int numBins = conSeqs.size();
547 if (numBins < processors) { processors = numBins; }
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));
559 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
561 //loop through and create all the processes you want
562 while (process != processors) {
566 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
569 //clear old file because we append in driver
570 m->mothurRemove(newSummaryFile + toString(getpid()) + ".temp");
572 otusToRemove = driver(newSummaryFile + toString(getpid()) + ".temp", minTms, maxTms, primers, conSeqs, lines[process].start, lines[process].end, numBinsProcessed, binIndex);
574 string tempFile = toString(getpid()) + ".otus2Remove.temp";
576 m->openOutputFile(tempFile, outTemp);
578 outTemp << numBinsProcessed << endl;
579 outTemp << otusToRemove.size() << endl;
580 for (set<int>::iterator it = otusToRemove.begin(); it != otusToRemove.end(); it++) { outTemp << *it << endl; }
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); }
592 otusToRemove = driver(newSummaryFile, minTms, maxTms, primers, conSeqs, lines[0].start, lines[0].end, numBinsProcessed, binIndex);
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];
600 for (int i = 0; i < processIDS.size(); i++) {
601 string tempFile = toString(processIDS[i]) + ".otus2Remove.temp";
603 m->openInputFile(tempFile, intemp);
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++) {
611 intemp >> otu; m->gobble(intemp);
612 otusToRemove.insert(otu);
615 m->mothurRemove(tempFile);
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 //////////////////////////////////////////////////////////////////////////////////////////////////////
626 vector<primerDesignData*> pDataArray;
627 DWORD dwThreadIdArray[processors-1];
628 HANDLE hThreadArray[processors-1];
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);
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);
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]);
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);
649 //Wait until all threads have terminated.
650 WaitForMultipleObjects(processors-1, hThreadArray, TRUE, INFINITE);
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);
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];
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"));
674 catch(exception& e) {
675 m->errorOut(e, "PrimerDesignCommand", "createProcesses");
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){
682 set<int> otuToRemove;
685 m->openOutputFileAppend(summaryFileName, outSum);
687 for (int i = start; i < end; i++) {
689 if (m->control_pressed) { break; }
691 if (i != (binIndex)) {
693 for (set<string>::iterator it = primers.begin(); it != primers.end(); it++) {
694 vector<int> primerStarts;
695 vector<int> primerEnds;
696 vector<int> mismatches;
698 bool found = findPrimer(conSeqs[i].getUnaligned(), (*it), primerStarts, primerEnds, mismatches);
700 //if we found it report to the table
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;
705 otuToRemove.insert(i);
717 catch(exception& e) {
718 m->errorOut(e, "PrimerDesignCommand", "driver");
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){
725 vector< vector< vector<unsigned int> > > counts;
726 map<string, int> seq2Bin;
730 m->openInputFile(fastafile, in);
738 if (m->control_pressed) { in.close(); return counts; }
740 Sequence seq(in); m->gobble(in);
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;
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; }
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; }
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;
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;
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; }
787 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
788 unsigned long long pos = in.tellg();
789 if ((pos == -1) || (pos >= end)) { break; }
791 if (in.eof()) { break; }
799 catch(exception& e) {
800 m->errorOut(e, "PrimerDesignCommand", "driverGetCounts");
804 /**************************************************************************************************/
805 vector<Sequence> PrimerDesignCommand::createProcessesConSeqs(map<string, int>& nameMap, unsigned long int& numSeqs) {
807 vector< vector< vector<unsigned int> > > counts;
808 vector<unsigned int> otuCounts;
809 vector<int> processIDS;
811 unsigned long int fastaCount = 0;
813 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
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)])); }
820 //loop through and create all the processes you want
821 while (process != processors) {
825 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
828 counts = driverGetCounts(nameMap, fastaCount, otuCounts, lines[process].start, lines[process].end);
830 string tempFile = toString(getpid()) + ".cons_counts.temp";
832 m->openOutputFile(tempFile, outTemp);
834 outTemp << fastaCount << endl;
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'; }
845 outTemp << otuCounts.size() << endl;
846 for (int i = 0; i < otuCounts.size(); i++) { outTemp << otuCounts[i] << '\t'; }
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); }
859 counts = driverGetCounts(nameMap, fastaCount, otuCounts, lines[0].start, lines[0].end);
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];
867 for (int i = 0; i < processIDS.size(); i++) {
868 string tempFile = toString(processIDS[i]) + ".cons_counts.temp";
870 m->openInputFile(tempFile, intemp);
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; }
878 for (int k = 0; k < num; k++) {
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; }
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; }
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;
896 m->mothurRemove(tempFile);
901 unsigned long long start = 0;
902 unsigned long long end = 1000;
903 counts = driverGetCounts(nameMap, fastaCount, otuCounts, start, end);
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;
912 vector<Sequence> conSeqs;
914 if (m->control_pressed) { return conSeqs; }
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; }
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"; }
927 otuLabel += sbinNumber;
930 for (int j = 0; j < counts[i].size(); j++) {
931 cons += getBase(counts[i][j], otuCounts[i]);
933 Sequence consSeq(otuLabel, cons);
934 conSeqs.push_back(consSeq);
937 if (m->control_pressed) { conSeqs.clear(); return conSeqs; }
943 catch(exception& e) {
944 m->errorOut(e, "PrimerDesignCommand", "createProcessesConSeqs");
948 //***************************************************************************************************************
950 char PrimerDesignCommand::getBase(vector<unsigned int> counts, int size){ //A,T,G,C,Gap
957 * Y = T C (pyrimidine)
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) */
970 //zero out counts that don't make the cutoff
971 float percentage = (100.0 - cutoff) / 100.0;
973 for (int i = 0; i < counts.size(); i++) {
974 float countPercentage = counts[i] / (float) size;
975 if (countPercentage < percentage) { counts[i] = 0; }
979 if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] != 0)) { conBase = 'n'; }
981 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] == 0)) { conBase = 'N'; }
983 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] != 0)) { conBase = 'v'; }
985 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] == 0)) { conBase = 'V'; }
987 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] != 0)) { conBase = 'h'; }
989 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] == 0)) { conBase = 'H'; }
991 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] != 0)) { conBase = 'd'; }
993 else if ((counts[0] != 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] == 0)) { conBase = 'D'; }
995 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] != 0) && (counts[3] != 0) && (counts[4] != 0)) { conBase = 'b'; }
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'; }
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'; }
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'; }
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'; }
1023 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] != 0)) { conBase = 'a'; }
1025 else if ((counts[0] != 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] == 0)) { conBase = 'A'; }
1027 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] != 0)) { conBase = 't'; }
1029 else if ((counts[0] == 0) && (counts[1] != 0) && (counts[2] == 0) && (counts[3] == 0) && (counts[4] == 0)) { conBase = 'T'; }
1031 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] != 0)) { conBase = 'g'; }
1033 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] != 0) && (counts[3] == 0) && (counts[4] == 0)) { conBase = 'G'; }
1035 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] != 0)) { conBase = 'c'; }
1037 else if ((counts[0] == 0) && (counts[1] == 0) && (counts[2] == 0) && (counts[3] != 0) && (counts[4] == 0)) { conBase = 'C'; }
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(); }
1047 catch(exception& e) {
1048 m->errorOut(e, "PrimerDesignCommand", "getBase");
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){
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);
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;
1075 }else { //map everyone
1076 for (int j = 0; j < names.size(); j++) { seq2Bin[names[j]] = i; }
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);
1084 counts.push_back(temp2);
1089 catch(exception& e) {
1090 m->errorOut(e, "PrimerDesignCommand", "initializeCounts");
1094 //**********************************************************************************************************************
1095 map<string, int> PrimerDesignCommand::readCount(unsigned long int& numSeqs){
1097 map<string, int> nameMap;
1100 ct.readTable(countfile, false, false);
1101 vector<string> namesOfSeqs = ct.getNamesOfSeqs();
1102 numSeqs = ct.getNumUniqueSeqs();
1104 for (int i = 0; i < namesOfSeqs.size(); i++) {
1105 if (m->control_pressed) { break; }
1107 nameMap[namesOfSeqs[i]] = ct.getNumSeqs(namesOfSeqs[i]);
1112 catch(exception& e) {
1113 m->errorOut(e, "PrimerDesignCommand", "readCount");
1117 //**********************************************************************************************************************
1118 int PrimerDesignCommand::getListVector(){
1120 InputData input(listfile, "list");
1121 list = input.getListVector();
1122 string lastLabel = list->getLabel();
1124 if (label == "") { label = lastLabel; return 0; }
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;
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; }
1135 if(labels.count(list->getLabel()) == 1){
1136 processedLabels.insert(list->getLabel());
1137 userLabels.erase(list->getLabel());
1141 if ((m->anyLabelsToProcess(list->getLabel(), userLabels, "") == true) && (processedLabels.count(lastLabel) != 1)) {
1142 string saveLabel = list->getLabel();
1145 list = input.getListVector(lastLabel);
1147 processedLabels.insert(list->getLabel());
1148 userLabels.erase(list->getLabel());
1150 //restore real lastlabel to save below
1151 list->setLabel(saveLabel);
1155 lastLabel = list->getLabel();
1157 //get next line to process
1158 //prevent memory leak
1160 list = input.getListVector();
1164 if (m->control_pressed) { return 0; }
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();
1175 m->mothurOut(". Please refer to " + lastLabel + "."); m->mothurOutEndLine();
1179 //run last label if you need to
1180 if (needToRun == true) {
1182 list = input.getListVector(lastLabel);
1187 catch(exception& e) {
1188 m->errorOut(e, "PrimerDesignCommand", "getListVector");
1192 //********************************************************************/
1198 * Y = T C (pyrimidine)
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){
1211 int length = oligo.length();
1214 for(int i=0;i<length;i++){
1216 oligo[i] = toupper(oligo[i]);
1217 seq[i] = toupper(seq[i]);
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++; }
1242 catch(exception& e) {
1243 m->errorOut(e, "PrimerDesignCommand", "countDiffs");
1247 //**********************************************************************************************************************
1248 int PrimerDesignCommand::findIndex(string binLabel, vector<string> binLabels){
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; }
1257 catch(exception& e) {
1258 m->errorOut(e, "PrimerDesignCommand", "findIndex");
1262 //**********************************************************************************************************************