5 * Created by Pat Schloss on 12/22/10.
6 * Copyright 2010 Schloss Lab. All rights reserved.
10 #include "trimflowscommand.h"
11 #include "needlemanoverlap.hpp"
13 //**********************************************************************************************************************
15 vector<string> TrimFlowsCommand::getValidParameters(){
17 string Array[] = {"flow", "maxflows", "minflows",
18 "fasta", "minlength", "maxlength", "maxhomop", "signal", "noise"
19 "oligos", "pdiffs", "bdiffs", "tdiffs", "order",
20 "allfiles", "processors",
21 "outputdir","inputdir"
24 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
28 m->errorOut(e, "TrimFlowsCommand", "getValidParameters");
33 //**********************************************************************************************************************
35 vector<string> TrimFlowsCommand::getRequiredParameters(){
37 string Array[] = {"flow"};
38 vector<string> myArray (Array, Array+(sizeof(Array)/sizeof(string)));
42 m->errorOut(e, "TrimFlowsCommand", "getRequiredParameters");
47 //**********************************************************************************************************************
49 vector<string> TrimFlowsCommand::getRequiredFiles(){
51 vector<string> myArray;
55 m->errorOut(e, "TrimFlowsCommand", "getRequiredFiles");
60 //**********************************************************************************************************************
62 TrimFlowsCommand::TrimFlowsCommand(){
64 abort = true; calledHelp = true;
65 vector<string> tempOutNames;
66 outputTypes["flow"] = tempOutNames;
67 outputTypes["fasta"] = tempOutNames;
70 m->errorOut(e, "TrimFlowsCommand", "TrimFlowsCommand");
75 //***************************************************************************************************************
77 TrimFlowsCommand::~TrimFlowsCommand(){ /* do nothing */ }
79 //***************************************************************************************************************
81 void TrimFlowsCommand::help(){
83 m->mothurOut("The trim.flows command reads a flowgram file and creates .....\n");
84 m->mothurOut("Note: No spaces between parameter labels (i.e. fasta), '=' and parameters (i.e.yourFasta).\n");
85 m->mothurOut("For more details please check out the wiki http://www.mothur.org/wiki/Trim.flows.\n\n");
89 m->errorOut(e, "TrimFlowsCommand", "help");
94 //**********************************************************************************************************************
96 TrimFlowsCommand::TrimFlowsCommand(string option) {
99 abort = false; calledHelp = false;
102 //allow user to run help
103 if(option == "help") { help(); abort = true; calledHelp = true; }
106 //valid paramters for this command
107 string AlignArray[] = {"flow", "maxflows", "minflows",
108 "fasta", "minlength", "maxlength", "maxhomop", "signal", "noise",
109 "oligos", "pdiffs", "bdiffs", "tdiffs", "order",
110 "allfiles", "processors",
113 "outputdir","inputdir"
117 vector<string> myArray (AlignArray, AlignArray+(sizeof(AlignArray)/sizeof(string)));
119 OptionParser parser(option);
120 map<string,string> parameters = parser.getParameters();
122 ValidParameters validParameter;
123 map<string,string>::iterator it;
125 //check to make sure all parameters are valid for command
126 for (it = parameters.begin(); it != parameters.end(); it++) {
127 if (validParameter.isValidParameter(it->first, myArray, it->second) != true) { abort = true; }
130 //initialize outputTypes
131 vector<string> tempOutNames;
132 outputTypes["flow"] = tempOutNames;
133 outputTypes["fasta"] = tempOutNames;
135 //if the user changes the input directory command factory will send this info to us in the output parameter
136 string inputDir = validParameter.validFile(parameters, "inputdir", false);
137 if (inputDir == "not found"){ inputDir = ""; }
140 it = parameters.find("flow");
141 //user has given a template file
142 if(it != parameters.end()){
143 path = m->hasPath(it->second);
144 //if the user has not given a path then, add inputdir. else leave path alone.
145 if (path == "") { parameters["flow"] = inputDir + it->second; }
148 it = parameters.find("oligos");
149 //user has given a template file
150 if(it != parameters.end()){
151 path = m->hasPath(it->second);
152 //if the user has not given a path then, add inputdir. else leave path alone.
153 if (path == "") { parameters["oligos"] = inputDir + it->second; }
157 // it = parameters.find("group");
158 // //user has given a template file
159 // if(it != parameters.end()){
160 // path = m->hasPath(it->second);
161 // //if the user has not given a path then, add inputdir. else leave path alone.
162 // if (path == "") { parameters["group"] = inputDir + it->second; }
167 //check for required parameters
168 flowFileName = validParameter.validFile(parameters, "flow", true);
169 if (flowFileName == "not found") { m->mothurOut("flow is a required parameter for the trim.flows command."); m->mothurOutEndLine(); abort = true; }
170 else if (flowFileName == "not open") { abort = true; }
172 //if the user changes the output directory command factory will send this info to us in the output parameter
173 outputDir = validParameter.validFile(parameters, "outputdir", false); if (outputDir == "not found"){
175 outputDir += m->hasPath(flowFileName); //if user entered a file with a path then preserve it
179 //check for optional parameter and set defaults
180 // ...at some point should added some additional type checking...
183 temp = validParameter.validFile(parameters, "minflows", false); if (temp == "not found") { temp = "360"; }
184 convert(temp, minFlows);
186 temp = validParameter.validFile(parameters, "maxflows", false); if (temp == "not found") { temp = "720"; }
187 convert(temp, maxFlows);
190 temp = validParameter.validFile(parameters, "oligos", true);
191 if (temp == "not found") { oligoFileName = ""; }
192 else if(temp == "not open") { abort = true; }
193 else { oligoFileName = temp; }
195 temp = validParameter.validFile(parameters, "fasta", false); if (temp == "not found"){ fasta = 0; }
196 else if(m->isTrue(temp)) { fasta = 1; }
198 temp = validParameter.validFile(parameters, "maxhomop", false); if (temp == "not found"){ temp = "9"; }
199 convert(temp, maxHomoP);
201 temp = validParameter.validFile(parameters, "signal", false); if (temp == "not found"){ temp = "0.50"; }
202 convert(temp, signal);
204 temp = validParameter.validFile(parameters, "noise", false); if (temp == "not found"){ temp = "0.70"; }
205 convert(temp, noise);
207 temp = validParameter.validFile(parameters, "minlength", false); if (temp == "not found"){ temp = "0"; }
208 convert(temp, minLength);
210 temp = validParameter.validFile(parameters, "maxlength", false); if (temp == "not found"){ temp = "0"; }
211 convert(temp, maxLength);
213 temp = validParameter.validFile(parameters, "bdiffs", false); if (temp == "not found"){ temp = "0"; }
214 convert(temp, bdiffs);
216 temp = validParameter.validFile(parameters, "pdiffs", false); if (temp == "not found"){ temp = "0"; }
217 convert(temp, pdiffs);
219 temp = validParameter.validFile(parameters, "tdiffs", false);
220 if (temp == "not found"){ int tempTotal = pdiffs + bdiffs; temp = toString(tempTotal); }
221 convert(temp, tdiffs);
222 if(tdiffs == 0){ tdiffs = bdiffs + pdiffs; }
224 temp = validParameter.validFile(parameters, "allfiles", false); if (temp == "not found"){ temp = "T"; }
225 allFiles = m->isTrue(temp);
227 temp = validParameter.validFile(parameters, "processors", false); if (temp == "not found"){ temp = "1"; }
228 convert(temp, processors);
230 flowOrder = validParameter.validFile(parameters, "order", false);
231 if (flowOrder == "not found"){ flowOrder = "TACG"; }
232 else if(flowOrder.length() != 4){
233 m->mothurOut("The value of the order option must be four bases long\n");
236 if(oligoFileName == ""){ allFiles = 0; }
243 catch(exception& e) {
244 m->errorOut(e, "TrimFlowsCommand", "TrimSeqsCommand");
249 //***************************************************************************************************************
251 int TrimFlowsCommand::execute(){
254 if (abort == true) { if (calledHelp) { return 0; } return 2; }
256 string trimFlowFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "trim.flow";
257 outputNames.push_back(trimFlowFileName); outputTypes["flow"].push_back(trimFlowFileName);
259 string scrapFlowFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "scrap.flow";
260 outputNames.push_back(scrapFlowFileName); outputTypes["flow"].push_back(scrapFlowFileName);
262 string fastaFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.fasta";
264 outputNames.push_back(fastaFileName); outputTypes["fasta"].push_back(fastaFileName);
267 vector<unsigned long int> flowFilePos = getFlowFileBreaks();
268 for (int i = 0; i < (flowFilePos.size()-1); i++) {
269 lines.push_back(new linePair(flowFilePos[i], flowFilePos[(i+1)]));
272 vector<vector<string> > barcodePrimerComboFileNames;
273 if(oligoFileName != ""){
274 getOligos(barcodePrimerComboFileNames);
277 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
279 driverCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames, lines[0]);
281 createProcessesCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames);
284 driverCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames, lines[0]);
287 if (m->control_pressed) { return 0; }
289 string flowFilesFileName;
294 flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files";
295 m->openOutputFile(flowFilesFileName, output);
297 for(int i=0;i<barcodePrimerComboFileNames.size();i++){
298 for(int j=0;j<barcodePrimerComboFileNames[0].size();j++){
301 unsigned long int size;
303 //get num bytes in file
304 pFile = fopen (barcodePrimerComboFileNames[i][j].c_str(),"rb");
305 if (pFile==NULL) perror ("Error opening file");
307 fseek (pFile, 0, SEEK_END);
313 remove(barcodePrimerComboFileNames[i][j].c_str());
316 output << barcodePrimerComboFileNames[i][j] << endl;
323 flowFilesFileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + "flow.files";
324 m->openOutputFile(flowFilesFileName, output);
326 output << trimFlowFileName << endl;
330 outputTypes["flow.files"].push_back(flowFilesFileName);
331 outputNames.push_back(flowFileName);
333 //set fasta file as new current fastafile
335 itTypes = outputTypes.find("fasta");
336 if (itTypes != outputTypes.end()) {
337 if ((itTypes->second).size() != 0) { current = (itTypes->second)[0]; m->setFastaFile(current); }
340 m->mothurOutEndLine();
341 m->mothurOut("Output File Names: "); m->mothurOutEndLine();
342 for (int i = 0; i < outputNames.size(); i++) { m->mothurOut(outputNames[i]); m->mothurOutEndLine(); }
343 m->mothurOutEndLine();
347 catch(exception& e) {
348 m->errorOut(e, "TrimSeqsCommand", "execute");
353 //***************************************************************************************************************
355 int TrimFlowsCommand::driverCreateTrim(string flowFileName, string trimFlowFileName, string scrapFlowFileName, string fastaFileName, vector<vector<string> > barcodePrimerComboFileNames, linePair* line){
359 ofstream trimFlowFile;
360 m->openOutputFile(trimFlowFileName, trimFlowFile);
361 trimFlowFile.setf(ios::fixed, ios::floatfield); trimFlowFile.setf(ios::showpoint);
363 ofstream scrapFlowFile;
364 m->openOutputFile(scrapFlowFileName, scrapFlowFile);
365 scrapFlowFile.setf(ios::fixed, ios::floatfield); scrapFlowFile.setf(ios::showpoint);
367 if(line->start == 4){
368 scrapFlowFile << maxFlows << endl;
369 trimFlowFile << maxFlows << endl;
371 for(int i=0;i<barcodePrimerComboFileNames.size();i++){
372 for(int j=0;j<barcodePrimerComboFileNames[0].size();j++){
373 // barcodePrimerComboFileNames[i][j] += toString(getpid()) + ".temp";
375 m->openOutputFile(barcodePrimerComboFileNames[i][j], temp);
376 temp << maxFlows << endl;
383 FlowData flowData(numFlows, signal, noise, maxHomoP, flowOrder);
386 if(fasta){ m->openOutputFile(fastaFileName, fastaFile); }
389 m->openInputFile(flowFileName, flowFile);
391 flowFile.seekg(line->start);
399 int currentSeqDiffs = 0;
400 string trashCode = "";
402 flowData.getNext(flowFile);
403 flowData.capFlows(maxFlows);
405 Sequence currSeq = flowData.getSequence();
406 if(!flowData.hasMinFlows(minFlows)){ //screen to see if sequence is of a minimum number of flows
411 if(minLength > 0 || maxLength > 0){ //screen to see if sequence is above and below a specific number of bases
412 int seqLength = currSeq.getNumBases();
413 if(seqLength < minLength || seqLength > maxLength){
420 int barcodeIndex = 0;
422 if(barcodes.size() != 0){
423 success = stripBarcode(currSeq, barcodeIndex);
424 if(success > bdiffs) { trashCode += 'b'; }
425 else{ currentSeqDiffs += success; }
428 if(numFPrimers != 0){
429 success = stripForward(currSeq, primerIndex);
430 if(success > pdiffs) { trashCode += 'f'; }
431 else{ currentSeqDiffs += success; }
434 if (currentSeqDiffs > tdiffs) { trashCode += 't'; }
436 if(numRPrimers != 0){
437 success = stripReverse(currSeq);
438 if(!success) { trashCode += 'r'; }
441 if(trashCode.length() == 0){
443 flowData.printFlows(trimFlowFile);
445 if(fasta) { currSeq.printSequence(fastaFile); }
449 m->openOutputFileAppend(barcodePrimerComboFileNames[barcodeIndex][primerIndex], output);
450 output.setf(ios::fixed, ios::floatfield); trimFlowFile.setf(ios::showpoint);
452 flowData.printFlows(output);
457 flowData.printFlows(scrapFlowFile, trashCode);
463 if((count) % 10000 == 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
465 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
466 unsigned long int pos = flowFile.tellg();
468 if ((pos == -1) || (pos >= line->end)) { break; }
470 if (flowFile.eof()) { break; }
475 if((count) % 10000 != 0){ m->mothurOut(toString(count)); m->mothurOutEndLine(); }
477 trimFlowFile.close();
478 scrapFlowFile.close();
480 if(fasta){ fastaFile.close(); }
484 catch(exception& e) {
485 m->errorOut(e, "TrimSeqsCommand", "driverCreateTrim");
490 //***************************************************************************************************************
492 void TrimFlowsCommand::getOligos(vector<vector<string> >& outFlowFileNames){
495 m->openInputFile(oligoFileName, oligosFile);
497 string type, oligo, group;
500 int indexBarcode = 0;
502 while(!oligosFile.eof()){
504 oligosFile >> type; m->gobble(oligosFile); //get the first column value of the row - is it a comment or a feature we are interested in?
506 if(type[0] == '#'){ //igore the line because there's a comment
507 while (!oligosFile.eof()) { char c = oligosFile.get(); if (c == 10 || c == 13){ break; } } // get rest of line if there's any crap there
509 else{ //there's a feature we're interested in
511 for(int i=0;i<type.length();i++){ type[i] = toupper(type[i]); } //make type case insensitive
513 oligosFile >> oligo; //get the DNA sequence for the feature
515 for(int i=0;i<oligo.length();i++){ //make type case insensitive and change any U's to T's
516 oligo[i] = toupper(oligo[i]);
517 if(oligo[i] == 'U') { oligo[i] = 'T'; }
520 if(type == "FORWARD"){ //if the feature is a forward primer...
523 while (!oligosFile.eof()) { // get rest of line in case there is a primer name = will have the name of the primer
524 char c = oligosFile.get();
525 if (c == 10 || c == 13){ break; }
526 else if (c == 32 || c == 9){;} //space or tab
530 //have we seen this primer already?
531 map<string, int>::iterator itPrimer = primers.find(oligo);
532 if (itPrimer != primers.end()) { m->mothurOut("primer " + oligo + " is in your oligos file already."); m->mothurOutEndLine(); }
534 primers[oligo]=indexPrimer; indexPrimer++;
535 primerNameVector.push_back(group);
538 else if(type == "REVERSE"){
539 Sequence oligoRC("reverse", oligo);
540 oligoRC.reverseComplement();
541 revPrimer.push_back(oligoRC.getUnaligned());
543 else if(type == "BARCODE"){
546 //check for repeat barcodes
547 map<string, int>::iterator itBar = barcodes.find(oligo);
548 if (itBar != barcodes.end()) { m->mothurOut("barcode " + oligo + " is in your oligos file already."); m->mothurOutEndLine(); }
550 barcodes[oligo]=indexBarcode; indexBarcode++;
551 barcodeNameVector.push_back(group);
554 m->mothurOut(type + " is not recognized as a valid type. Choices are forward, reverse, and barcode. Ignoring " + oligo + "."); m->mothurOutEndLine();
558 m->gobble(oligosFile);
562 if(barcodeNameVector.size() == 0 && primerNameVector[0] == ""){ allFiles = 0; }
564 //add in potential combos
565 if(barcodeNameVector.size() == 0){
567 barcodeNameVector.push_back("");
570 if(primerNameVector.size() == 0){
572 primerNameVector.push_back("");
576 outFlowFileNames.resize(barcodeNameVector.size());
577 for(int i=0;i<outFlowFileNames.size();i++){
578 outFlowFileNames[i].assign(primerNameVector.size(), "");
583 for(map<string, int>::iterator itBar = barcodes.begin();itBar != barcodes.end();itBar++){
584 for(map<string, int>::iterator itPrimer = primers.begin();itPrimer != primers.end(); itPrimer++){
586 string primerName = primerNameVector[itPrimer->second];
587 string barcodeName = barcodeNameVector[itBar->second];
589 string comboGroupName = "";
590 string fileName = "";
592 if(primerName == ""){
593 comboGroupName = barcodeNameVector[itBar->second];
594 fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + comboGroupName + ".flow";
597 if(barcodeName == ""){
598 comboGroupName = primerNameVector[itPrimer->second];
601 comboGroupName = barcodeNameVector[itBar->second] + "." + primerNameVector[itPrimer->second];
603 fileName = outputDir + m->getRootName(m->getSimpleName(flowFileName)) + comboGroupName + ".flow";
606 outputNames.push_back(fileName);
607 outputTypes["flow"].push_back(fileName);
608 outFlowFileNames[itBar->second][itPrimer->second] = fileName;
611 m->openOutputFile(fileName, temp);
617 numFPrimers = primers.size();
618 numRPrimers = revPrimer.size();
621 catch(exception& e) {
622 m->errorOut(e, "TrimSeqsCommand", "getOligos");
627 //***************************************************************************************************************
629 int TrimFlowsCommand::stripBarcode(Sequence& seq, int& group){
632 string rawSequence = seq.getUnaligned();
633 int success = bdiffs + 1; //guilty until proven innocent
635 //can you find the barcode
636 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
637 string oligo = it->first;
638 if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length
639 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
643 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
645 seq.setUnaligned(rawSequence.substr(oligo.length()));
652 //if you found the barcode or if you don't want to allow for diffs
653 if ((bdiffs == 0) || (success == 0)) { return success; }
655 else { //try aligning and see if you can find it
659 Alignment* alignment;
660 if (barcodes.size() > 0) {
661 map<string,int>::iterator it=barcodes.begin();
663 for(it;it!=barcodes.end();it++){
664 if(it->first.length() > maxLength){
665 maxLength = it->first.length();
668 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+bdiffs+1));
670 }else{ alignment = NULL; }
672 //can you find the barcode
678 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
679 string oligo = it->first;
680 // int length = oligo.length();
682 if(rawSequence.length() < maxLength){ //let's just assume that the barcodes are the same length
683 success = bdiffs + 10;
687 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
688 alignment->align(oligo, rawSequence.substr(0,oligo.length()+bdiffs));
689 oligo = alignment->getSeqAAln();
690 string temp = alignment->getSeqBAln();
692 int alnLength = oligo.length();
694 for(int i=oligo.length()-1;i>=0;i--){
695 if(oligo[i] != '-'){ alnLength = i+1; break; }
697 oligo = oligo.substr(0,alnLength);
698 temp = temp.substr(0,alnLength);
700 int numDiff = countDiffs(oligo, temp);
702 if(numDiff < minDiff){
705 minGroup = it->second;
707 for(int i=0;i<alnLength;i++){
713 else if(numDiff == minDiff){
719 if(minDiff > bdiffs) { success = minDiff; } //no good matches
720 else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes
721 else{ //use the best match
723 seq.setUnaligned(rawSequence.substr(minPos));
727 if (alignment != NULL) { delete alignment; }
734 catch(exception& e) {
735 m->errorOut(e, "TrimFlowsCommand", "stripBarcode");
741 //***************************************************************************************************************
743 int TrimFlowsCommand::stripForward(Sequence& seq, int& group){
746 string rawSequence = seq.getUnaligned();
747 int success = pdiffs + 1; //guilty until proven innocent
749 //can you find the primer
750 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
751 string oligo = it->first;
752 if(rawSequence.length() < oligo.length()){ //let's just assume that the primers are the same length
753 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
757 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
759 seq.setUnaligned(rawSequence.substr(oligo.length()));
765 //if you found the barcode or if you don't want to allow for diffs
766 if ((pdiffs == 0) || (success == 0)) { return success; }
768 else { //try aligning and see if you can find it
772 Alignment* alignment;
773 if (primers.size() > 0) {
774 map<string,int>::iterator it=primers.begin();
776 for(it;it!=primers.end();it++){
777 if(it->first.length() > maxLength){
778 maxLength = it->first.length();
781 alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+pdiffs+1));
783 }else{ alignment = NULL; }
785 //can you find the barcode
791 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
792 string oligo = it->first;
793 // int length = oligo.length();
795 if(rawSequence.length() < maxLength){
796 success = pdiffs + 100;
800 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
801 alignment->align(oligo, rawSequence.substr(0,oligo.length()+pdiffs));
802 oligo = alignment->getSeqAAln();
803 string temp = alignment->getSeqBAln();
805 int alnLength = oligo.length();
807 for(int i=oligo.length()-1;i>=0;i--){
808 if(oligo[i] != '-'){ alnLength = i+1; break; }
810 oligo = oligo.substr(0,alnLength);
811 temp = temp.substr(0,alnLength);
813 int numDiff = countDiffs(oligo, temp);
815 if(numDiff < minDiff){
818 minGroup = it->second;
820 for(int i=0;i<alnLength;i++){
826 else if(numDiff == minDiff){
832 if(minDiff > pdiffs) { success = minDiff; } //no good matches
833 else if(minCount > 1) { success = pdiffs + 10; } //can't tell the difference between multiple primers
834 else{ //use the best match
836 seq.setUnaligned(rawSequence.substr(minPos));
840 if (alignment != NULL) { delete alignment; }
847 catch(exception& e) {
848 m->errorOut(e, "TrimFlowsCommand", "stripForward");
853 //***************************************************************************************************************
855 bool TrimFlowsCommand::stripReverse(Sequence& seq){
858 string rawSequence = seq.getUnaligned();
859 bool success = 0; //guilty until proven innocent
861 for(int i=0;i<numRPrimers;i++){
862 string oligo = revPrimer[i];
864 if(rawSequence.length() < oligo.length()){
869 if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
870 seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
879 catch(exception& e) {
880 m->errorOut(e, "TrimFlowsCommand", "stripReverse");
886 //***************************************************************************************************************
888 bool TrimFlowsCommand::compareDNASeq(string oligo, string seq){
891 int length = oligo.length();
893 for(int i=0;i<length;i++){
895 if(oligo[i] != seq[i]){
896 if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C') { success = 0; }
897 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N')) { success = 0; }
898 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G')) { success = 0; }
899 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T')) { success = 0; }
900 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A')) { success = 0; }
901 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G')) { success = 0; }
902 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A')) { success = 0; }
903 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G')) { success = 0; }
904 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G')) { success = 0; }
905 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G')) { success = 0; }
906 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C')) { success = 0; }
907 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G')) { success = 0; }
909 if(success == 0) { break; }
918 catch(exception& e) {
919 m->errorOut(e, "TrimFlowsCommand", "compareDNASeq");
925 //***************************************************************************************************************
927 int TrimFlowsCommand::countDiffs(string oligo, string seq){
930 int length = oligo.length();
933 for(int i=0;i<length;i++){
935 if(oligo[i] != seq[i]){
936 if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C' || oligo[i] == '-' || oligo[i] == '.') { countDiffs++; }
937 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N')) { countDiffs++; }
938 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G')) { countDiffs++; }
939 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T')) { countDiffs++; }
940 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A')) { countDiffs++; }
941 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G')) { countDiffs++; }
942 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A')) { countDiffs++; }
943 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G')) { countDiffs++; }
944 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G')) { countDiffs++; }
945 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G')) { countDiffs++; }
946 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C')) { countDiffs++; }
947 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G')) { countDiffs++; }
954 catch(exception& e) {
955 m->errorOut(e, "TrimFlowsCommand", "countDiffs");
961 /**************************************************************************************************/
963 vector<unsigned long int> TrimFlowsCommand::getFlowFileBreaks() {
967 vector<unsigned long int> filePos;
968 filePos.push_back(0);
971 unsigned long int size;
973 //get num bytes in file
974 pFile = fopen (flowFileName.c_str(),"rb");
975 if (pFile==NULL) perror ("Error opening file");
977 fseek (pFile, 0, SEEK_END);
982 //estimate file breaks
983 unsigned long int chunkSize = 0;
984 chunkSize = size / processors;
986 //file too small to divide by processors
987 if (chunkSize == 0) { processors = 1; filePos.push_back(size); return filePos; }
989 //for each process seekg to closest file break and search for next '>' char. make that the filebreak
990 for (int i = 0; i < processors; i++) {
991 unsigned long int spot = (i+1) * chunkSize;
994 m->openInputFile(flowFileName, in);
997 string dummy = m->getline(in);
999 //there was not another sequence before the end of the file
1000 unsigned long int sanityPos = in.tellg();
1002 // if (sanityPos == -1) { break; }
1003 // else { filePos.push_back(newSpot); }
1004 if (sanityPos == -1) { break; }
1005 else { filePos.push_back(sanityPos); }
1011 filePos.push_back(size);
1013 //sanity check filePos
1014 for (int i = 0; i < (filePos.size()-1); i++) {
1015 if (filePos[(i+1)] <= filePos[i]) { filePos.erase(filePos.begin()+(i+1)); i--; }
1019 m->openInputFile(flowFileName, in);
1022 unsigned long int spot = in.tellg();
1026 processors = (filePos.size() - 1);
1030 catch(exception& e) {
1031 m->errorOut(e, "TrimSeqsCommand", "getFlowFileBreaks");
1036 /**************************************************************************************************/
1038 int TrimFlowsCommand::createProcessesCreateTrim(string flowFileName, string trimFlowFileName, string scrapFlowFileName, string fastaFileName, vector<vector<string> > barcodePrimerComboFileNames){
1041 #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
1043 int exitCommand = 1;
1046 //loop through and create all the processes you want
1047 while (process != processors) {
1051 processIDS.push_back(pid); //create map from line number to pid so you can append files in correct order later
1053 }else if (pid == 0){
1055 vector<vector<string> > tempBarcodePrimerComboFileNames = barcodePrimerComboFileNames;
1057 for(int i=0;i<tempBarcodePrimerComboFileNames.size();i++){
1058 for(int j=0;j<tempBarcodePrimerComboFileNames[0].size();j++){
1059 tempBarcodePrimerComboFileNames[i][j] += toString(getpid()) + ".temp";
1061 m->openOutputFile(tempBarcodePrimerComboFileNames[i][j], temp);
1067 driverCreateTrim(flowFileName,
1068 (trimFlowFileName + toString(getpid()) + ".temp"),
1069 (scrapFlowFileName + toString(getpid()) + ".temp"),
1070 (fastaFileName + toString(getpid()) + ".temp"),
1071 tempBarcodePrimerComboFileNames, lines[process]);
1075 m->mothurOut("[ERROR]: unable to spawn the necessary processes."); m->mothurOutEndLine();
1076 for (int i = 0; i < processIDS.size(); i++) { kill (processIDS[i], SIGINT); }
1083 m->openOutputFile(trimFlowFileName, temp);
1086 m->openOutputFile(scrapFlowFileName, temp);
1090 m->openOutputFile(fastaFileName, temp);
1094 driverCreateTrim(flowFileName, trimFlowFileName, scrapFlowFileName, fastaFileName, barcodePrimerComboFileNames, lines[0]);
1096 //force parent to wait until all the processes are done
1097 for (int i=0;i<processIDS.size();i++) {
1098 int temp = processIDS[i];
1103 m->mothurOutEndLine();
1104 for(int i=0;i<processIDS.size();i++){
1106 m->mothurOut("Appending files from process " + toString(processIDS[i])); m->mothurOutEndLine();
1108 m->appendFiles((trimFlowFileName + toString(processIDS[i]) + ".temp"), trimFlowFileName);
1109 remove((trimFlowFileName + toString(processIDS[i]) + ".temp").c_str());
1110 // m->mothurOut("\tDone with trim.flow file"); m->mothurOutEndLine();
1112 m->appendFiles((scrapFlowFileName + toString(processIDS[i]) + ".temp"), scrapFlowFileName);
1113 remove((scrapFlowFileName + toString(processIDS[i]) + ".temp").c_str());
1114 // m->mothurOut("\tDone with scrap.flow file"); m->mothurOutEndLine();
1117 m->appendFiles((fastaFileName + toString(processIDS[i]) + ".temp"), fastaFileName);
1118 remove((fastaFileName + toString(processIDS[i]) + ".temp").c_str());
1119 // m->mothurOut("\tDone with flow.fasta file"); m->mothurOutEndLine();
1122 for (int j = 0; j < barcodePrimerComboFileNames.size(); j++) {
1123 for (int k = 0; k < barcodePrimerComboFileNames[0].size(); k++) {
1124 m->appendFiles((barcodePrimerComboFileNames[j][k] + toString(processIDS[i]) + ".temp"), barcodePrimerComboFileNames[j][k]);
1125 remove((barcodePrimerComboFileNames[j][k] + toString(processIDS[i]) + ".temp").c_str());
1134 catch(exception& e) {
1135 m->errorOut(e, "TrimFlowsCommand", "createProcessesCreateTrim");
1140 //***************************************************************************************************************