5 * Created by westcott on 9/1/11.
6 * Copyright 2011 Schloss Lab. All rights reserved.
10 #include "trimoligos.h"
11 #include "alignment.hpp"
12 #include "needlemanoverlap.hpp"
15 /********************************************************************/
16 //strip, pdiffs, bdiffs, primers, barcodes, revPrimers
17 TrimOligos::TrimOligos(int p, int b, int l, int s, map<string, int> pr, map<string, int> br, vector<string> r, vector<string> lk, vector<string> sp){
19 m = MothurOut::getInstance();
32 maxFBarcodeLength = 0;
33 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
34 if(it->first.length() > maxFBarcodeLength){
35 maxFBarcodeLength = it->first.length();
39 map<string,int>::iterator it;
40 for(it=primers.begin();it!=primers.end();it++){
41 if(it->first.length() > maxFPrimerLength){
42 maxFPrimerLength = it->first.length();
47 for(int i = 0; i < linker.size(); i++){
48 if(linker[i].length() > maxLinkerLength){
49 maxLinkerLength = linker[i].length();
54 for(int i = 0; i < spacer.size(); i++){
55 if(spacer[i].length() > maxSpacerLength){
56 maxSpacerLength = spacer[i].length();
61 m->errorOut(e, "TrimOligos", "TrimOligos");
65 /********************************************************************/
66 //strip, pdiffs, bdiffs, primers, barcodes, revPrimers
67 TrimOligos::TrimOligos(int p, int b, int l, int s, map<int, oligosPair> pr, map<int, oligosPair> br){
69 m = MothurOut::getInstance();
77 maxFBarcodeLength = 0;
78 for(map<int,oligosPair>::iterator it=br.begin();it!=br.end();it++){
79 string forward = it->second.forward;
80 map<string, vector<int> >::iterator itForward = ifbarcodes.find(forward);
82 if(forward.length() > maxFBarcodeLength){ maxFBarcodeLength = forward.length(); }
84 if (itForward == ifbarcodes.end()) {
85 vector<int> temp; temp.push_back(it->first);
86 ifbarcodes[forward] = temp;
87 }else { itForward->second.push_back(it->first); }
91 for(map<int,oligosPair>::iterator it=pr.begin();it!=pr.end();it++){
92 string forward = it->second.forward;
93 map<string, vector<int> >::iterator itForward = ifprimers.find(forward);
95 if(forward.length() > maxFPrimerLength){ maxFPrimerLength = forward.length(); }
97 if (itForward == ifprimers.end()) {
98 vector<int> temp; temp.push_back(it->first);
99 ifprimers[forward] = temp;
100 }else { itForward->second.push_back(it->first); }
103 maxRBarcodeLength = 0;
104 for(map<int,oligosPair>::iterator it=br.begin();it!=br.end();it++){
105 string forward = it->second.reverse;
106 map<string, vector<int> >::iterator itForward = irbarcodes.find(forward);
108 if(forward.length() > maxRBarcodeLength){ maxRBarcodeLength = forward.length(); }
110 if (itForward == irbarcodes.end()) {
111 vector<int> temp; temp.push_back(it->first);
112 irbarcodes[forward] = temp;
113 }else { itForward->second.push_back(it->first); }
116 maxRPrimerLength = 0;
117 for(map<int,oligosPair>::iterator it=pr.begin();it!=pr.end();it++){
118 string forward = it->second.reverse;
119 map<string, vector<int> >::iterator itForward = irprimers.find(forward);
121 if(forward.length() > maxRPrimerLength){ maxRPrimerLength = forward.length(); }
123 if (itForward == irprimers.end()) {
124 vector<int> temp; temp.push_back(it->first);
125 irprimers[forward] = temp;
126 }else { itForward->second.push_back(it->first); }
132 catch(exception& e) {
133 m->errorOut(e, "TrimOligos", "TrimOligos");
137 /********************************************************************/
138 //strip, pdiffs, bdiffs, primers, barcodes, revPrimers
139 TrimOligos::TrimOligos(int p, int b, map<string, int> pr, map<string, int> br, vector<string> r){
141 m = MothurOut::getInstance();
151 maxFBarcodeLength = 0;
152 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
153 string oligo = it->first;
154 if(oligo.length() > maxFBarcodeLength){
155 maxFBarcodeLength = oligo.length();
158 maxRBarcodeLength = maxFBarcodeLength;
160 maxFPrimerLength = 0;
161 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
162 string oligo = it->first;
163 if(oligo.length() > maxFPrimerLength){
164 maxFPrimerLength = oligo.length();
167 maxRPrimerLength = maxFPrimerLength;
169 catch(exception& e) {
170 m->errorOut(e, "TrimOligos", "TrimOligos");
174 /********************************************************************/
175 TrimOligos::~TrimOligos() {}
176 //********************************************************************/
177 bool TrimOligos::findForward(Sequence& seq, int& primerStart, int& primerEnd){
180 string rawSequence = seq.getUnaligned();
182 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
183 string oligo = it->first;
185 if(rawSequence.length() < oligo.length()) { break; }
188 int olength = oligo.length();
189 for (int j = 0; j < rawSequence.length()-olength; j++){
190 if (m->control_pressed) { primerStart = 0; primerEnd = 0; return false; }
191 string rawChunk = rawSequence.substr(j, olength);
192 if(compareDNASeq(oligo, rawChunk)) {
194 primerEnd = primerStart + olength;
201 primerStart = 0; primerEnd = 0;
202 //if you don't want to allow for diffs
203 if (pdiffs == 0) { return false; }
204 else { //try aligning and see if you can find it
206 //can you find the barcode
210 Alignment* alignment;
211 if (primers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFPrimerLength+pdiffs+1)); }
212 else{ alignment = NULL; }
214 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
216 string prim = it->first;
218 int olength = prim.length();
219 if (rawSequence.length() < olength) {} //ignore primers too long for this seq
221 for (int j = 0; j < rawSequence.length()-olength; j++){
223 string oligo = it->first;
225 if (m->control_pressed) { primerStart = 0; primerEnd = 0; return false; }
227 string rawChunk = rawSequence.substr(j, olength+pdiffs);
229 //use needleman to align first primer.length()+numdiffs of sequence to each barcode
230 alignment->alignPrimer(oligo, rawChunk);
231 oligo = alignment->getSeqAAln();
232 string temp = alignment->getSeqBAln();
234 int alnLength = oligo.length();
236 for(int i=oligo.length()-1;i>=0;i--){
237 if(oligo[i] != '-'){ alnLength = i+1; break; }
239 oligo = oligo.substr(0,alnLength);
240 temp = temp.substr(0,alnLength);
242 int numDiff = countDiffs(oligo, temp);
244 if(numDiff < minDiff){
248 primerEnd = primerStart + alnLength;
249 }else if(numDiff == minDiff){ minCount++; }
254 if (alignment != NULL) { delete alignment; }
256 if(minDiff > pdiffs) { primerStart = 0; primerEnd = 0; return false; } //no good matches
257 else if(minCount > 1) { primerStart = 0; primerEnd = 0; return false; } //can't tell the difference between multiple primers
261 primerStart = 0; primerEnd = 0;
265 catch(exception& e) {
266 m->errorOut(e, "TrimOligos", "stripForward");
270 //******************************************************************/
271 bool TrimOligos::findReverse(Sequence& seq, int& primerStart, int& primerEnd){
274 string rawSequence = seq.getUnaligned();
276 for(int i=0;i<revPrimer.size();i++){
277 string oligo = revPrimer[i];
278 if(rawSequence.length() < oligo.length()) { break; }
281 int olength = oligo.length();
282 for (int j = rawSequence.length()-olength; j >= 0; j--){
283 if (m->control_pressed) { primerStart = 0; primerEnd = 0; return false; }
284 string rawChunk = rawSequence.substr(j, olength);
286 if(compareDNASeq(oligo, rawChunk)) {
288 primerEnd = primerStart + olength;
295 primerStart = 0; primerEnd = 0;
298 catch(exception& e) {
299 m->errorOut(e, "PcrSeqsCommand", "findReverse");
303 //*******************************************************************/
305 int TrimOligos::stripBarcode(Sequence& seq, QualityScores& qual, int& group){
308 if (paired) { int success = stripPairedBarcode(seq, qual, group); return success; }
310 string rawSequence = seq.getUnaligned();
311 int success = bdiffs + 1; //guilty until proven innocent
313 //can you find the barcode
314 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
315 string oligo = it->first;
316 if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length
317 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
321 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
323 seq.setUnaligned(rawSequence.substr(oligo.length()));
325 if(qual.getName() != ""){
326 qual.trimQScores(oligo.length(), -1);
334 //if you found the barcode or if you don't want to allow for diffs
335 if ((bdiffs == 0) || (success == 0)) { return success; }
337 else { //try aligning and see if you can find it
338 Alignment* alignment;
339 if (barcodes.size() > 0) {alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFBarcodeLength+bdiffs+1)); }
340 else{ alignment = NULL; }
342 //can you find the barcode
348 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
349 string oligo = it->first;
350 // int length = oligo.length();
352 if(rawSequence.length() < maxFBarcodeLength){ //let's just assume that the barcodes are the same length
353 success = bdiffs + 10;
357 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
358 alignment->alignPrimer(oligo, rawSequence.substr(0,oligo.length()+bdiffs));
359 oligo = alignment->getSeqAAln();
360 string temp = alignment->getSeqBAln();
362 int alnLength = oligo.length();
364 for(int i=oligo.length()-1;i>=0;i--){
365 if(oligo[i] != '-'){ alnLength = i+1; break; }
367 oligo = oligo.substr(0,alnLength);
368 temp = temp.substr(0,alnLength);
369 int numDiff = countDiffs(oligo, temp);
371 if (m->debug) { m->mothurOut("[DEBUG]: " + seq.getName() + " aligned fragment =" + temp + ", barcode =" + oligo + ", numDiffs = " + toString(numDiff) + "\n"); }
373 if(numDiff < minDiff){
376 minGroup = it->second;
378 for(int i=0;i<alnLength;i++){
384 else if(numDiff == minDiff){
390 if(minDiff > bdiffs) { success = minDiff; } //no good matches
391 else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes
392 else{ //use the best match
394 seq.setUnaligned(rawSequence.substr(minPos));
396 if(qual.getName() != ""){
397 qual.trimQScores(minPos, -1);
402 if (alignment != NULL) { delete alignment; }
409 catch(exception& e) {
410 m->errorOut(e, "TrimOligos", "stripBarcode");
414 //*******************************************************************/
415 int TrimOligos::stripBarcode(Sequence& forwardSeq, Sequence& reverseSeq, int& group){
417 //look for forward barcode
418 string rawFSequence = forwardSeq.getUnaligned();
419 string rawRSequence = reverseSeq.getUnaligned();
420 int success = bdiffs + 1; //guilty until proven innocent
422 //can you find the forward barcode
423 for(map<int,oligosPair>::iterator it=ipbarcodes.begin();it!=ipbarcodes.end();it++){
424 string foligo = it->second.forward;
425 string roligo = it->second.reverse;
427 if(rawFSequence.length() < foligo.length()){ //let's just assume that the barcodes are the same length
428 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
431 if(rawRSequence.length() < roligo.length()){ //let's just assume that the barcodes are the same length
432 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
436 if((compareDNASeq(foligo, rawFSequence.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawRSequence.substr(0,roligo.length())))) {
438 forwardSeq.setUnaligned(rawFSequence.substr(foligo.length()));
439 reverseSeq.setUnaligned(rawRSequence.substr(roligo.length()));
445 //if you found the barcode or if you don't want to allow for diffs
446 if ((bdiffs == 0) || (success == 0)) { return success; }
447 else { //try aligning and see if you can find it
448 Alignment* alignment;
449 if (ifbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFBarcodeLength+bdiffs+1)); }
450 else{ alignment = NULL; }
452 //can you find the barcode
455 vector< vector<int> > minFGroup;
458 //the pair can have duplicates, but we only want to search for the unique forward and reverses, example
467 only want to look for forward = Sarah, John, Anna, Pat, Gail
468 reverse = Westcott, Schloss, Brown, Moore
469 but if best match forward = 4, and reverse = 1, we want to count as a valid match because forward 1 and forward 4 are the same. so both barcodes map to same group.
471 //cout << endl << forwardSeq.getName() << endl;
472 for(map<string, vector<int> >::iterator it=ifbarcodes.begin();it!=ifbarcodes.end();it++){
473 string oligo = it->first;
475 if(rawFSequence.length() < maxFBarcodeLength){ //let's just assume that the barcodes are the same length
476 success = bdiffs + 10;
479 //cout << "before = " << oligo << '\t' << rawFSequence.substr(0,oligo.length()+bdiffs) << endl;
480 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
481 alignment->alignPrimer(oligo, rawFSequence.substr(0,oligo.length()+bdiffs));
482 oligo = alignment->getSeqAAln();
483 string temp = alignment->getSeqBAln();
485 int alnLength = oligo.length();
487 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
488 oligo = oligo.substr(0,alnLength);
489 temp = temp.substr(0,alnLength);
490 int numDiff = countDiffs(oligo, temp);
492 if (m->debug) { m->mothurOut("[DEBUG]: forward " + forwardSeq.getName() + " aligned fragment=" + temp + ", barcode=" + oligo + ", numDiffs=" + toString(numDiff) + ".\n"); }
494 if (alnLength == 0) { numDiff = bdiffs + 100; }
495 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
497 if(numDiff < minDiff){
501 minFGroup.push_back(it->second);
504 for(int i=0;i<alnLength;i++){
509 minFPos.push_back(tempminFPos);
510 }else if(numDiff == minDiff){
511 minFGroup.push_back(it->second);
513 for(int i=0;i<alnLength;i++){
518 minFPos.push_back(tempminFPos);
522 //cout << minDiff << '\t' << minCount << '\t' << endl;
523 if(minDiff > bdiffs) { success = minDiff; } //no good matches
525 //check for reverse match
526 if (alignment != NULL) { delete alignment; }
528 if (irbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxRBarcodeLength+bdiffs+1)); }
529 else{ alignment = NULL; }
531 //can you find the barcode
534 vector< vector<int> > minRGroup;
537 for(map<string, vector<int> >::iterator it=irbarcodes.begin();it!=irbarcodes.end();it++){
538 string oligo = it->first;
539 //cout << "before = " << oligo << '\t' << rawRSequence.substr(0,oligo.length()+bdiffs) << endl;
540 if(rawRSequence.length() < maxRBarcodeLength){ //let's just assume that the barcodes are the same length
541 success = bdiffs + 10;
545 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
546 alignment->alignPrimer(oligo, rawRSequence.substr(0,oligo.length()+bdiffs));
547 oligo = alignment->getSeqAAln();
548 string temp = alignment->getSeqBAln();
550 int alnLength = oligo.length();
551 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
552 oligo = oligo.substr(0,alnLength);
553 temp = temp.substr(0,alnLength);
554 int numDiff = countDiffs(oligo, temp);
556 if (m->debug) { m->mothurOut("[DEBUG]: reverse " + forwardSeq.getName() + " aligned fragment=" + temp + ", barcode=" + oligo + ", numDiffs=" + toString(numDiff) + ".\n"); }
558 if (alnLength == 0) { numDiff = bdiffs + 100; }
560 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
561 if(numDiff < minDiff){
565 minRGroup.push_back(it->second);
568 for(int i=0;i<alnLength;i++){
573 minRPos.push_back(tempminRPos);
574 }else if(numDiff == minDiff){
576 for(int i=0;i<alnLength;i++){
581 minRPos.push_back(tempminRPos);
582 minRGroup.push_back(it->second);
587 if(minDiff > bdiffs) { success = minDiff; } //no good matches
589 bool foundMatch = false;
591 for (int i = 0; i < minFGroup.size(); i++) {
592 for (int j = 0; j < minFGroup[i].size(); j++) {
593 for (int k = 0; k < minRGroup.size(); k++) {
594 if (m->inUsersGroups(minFGroup[i][j], minRGroup[k])) { matches.push_back(minFGroup[i][j]); k+= minRGroup.size(); }
601 if (matches.size() == 1) {
608 //we have an acceptable match for the forward and reverse, but do they match?
610 forwardSeq.setUnaligned(rawFSequence.substr(fStart));
611 reverseSeq.setUnaligned(rawRSequence.substr(rStart));
613 }else { success = bdiffs + 100; }
617 if (alignment != NULL) { delete alignment; }
623 catch(exception& e) {
624 m->errorOut(e, "TrimOligos", "stripIBarcode");
629 //*******************************************************************/
630 int TrimOligos::stripBarcode(Sequence& forwardSeq, Sequence& reverseSeq, QualityScores& forwardQual, QualityScores& reverseQual, int& group){
632 //look for forward barcode
633 string rawFSequence = forwardSeq.getUnaligned();
634 string rawRSequence = reverseSeq.getUnaligned();
635 int success = bdiffs + 1; //guilty until proven innocent
637 //can you find the forward barcode
638 for(map<int,oligosPair>::iterator it=ipbarcodes.begin();it!=ipbarcodes.end();it++){
639 string foligo = it->second.forward;
640 string roligo = it->second.reverse;
642 if(rawFSequence.length() < foligo.length()){ //let's just assume that the barcodes are the same length
643 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
646 if(rawRSequence.length() < roligo.length()){ //let's just assume that the barcodes are the same length
647 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
651 if((compareDNASeq(foligo, rawFSequence.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawRSequence.substr(0,roligo.length())))) {
653 forwardSeq.setUnaligned(rawFSequence.substr(foligo.length()));
654 reverseSeq.setUnaligned(rawRSequence.substr(roligo.length()));
655 forwardQual.trimQScores(foligo.length(), -1);
656 reverseQual.trimQScores(roligo.length(), -1);
662 //if you found the barcode or if you don't want to allow for diffs
663 if ((bdiffs == 0) || (success == 0)) { return success; }
664 else { //try aligning and see if you can find it
665 Alignment* alignment;
666 if (ifbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFBarcodeLength+bdiffs+1)); }
667 else{ alignment = NULL; }
669 //can you find the barcode
672 vector< vector<int> > minFGroup;
675 //the pair can have duplicates, but we only want to search for the unique forward and reverses, example
684 only want to look for forward = Sarah, John, Anna, Pat, Gail
685 reverse = Westcott, Schloss, Brown, Moore
686 but if best match forward = 4, and reverse = 1, we want to count as a valid match because forward 1 and forward 4 are the same. so both barcodes map to same group.
688 //cout << endl << forwardSeq.getName() << endl;
689 for(map<string, vector<int> >::iterator it=ifbarcodes.begin();it!=ifbarcodes.end();it++){
690 string oligo = it->first;
692 if(rawFSequence.length() < maxFBarcodeLength){ //let's just assume that the barcodes are the same length
693 success = bdiffs + 10;
696 //cout << "before = " << oligo << '\t' << rawFSequence.substr(0,oligo.length()+bdiffs) << endl;
697 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
698 alignment->alignPrimer(oligo, rawFSequence.substr(0,oligo.length()+bdiffs));
699 oligo = alignment->getSeqAAln();
700 string temp = alignment->getSeqBAln();
702 int alnLength = oligo.length();
704 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
705 oligo = oligo.substr(0,alnLength);
706 temp = temp.substr(0,alnLength);
707 int numDiff = countDiffs(oligo, temp);
709 if (alnLength == 0) { numDiff = bdiffs + 100; }
710 if (m->debug) { m->mothurOut("[DEBUG]: forward " + forwardSeq.getName() + " aligned fragment=" + temp + ", barcode=" + oligo + ", numDiffs=" + toString(numDiff) + ".\n"); }
712 if(numDiff < minDiff){
716 minFGroup.push_back(it->second);
719 for(int i=0;i<alnLength;i++){
724 minFPos.push_back(tempminFPos);
725 }else if(numDiff == minDiff){
726 minFGroup.push_back(it->second);
728 for(int i=0;i<alnLength;i++){
733 minFPos.push_back(tempminFPos);
737 //cout << minDiff << '\t' << minCount << '\t' << endl;
738 if(minDiff > bdiffs) { success = minDiff; } //no good matches
740 //check for reverse match
741 if (alignment != NULL) { delete alignment; }
743 if (irbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxRBarcodeLength+bdiffs+1)); }
744 else{ alignment = NULL; }
746 //can you find the barcode
749 vector< vector<int> > minRGroup;
752 for(map<string, vector<int> >::iterator it=irbarcodes.begin();it!=irbarcodes.end();it++){
753 string oligo = it->first;
754 //cout << "before = " << oligo << '\t' << rawRSequence.substr(0,oligo.length()+bdiffs) << endl;
755 if(rawRSequence.length() < maxRBarcodeLength){ //let's just assume that the barcodes are the same length
756 success = bdiffs + 10;
760 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
761 alignment->alignPrimer(oligo, rawRSequence.substr(0,oligo.length()+bdiffs));
762 oligo = alignment->getSeqAAln();
763 string temp = alignment->getSeqBAln();
765 int alnLength = oligo.length();
766 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
767 oligo = oligo.substr(0,alnLength);
768 temp = temp.substr(0,alnLength);
769 int numDiff = countDiffs(oligo, temp);
770 if (alnLength == 0) { numDiff = bdiffs + 100; }
772 if (m->debug) { m->mothurOut("[DEBUG]: reverse " + reverseSeq.getName() + " aligned fragment=" + temp + ", barcode=" + oligo + ", numDiffs=" + toString(numDiff) + ".\n"); }
774 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
775 if(numDiff < minDiff){
779 minRGroup.push_back(it->second);
782 for(int i=0;i<alnLength;i++){
787 minRPos.push_back(tempminRPos);
788 }else if(numDiff == minDiff){
790 for(int i=0;i<alnLength;i++){
795 minRPos.push_back(tempminRPos);
796 minRGroup.push_back(it->second);
801 if(minDiff > bdiffs) { success = minDiff; } //no good matches
804 string fMatches = "";
805 for (int i = 0; i < minFGroup.size(); i++) { for (int j = 0; j < minFGroup[i].size(); j++) { fMatches += toString(minFGroup[i][j]) + " "; } }
806 m->mothurOut("[DEBUG]: forward matches = " + fMatches + "\n");
807 string rMatches = "";
808 for (int i = 0; i < minRGroup.size(); i++) { for (int j = 0; j < minRGroup[i].size(); j++) { rMatches += toString(minRGroup[i][j]) + " "; } }
809 m->mothurOut("[DEBUG]: reverse matches = " + rMatches + "\n");
811 bool foundMatch = false;
813 for (int i = 0; i < minFGroup.size(); i++) {
814 for (int j = 0; j < minFGroup[i].size(); j++) {
815 for (int k = 0; k < minRGroup.size(); k++) {
816 if (m->inUsersGroups(minFGroup[i][j], minRGroup[k])) { matches.push_back(minFGroup[i][j]); k+= minRGroup.size(); }
823 if (matches.size() == 1) {
828 if (m->debug) { m->mothurOut("[DEBUG]: found match = " + toString(matches[0]) + ".\n"); }
830 if (m->debug) { m->mothurOut("[DEBUG]: failed too many matches.\n"); }
833 //we have an acceptable match for the forward and reverse, but do they match?
835 forwardSeq.setUnaligned(rawFSequence.substr(fStart));
836 reverseSeq.setUnaligned(rawRSequence.substr(rStart));
837 forwardQual.trimQScores(fStart, -1);
838 reverseQual.trimQScores(rStart, -1);
840 }else { success = bdiffs + 100; }
844 if (alignment != NULL) { delete alignment; }
850 catch(exception& e) {
851 m->errorOut(e, "TrimOligos", "stripIBarcode");
856 //*******************************************************************/
857 int TrimOligos::stripPairedBarcode(Sequence& seq, QualityScores& qual, int& group){
859 //look for forward barcode
860 string rawSeq = seq.getUnaligned();
861 int success = bdiffs + 1; //guilty until proven innocent
862 //cout << seq.getName() << endl;
863 //can you find the forward barcode
864 for(map<int,oligosPair>::iterator it=ipbarcodes.begin();it!=ipbarcodes.end();it++){
865 string foligo = it->second.forward;
866 string roligo = it->second.reverse;
867 //cout << it->first << '\t' << foligo << '\t' << roligo << endl;
868 //cout << it->first << '\t' << rawSeq.substr(0,foligo.length()) << '\t' << rawSeq.substr(rawSeq.length()-roligo.length(),roligo.length()) << endl;
869 if(rawSeq.length() < foligo.length()){ //let's just assume that the barcodes are the same length
870 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
873 if(rawSeq.length() < roligo.length()){ //let's just assume that the barcodes are the same length
874 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
878 if (rawSeq.length() < (foligo.length() + roligo.length())) { success = pdiffs + 10; break; }
880 if((compareDNASeq(foligo, rawSeq.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawSeq.substr(rawSeq.length()-roligo.length(),roligo.length())))) {
882 string trimmedSeq = rawSeq.substr(foligo.length()); //trim forward barcode
883 seq.setUnaligned(trimmedSeq.substr(0,(trimmedSeq.length()-roligo.length()))); //trim reverse barcode
884 if(qual.getName() != ""){
885 qual.trimQScores(-1, rawSeq.length()-roligo.length());
886 qual.trimQScores(foligo.length(), -1);
892 //cout << "success=" << success << endl;
893 //if you found the barcode or if you don't want to allow for diffs
894 if ((bdiffs == 0) || (success == 0)) { return success; }
895 else { //try aligning and see if you can find it
896 Alignment* alignment;
897 if (ifbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFBarcodeLength+bdiffs+1)); }
898 else{ alignment = NULL; }
900 //can you find the barcode
903 vector< vector<int> > minFGroup;
906 //the pair can have duplicates, but we only want to search for the unique forward and reverses, example
915 only want to look for forward = Sarah, John, Anna, Pat, Gail
916 reverse = Westcott, Schloss, Brown, Moore
917 but if best match forward = 4, and reverse = 1, we want to count as a valid match because forward 1 and forward 4 are the same. so both barcodes map to same group.
919 //cout << endl << forwardSeq.getName() << endl;
920 for(map<string, vector<int> >::iterator it=ifbarcodes.begin();it!=ifbarcodes.end();it++){
921 string oligo = it->first;
923 if(rawSeq.length() < maxFBarcodeLength){ //let's just assume that the barcodes are the same length
924 success = bdiffs + 10;
927 //cout << "before = " << oligo << '\t' << rawSeq.substr(0,oligo.length()+bdiffs) << endl;
928 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
929 alignment->alignPrimer(oligo, rawSeq.substr(0,oligo.length()+bdiffs));
930 oligo = alignment->getSeqAAln();
931 string temp = alignment->getSeqBAln();
933 int alnLength = oligo.length();
935 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
936 oligo = oligo.substr(0,alnLength);
937 temp = temp.substr(0,alnLength);
938 int numDiff = countDiffs(oligo, temp);
940 if (alnLength == 0) { numDiff = bdiffs + 100; }
941 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
943 if(numDiff < minDiff){
947 minFGroup.push_back(it->second);
950 for(int i=0;i<alnLength;i++){
955 minFPos.push_back(tempminFPos);
956 }else if(numDiff == minDiff){
957 minFGroup.push_back(it->second);
959 for(int i=0;i<alnLength;i++){
964 minFPos.push_back(tempminFPos);
968 //cout << minDiff << '\t' << minCount << '\t' << endl;
969 if(minDiff > bdiffs) { success = minDiff; } //no good matches
971 //check for reverse match
972 if (alignment != NULL) { delete alignment; }
974 if (irbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxRBarcodeLength+bdiffs+1)); }
975 else{ alignment = NULL; }
977 //can you find the barcode
980 vector< vector<int> > minRGroup;
983 string rawRSequence = reverseOligo(seq.getUnaligned());
984 //cout << irbarcodes.size() << '\t' << maxRBarcodeLength << endl;
985 for(map<string, vector<int> >::iterator it=irbarcodes.begin();it!=irbarcodes.end();it++){
986 string oligo = reverseOligo(it->first);
987 //cout << "r before = " << reverseOligo(oligo) << '\t' << reverseOligo(rawRSequence.substr(0,oligo.length()+bdiffs)) << endl;
988 if(rawRSequence.length() < maxRBarcodeLength){ //let's just assume that the barcodes are the same length
989 success = bdiffs + 10;
993 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
994 alignment->alignPrimer(oligo, rawRSequence.substr(0,oligo.length()+bdiffs));
995 oligo = alignment->getSeqAAln();
996 string temp = alignment->getSeqBAln();
998 int alnLength = oligo.length();
999 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1000 oligo = oligo.substr(0,alnLength);
1001 temp = temp.substr(0,alnLength);
1002 int numDiff = countDiffs(oligo, temp);
1003 if (alnLength == 0) { numDiff = bdiffs + 100; }
1005 //cout << "r after = " << reverseOligo(oligo) << '\t' << reverseOligo(temp) << '\t' << numDiff << endl;
1006 if(numDiff < minDiff){
1010 minRGroup.push_back(it->second);
1011 int tempminRPos = 0;
1013 for(int i=0;i<alnLength;i++){
1018 minRPos.push_back(tempminRPos);
1019 }else if(numDiff == minDiff){
1020 int tempminRPos = 0;
1021 for(int i=0;i<alnLength;i++){
1026 minRPos.push_back(tempminRPos);
1027 minRGroup.push_back(it->second);
1032 if(minDiff > bdiffs) { success = minDiff; } //no good matches
1034 bool foundMatch = false;
1035 vector<int> matches;
1036 for (int i = 0; i < minFGroup.size(); i++) {
1037 for (int j = 0; j < minFGroup[i].size(); j++) {
1038 for (int k = 0; k < minRGroup.size(); k++) {
1039 if (m->inUsersGroups(minFGroup[i][j], minRGroup[k])) { matches.push_back(minFGroup[i][j]); k+= minRGroup.size(); }
1046 if (matches.size() == 1) {
1049 fStart = minFPos[0];
1050 rStart = rawSeq.length() - minRPos[0];
1051 if (fStart > rStart) { foundMatch = false; } //only barcodes not a good sequence
1054 //we have an acceptable match for the forward and reverse, but do they match?
1056 string trimmedSeq = rawSeq.substr(0, rStart); //trim reverse barcode
1057 seq.setUnaligned(trimmedSeq.substr(fStart)); //trim forward barcode
1058 if(qual.getName() != ""){
1059 qual.trimQScores(-1, rStart);
1060 qual.trimQScores(fStart, -1);
1063 //cout << "barcode = " << ipbarcodes[group].forward << '\t' << ipbarcodes[group].reverse << endl;
1064 }else { success = bdiffs + 100; }
1068 if (alignment != NULL) { delete alignment; }
1074 catch(exception& e) {
1075 m->errorOut(e, "TrimOligos", "stripPairedBarcode");
1080 //*******************************************************************/
1081 int TrimOligos::stripPairedBarcode(Sequence& seq, int& group){
1083 //look for forward barcode
1084 string rawSeq = seq.getUnaligned();
1085 int success = bdiffs + 1; //guilty until proven innocent
1086 //cout << seq.getName() << endl;
1087 //can you find the forward barcode
1088 for(map<int,oligosPair>::iterator it=ipbarcodes.begin();it!=ipbarcodes.end();it++){
1089 string foligo = it->second.forward;
1090 string roligo = it->second.reverse;
1091 //cout << it->first << '\t' << foligo << '\t' << roligo << endl;
1092 //cout << it->first << '\t' << rawSeq.substr(0,foligo.length()) << '\t' << rawSeq.substr(rawSeq.length()-roligo.length(),roligo.length()) << endl;
1093 if(rawSeq.length() < foligo.length()){ //let's just assume that the barcodes are the same length
1094 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
1097 if(rawSeq.length() < roligo.length()){ //let's just assume that the barcodes are the same length
1098 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
1102 if (rawSeq.length() < (foligo.length() + roligo.length())) { success = pdiffs + 10; break; }
1104 if((compareDNASeq(foligo, rawSeq.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawSeq.substr(rawSeq.length()-roligo.length(),roligo.length())))) {
1106 string trimmedSeq = rawSeq.substr(foligo.length()); //trim forward barcode
1107 seq.setUnaligned(trimmedSeq.substr(0,(trimmedSeq.length()-roligo.length()))); //trim reverse barcode
1112 //cout << "success=" << success << endl;
1113 //if you found the barcode or if you don't want to allow for diffs
1114 if ((bdiffs == 0) || (success == 0)) { return success; }
1115 else { //try aligning and see if you can find it
1116 Alignment* alignment;
1117 if (ifbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFBarcodeLength+bdiffs+1)); }
1118 else{ alignment = NULL; }
1120 //can you find the barcode
1123 vector< vector<int> > minFGroup;
1124 vector<int> minFPos;
1126 //the pair can have duplicates, but we only want to search for the unique forward and reverses, example
1135 only want to look for forward = Sarah, John, Anna, Pat, Gail
1136 reverse = Westcott, Schloss, Brown, Moore
1137 but if best match forward = 4, and reverse = 1, we want to count as a valid match because forward 1 and forward 4 are the same. so both barcodes map to same group.
1139 //cout << endl << forwardSeq.getName() << endl;
1140 for(map<string, vector<int> >::iterator it=ifbarcodes.begin();it!=ifbarcodes.end();it++){
1141 string oligo = it->first;
1143 if(rawSeq.length() < maxFBarcodeLength){ //let's just assume that the barcodes are the same length
1144 success = bdiffs + 10;
1147 //cout << "before = " << oligo << '\t' << rawSeq.substr(0,oligo.length()+bdiffs) << endl;
1148 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1149 alignment->alignPrimer(oligo, rawSeq.substr(0,oligo.length()+bdiffs));
1150 oligo = alignment->getSeqAAln();
1151 string temp = alignment->getSeqBAln();
1153 int alnLength = oligo.length();
1155 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1156 oligo = oligo.substr(0,alnLength);
1157 temp = temp.substr(0,alnLength);
1158 int numDiff = countDiffs(oligo, temp);
1160 if (alnLength == 0) { numDiff = bdiffs + 100; }
1161 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
1163 if(numDiff < minDiff){
1167 minFGroup.push_back(it->second);
1168 int tempminFPos = 0;
1170 for(int i=0;i<alnLength;i++){
1175 minFPos.push_back(tempminFPos);
1176 }else if(numDiff == minDiff){
1177 minFGroup.push_back(it->second);
1178 int tempminFPos = 0;
1179 for(int i=0;i<alnLength;i++){
1184 minFPos.push_back(tempminFPos);
1188 //cout << minDiff << '\t' << minCount << '\t' << endl;
1189 if(minDiff > bdiffs) { success = minDiff; } //no good matches
1191 //check for reverse match
1192 if (alignment != NULL) { delete alignment; }
1194 if (irbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxRBarcodeLength+bdiffs+1)); }
1195 else{ alignment = NULL; }
1197 //can you find the barcode
1200 vector< vector<int> > minRGroup;
1201 vector<int> minRPos;
1203 string rawRSequence = reverseOligo(seq.getUnaligned());
1204 //cout << irbarcodes.size() << '\t' << maxRBarcodeLength << endl;
1205 for(map<string, vector<int> >::iterator it=irbarcodes.begin();it!=irbarcodes.end();it++){
1206 string oligo = reverseOligo(it->first);
1207 //cout << "r before = " << reverseOligo(oligo) << '\t' << reverseOligo(rawRSequence.substr(0,oligo.length()+bdiffs)) << endl;
1208 if(rawRSequence.length() < maxRBarcodeLength){ //let's just assume that the barcodes are the same length
1209 success = bdiffs + 10;
1213 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1214 alignment->alignPrimer(oligo, rawRSequence.substr(0,oligo.length()+bdiffs));
1215 oligo = alignment->getSeqAAln();
1216 string temp = alignment->getSeqBAln();
1218 int alnLength = oligo.length();
1219 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1220 oligo = oligo.substr(0,alnLength);
1221 temp = temp.substr(0,alnLength);
1222 int numDiff = countDiffs(oligo, temp);
1223 if (alnLength == 0) { numDiff = bdiffs + 100; }
1225 //cout << "r after = " << reverseOligo(oligo) << '\t' << reverseOligo(temp) << '\t' << numDiff << endl;
1226 if(numDiff < minDiff){
1230 minRGroup.push_back(it->second);
1231 int tempminRPos = 0;
1233 for(int i=0;i<alnLength;i++){
1238 minRPos.push_back(tempminRPos);
1239 }else if(numDiff == minDiff){
1240 int tempminRPos = 0;
1241 for(int i=0;i<alnLength;i++){
1246 minRPos.push_back(tempminRPos);
1247 minRGroup.push_back(it->second);
1252 if(minDiff > bdiffs) { success = minDiff; } //no good matches
1254 bool foundMatch = false;
1255 vector<int> matches;
1256 for (int i = 0; i < minFGroup.size(); i++) {
1257 for (int j = 0; j < minFGroup[i].size(); j++) {
1258 for (int k = 0; k < minRGroup.size(); k++) {
1259 if (m->inUsersGroups(minFGroup[i][j], minRGroup[k])) { matches.push_back(minFGroup[i][j]); k+= minRGroup.size(); }
1266 if (matches.size() == 1) {
1269 fStart = minFPos[0];
1270 rStart = rawSeq.length() - minRPos[0];
1271 if (fStart > rStart) { foundMatch = false; } //only barcodes not a good sequence
1274 //we have an acceptable match for the forward and reverse, but do they match?
1276 string trimmedSeq = rawSeq.substr(0, rStart); //trim reverse barcode
1277 seq.setUnaligned(trimmedSeq.substr(fStart)); //trim forward barcode
1279 //cout << "barcode = " << ipbarcodes[group].forward << '\t' << ipbarcodes[group].reverse << endl;
1280 }else { success = bdiffs + 100; }
1284 if (alignment != NULL) { delete alignment; }
1290 catch(exception& e) {
1291 m->errorOut(e, "TrimOligos", "stripPairedBarcode");
1296 //*******************************************************************/
1297 int TrimOligos::stripPairedPrimers(Sequence& seq, QualityScores& qual, int& group, bool keepForward){
1300 string rawSeq = seq.getUnaligned();
1301 int success = pdiffs + 1; //guilty until proven innocent
1302 //cout << seq.getName() << endl;
1303 //can you find the forward
1304 for(map<int,oligosPair>::iterator it=ipprimers.begin();it!=ipprimers.end();it++){
1305 string foligo = it->second.forward;
1306 string roligo = it->second.reverse;
1308 //cout << it->first << '\t' << foligo << '\t' << roligo << endl;
1309 //cout << it->first << '\t' << rawSeq.substr(0,foligo.length()) << '\t' << rawSeq.substr(rawSeq.length()-roligo.length(),roligo.length()) << endl;
1310 if(rawSeq.length() < foligo.length()){ //let's just assume that the barcodes are the same length
1311 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
1314 if(rawSeq.length() < roligo.length()){ //let's just assume that the barcodes are the same length
1315 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
1319 if (rawSeq.length() < (foligo.length() + roligo.length())) { success = pdiffs + 10; break; }
1321 if((compareDNASeq(foligo, rawSeq.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawSeq.substr(rawSeq.length()-roligo.length(),roligo.length())))) {
1324 string trimmedSeq = rawSeq.substr(foligo.length()); //trim forward barcode
1325 seq.setUnaligned(trimmedSeq.substr(0,(trimmedSeq.length()-roligo.length()))); //trim reverse barcode
1326 if(qual.getName() != ""){
1327 qual.trimQScores(-1, rawSeq.length()-roligo.length());
1328 qual.trimQScores(foligo.length(), -1);
1335 //cout << "success=" << success << endl;
1336 //if you found the barcode or if you don't want to allow for diffs
1337 if ((pdiffs == 0) || (success == 0)) { return success; }
1338 else { //try aligning and see if you can find it
1339 Alignment* alignment;
1340 if (ifprimers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFPrimerLength+pdiffs+1)); }
1341 else{ alignment = NULL; }
1343 //can you find the barcode
1346 vector< vector<int> > minFGroup;
1347 vector<int> minFPos;
1349 //the pair can have duplicates, but we only want to search for the unique forward and reverses, example
1358 only want to look for forward = Sarah, John, Anna, Pat, Gail
1359 reverse = Westcott, Schloss, Brown, Moore
1360 but if best match forward = 4, and reverse = 1, we want to count as a valid match because forward 1 and forward 4 are the same. so both barcodes map to same group.
1362 //cout << endl << forwardSeq.getName() << endl;
1363 for(map<string, vector<int> >::iterator it=ifprimers.begin();it!=ifprimers.end();it++){
1364 string oligo = it->first;
1366 if(rawSeq.length() < maxFPrimerLength){ //let's just assume that the barcodes are the same length
1367 success = pdiffs + 10;
1370 //cout << "before = " << oligo << '\t' << rawSeq.substr(0,oligo.length()+pdiffs) << endl;
1371 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1372 alignment->alignPrimer(oligo, rawSeq.substr(0,oligo.length()+pdiffs));
1373 oligo = alignment->getSeqAAln();
1374 string temp = alignment->getSeqBAln();
1376 int alnLength = oligo.length();
1378 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1379 oligo = oligo.substr(0,alnLength);
1380 temp = temp.substr(0,alnLength);
1381 int numDiff = countDiffs(oligo, temp);
1383 if (alnLength == 0) { numDiff = pdiffs + 100; }
1384 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
1386 if(numDiff < minDiff){
1390 minFGroup.push_back(it->second);
1391 int tempminFPos = 0;
1393 for(int i=0;i<alnLength;i++){
1398 minFPos.push_back(tempminFPos);
1399 }else if(numDiff == minDiff){
1400 minFGroup.push_back(it->second);
1401 int tempminFPos = 0;
1402 for(int i=0;i<alnLength;i++){
1407 minFPos.push_back(tempminFPos);
1411 //cout << minDiff << '\t' << minCount << '\t' << endl;
1412 if(minDiff > pdiffs) { success = minDiff; } //no good matches
1414 //check for reverse match
1415 if (alignment != NULL) { delete alignment; }
1417 if (irprimers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxRPrimerLength+pdiffs+1)); }
1418 else{ alignment = NULL; }
1420 //can you find the barcode
1423 vector< vector<int> > minRGroup;
1424 vector<int> minRPos;
1426 string rawRSequence = reverseOligo(seq.getUnaligned());
1428 for(map<string, vector<int> >::iterator it=irprimers.begin();it!=irprimers.end();it++){
1429 string oligo = reverseOligo(it->first);
1430 //cout << "r before = " << reverseOligo(oligo) << '\t' << reverseOligo(rawRSequence.substr(0,oligo.length()+pdiffs)) << endl;
1431 if(rawRSequence.length() < maxRPrimerLength){ //let's just assume that the barcodes are the same length
1432 success = pdiffs + 10;
1436 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1437 alignment->alignPrimer(oligo, rawRSequence.substr(0,oligo.length()+pdiffs));
1438 oligo = alignment->getSeqAAln();
1439 string temp = alignment->getSeqBAln();
1441 int alnLength = oligo.length();
1442 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1443 oligo = oligo.substr(0,alnLength);
1444 temp = temp.substr(0,alnLength);
1445 int numDiff = countDiffs(oligo, temp);
1446 if (alnLength == 0) { numDiff = pdiffs + 100; }
1448 //cout << "r after = " << reverseOligo(oligo) << '\t' << reverseOligo(temp) << '\t' << numDiff << endl;
1449 if(numDiff < minDiff){
1453 minRGroup.push_back(it->second);
1454 int tempminRPos = 0;
1456 for(int i=0;i<alnLength;i++){
1461 minRPos.push_back(tempminRPos);
1462 }else if(numDiff == minDiff){
1463 int tempminRPos = 0;
1464 for(int i=0;i<alnLength;i++){
1469 minRPos.push_back(tempminRPos);
1470 minRGroup.push_back(it->second);
1475 if(minDiff > pdiffs) { success = minDiff; } //no good matches
1477 bool foundMatch = false;
1478 vector<int> matches;
1479 for (int i = 0; i < minFGroup.size(); i++) {
1480 for (int j = 0; j < minFGroup[i].size(); j++) {
1481 for (int k = 0; k < minRGroup.size(); k++) {
1482 if (m->inUsersGroups(minFGroup[i][j], minRGroup[k])) { matches.push_back(minFGroup[i][j]); k+= minRGroup.size(); }
1489 if (matches.size() == 1) {
1492 fStart = minFPos[0];
1493 rStart = rawSeq.length() - minRPos[0];
1494 if (fStart > rStart) { foundMatch = false; } //only primers not a good sequence
1497 //we have an acceptable match for the forward and reverse, but do they match?
1500 string trimmedSeq = rawSeq.substr(0, rStart); //trim reverse barcode
1501 seq.setUnaligned(trimmedSeq.substr(fStart)); //trim forward barcode
1502 if(qual.getName() != ""){
1503 qual.trimQScores(-1, rStart);
1504 qual.trimQScores(fStart, -1);
1508 //cout << "primer = " << ipprimers[group].forward << '\t' << ipprimers[group].reverse << endl;
1509 }else { success = pdiffs + 100; }
1513 if (alignment != NULL) { delete alignment; }
1519 catch(exception& e) {
1520 m->errorOut(e, "TrimOligos", "stripPairedPrimers");
1525 //*******************************************************************/
1526 int TrimOligos::stripPairedPrimers(Sequence& seq, int& group){
1529 string rawSeq = seq.getUnaligned();
1530 int success = pdiffs + 1; //guilty until proven innocent
1531 //cout << seq.getName() << endl;
1532 //can you find the forward
1533 for(map<int,oligosPair>::iterator it=ipprimers.begin();it!=ipprimers.end();it++){
1534 string foligo = it->second.forward;
1535 string roligo = it->second.reverse;
1537 //cout << it->first << '\t' << foligo << '\t' << roligo << endl;
1538 //cout << it->first << '\t' << rawSeq.substr(0,foligo.length()) << '\t' << rawSeq.substr(rawSeq.length()-roligo.length(),roligo.length()) << endl;
1539 if(rawSeq.length() < foligo.length()){ //let's just assume that the barcodes are the same length
1540 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
1543 if(rawSeq.length() < roligo.length()){ //let's just assume that the barcodes are the same length
1544 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
1548 if (rawSeq.length() < (foligo.length() + roligo.length())) { success = pdiffs + 10; break; }
1550 if((compareDNASeq(foligo, rawSeq.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawSeq.substr(rawSeq.length()-roligo.length(),roligo.length())))) {
1552 string trimmedSeq = rawSeq.substr(foligo.length()); //trim forward barcode
1553 seq.setUnaligned(trimmedSeq.substr(0,(trimmedSeq.length()-roligo.length()))); //trim reverse barcode
1558 //cout << "success=" << success << endl;
1559 //if you found the barcode or if you don't want to allow for diffs
1560 if ((pdiffs == 0) || (success == 0)) { return success; }
1561 else { //try aligning and see if you can find it
1562 Alignment* alignment;
1563 if (ifprimers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFPrimerLength+pdiffs+1)); }
1564 else{ alignment = NULL; }
1566 //can you find the barcode
1569 vector< vector<int> > minFGroup;
1570 vector<int> minFPos;
1572 //the pair can have duplicates, but we only want to search for the unique forward and reverses, example
1581 only want to look for forward = Sarah, John, Anna, Pat, Gail
1582 reverse = Westcott, Schloss, Brown, Moore
1583 but if best match forward = 4, and reverse = 1, we want to count as a valid match because forward 1 and forward 4 are the same. so both barcodes map to same group.
1585 //cout << endl << forwardSeq.getName() << endl;
1586 for(map<string, vector<int> >::iterator it=ifprimers.begin();it!=ifprimers.end();it++){
1587 string oligo = it->first;
1589 if(rawSeq.length() < maxFPrimerLength){ //let's just assume that the barcodes are the same length
1590 success = pdiffs + 10;
1593 //cout << "before = " << oligo << '\t' << rawSeq.substr(0,oligo.length()+pdiffs) << endl;
1594 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1595 alignment->alignPrimer(oligo, rawSeq.substr(0,oligo.length()+pdiffs));
1596 oligo = alignment->getSeqAAln();
1597 string temp = alignment->getSeqBAln();
1599 int alnLength = oligo.length();
1601 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1602 oligo = oligo.substr(0,alnLength);
1603 temp = temp.substr(0,alnLength);
1604 int numDiff = countDiffs(oligo, temp);
1606 if (alnLength == 0) { numDiff = pdiffs + 100; }
1607 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
1609 if(numDiff < minDiff){
1613 minFGroup.push_back(it->second);
1614 int tempminFPos = 0;
1616 for(int i=0;i<alnLength;i++){
1621 minFPos.push_back(tempminFPos);
1622 }else if(numDiff == minDiff){
1623 minFGroup.push_back(it->second);
1624 int tempminFPos = 0;
1625 for(int i=0;i<alnLength;i++){
1630 minFPos.push_back(tempminFPos);
1634 //cout << minDiff << '\t' << minCount << '\t' << endl;
1635 if(minDiff > pdiffs) { success = minDiff; } //no good matches
1637 //check for reverse match
1638 if (alignment != NULL) { delete alignment; }
1640 if (irprimers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxRPrimerLength+pdiffs+1)); }
1641 else{ alignment = NULL; }
1643 //can you find the barcode
1646 vector< vector<int> > minRGroup;
1647 vector<int> minRPos;
1649 string rawRSequence = reverseOligo(seq.getUnaligned());
1651 for(map<string, vector<int> >::iterator it=irprimers.begin();it!=irprimers.end();it++){
1652 string oligo = reverseOligo(it->first);
1653 //cout << "r before = " << reverseOligo(oligo) << '\t' << reverseOligo(rawRSequence.substr(0,oligo.length()+pdiffs)) << endl;
1654 if(rawRSequence.length() < maxRPrimerLength){ //let's just assume that the barcodes are the same length
1655 success = pdiffs + 10;
1659 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1660 alignment->alignPrimer(oligo, rawRSequence.substr(0,oligo.length()+pdiffs));
1661 oligo = alignment->getSeqAAln();
1662 string temp = alignment->getSeqBAln();
1664 int alnLength = oligo.length();
1665 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1666 oligo = oligo.substr(0,alnLength);
1667 temp = temp.substr(0,alnLength);
1668 int numDiff = countDiffs(oligo, temp);
1669 if (alnLength == 0) { numDiff = pdiffs + 100; }
1671 //cout << "r after = " << reverseOligo(oligo) << '\t' << reverseOligo(temp) << '\t' << numDiff << endl;
1672 if(numDiff < minDiff){
1676 minRGroup.push_back(it->second);
1677 int tempminRPos = 0;
1679 for(int i=0;i<alnLength;i++){
1684 minRPos.push_back(tempminRPos);
1685 }else if(numDiff == minDiff){
1686 int tempminRPos = 0;
1687 for(int i=0;i<alnLength;i++){
1692 minRPos.push_back(tempminRPos);
1693 minRGroup.push_back(it->second);
1698 if(minDiff > pdiffs) { success = minDiff; } //no good matches
1700 bool foundMatch = false;
1701 vector<int> matches;
1702 for (int i = 0; i < minFGroup.size(); i++) {
1703 for (int j = 0; j < minFGroup[i].size(); j++) {
1704 for (int k = 0; k < minRGroup.size(); k++) {
1705 if (m->inUsersGroups(minFGroup[i][j], minRGroup[k])) { matches.push_back(minFGroup[i][j]); k+= minRGroup.size(); }
1712 if (matches.size() == 1) {
1715 fStart = minFPos[0];
1716 rStart = rawSeq.length() - minRPos[0];
1717 if (fStart > rStart) { foundMatch = false; } //only primers not a good sequence
1720 //we have an acceptable match for the forward and reverse, but do they match?
1722 string trimmedSeq = rawSeq.substr(0, rStart); //trim reverse barcode
1723 seq.setUnaligned(trimmedSeq.substr(fStart)); //trim forward barcode
1725 //cout << "primer = " << ipprimers[group].forward << '\t' << ipprimers[group].reverse << endl;
1726 }else { success = pdiffs + 100; }
1730 if (alignment != NULL) { delete alignment; }
1736 catch(exception& e) {
1737 m->errorOut(e, "TrimOligos", "stripPairedPrimers");
1742 //*******************************************************************/
1743 int TrimOligos::stripForward(Sequence& forwardSeq, Sequence& reverseSeq, QualityScores& forwardQual, QualityScores& reverseQual, int& group){
1745 //look for forward barcode
1746 string rawFSequence = forwardSeq.getUnaligned();
1747 string rawRSequence = reverseSeq.getUnaligned();
1748 int success = pdiffs + 1; //guilty until proven innocent
1750 //can you find the forward barcode
1751 for(map<int,oligosPair>::iterator it=ipprimers.begin();it!=ipprimers.end();it++){
1752 string foligo = it->second.forward;
1753 string roligo = it->second.reverse;
1755 if(rawFSequence.length() < foligo.length()){ //let's just assume that the barcodes are the same length
1756 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
1759 if(rawRSequence.length() < roligo.length()){ //let's just assume that the barcodes are the same length
1760 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
1764 if((compareDNASeq(foligo, rawFSequence.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawRSequence.substr(0,roligo.length())))) {
1766 forwardSeq.setUnaligned(rawFSequence.substr(foligo.length()));
1767 reverseSeq.setUnaligned(rawRSequence.substr(roligo.length()));
1768 forwardQual.trimQScores(foligo.length(), -1);
1769 reverseQual.trimQScores(roligo.length(), -1);
1775 //if you found the barcode or if you don't want to allow for diffs
1776 if ((pdiffs == 0) || (success == 0)) { return success; }
1777 else { //try aligning and see if you can find it
1778 Alignment* alignment;
1779 if (ifprimers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFPrimerLength+pdiffs+1)); }
1780 else{ alignment = NULL; }
1782 //can you find the barcode
1785 vector< vector<int> > minFGroup;
1786 vector<int> minFPos;
1788 //the pair can have duplicates, but we only want to search for the unique forward and reverses, example
1797 only want to look for forward = Sarah, John, Anna, Pat, Gail
1798 reverse = Westcott, Schloss, Brown, Moore
1799 but if best match forward = 4, and reverse = 1, we want to count as a valid match because forward 1 and forward 4 are the same. so both barcodes map to same group.
1801 //cout << endl << forwardSeq.getName() << endl;
1802 for(map<string, vector<int> >::iterator it=ifprimers.begin();it!=ifprimers.end();it++){
1803 string oligo = it->first;
1805 if(rawFSequence.length() < maxFPrimerLength){ //let's just assume that the barcodes are the same length
1806 success = pdiffs + 10;
1809 //cout << "before = " << oligo << '\t' << rawFSequence.substr(0,oligo.length()+pdiffs) << endl;
1810 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1811 alignment->alignPrimer(oligo, rawFSequence.substr(0,oligo.length()+pdiffs));
1812 oligo = alignment->getSeqAAln();
1813 string temp = alignment->getSeqBAln();
1815 int alnLength = oligo.length();
1817 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1818 oligo = oligo.substr(0,alnLength);
1819 temp = temp.substr(0,alnLength);
1820 int numDiff = countDiffs(oligo, temp);
1822 if (alnLength == 0) { numDiff = pdiffs + 100; }
1823 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
1825 if (m->debug) { m->mothurOut("[DEBUG]: forward " + forwardSeq.getName() + " aligned fragment=" + temp + ", primer=" + oligo + ", numDiffs=" + toString(numDiff) + ".\n"); }
1827 if(numDiff < minDiff){
1831 minFGroup.push_back(it->second);
1832 int tempminFPos = 0;
1834 for(int i=0;i<alnLength;i++){
1839 minFPos.push_back(tempminFPos);
1840 }else if(numDiff == minDiff){
1841 minFGroup.push_back(it->second);
1842 int tempminFPos = 0;
1843 for(int i=0;i<alnLength;i++){
1848 minFPos.push_back(tempminFPos);
1852 //cout << minDiff << '\t' << minCount << '\t' << endl;
1853 if(minDiff > pdiffs) { success = minDiff; } //no good matches
1855 //check for reverse match
1856 if (alignment != NULL) { delete alignment; }
1858 if (irbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxRPrimerLength+pdiffs+1)); }
1859 else{ alignment = NULL; }
1861 //can you find the barcode
1864 vector< vector<int> > minRGroup;
1865 vector<int> minRPos;
1867 for(map<string, vector<int> >::iterator it=irprimers.begin();it!=irprimers.end();it++){
1868 string oligo = it->first;
1869 //cout << "before = " << oligo << '\t' << rawRSequence.substr(0,oligo.length()+pdiffs) << endl;
1870 if(rawRSequence.length() < maxRPrimerLength){ //let's just assume that the barcodes are the same length
1871 success = pdiffs + 10;
1875 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
1876 alignment->alignPrimer(oligo, rawRSequence.substr(0,oligo.length()+pdiffs));
1877 oligo = alignment->getSeqAAln();
1878 string temp = alignment->getSeqBAln();
1880 int alnLength = oligo.length();
1881 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
1882 oligo = oligo.substr(0,alnLength);
1883 temp = temp.substr(0,alnLength);
1884 int numDiff = countDiffs(oligo, temp);
1885 if (alnLength == 0) { numDiff = pdiffs + 100; }
1887 if (m->debug) { m->mothurOut("[DEBUG]: reverse " + forwardSeq.getName() + " aligned fragment=" + temp + ", primer=" + oligo + ", numDiffs=" + toString(numDiff) + ".\n"); }
1889 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
1890 if(numDiff < minDiff){
1894 minRGroup.push_back(it->second);
1895 int tempminRPos = 0;
1897 for(int i=0;i<alnLength;i++){
1902 minRPos.push_back(tempminRPos);
1903 }else if(numDiff == minDiff){
1904 int tempminRPos = 0;
1905 for(int i=0;i<alnLength;i++){
1910 minRPos.push_back(tempminRPos);
1911 minRGroup.push_back(it->second);
1916 if(minDiff > pdiffs) { success = minDiff; } //no good matches
1918 bool foundMatch = false;
1919 vector<int> matches;
1920 for (int i = 0; i < minFGroup.size(); i++) {
1921 for (int j = 0; j < minFGroup[i].size(); j++) {
1922 for (int k = 0; k < minRGroup.size(); k++) {
1923 if (m->inUsersGroups(minFGroup[i][j], minRGroup[k])) { matches.push_back(minFGroup[i][j]); k+= minRGroup.size(); }
1930 if (matches.size() == 1) {
1933 fStart = minFPos[0];
1934 rStart = minRPos[0];
1937 //we have an acceptable match for the forward and reverse, but do they match?
1939 forwardSeq.setUnaligned(rawFSequence.substr(fStart));
1940 reverseSeq.setUnaligned(rawRSequence.substr(rStart));
1941 forwardQual.trimQScores(fStart, -1);
1942 reverseQual.trimQScores(rStart, -1);
1944 }else { success = pdiffs + 100; }
1948 if (alignment != NULL) { delete alignment; }
1954 catch(exception& e) {
1955 m->errorOut(e, "TrimOligos", "stripIForward");
1960 //*******************************************************************/
1961 int TrimOligos::stripForward(Sequence& forwardSeq, Sequence& reverseSeq, int& group){
1963 //look for forward barcode
1964 string rawFSequence = forwardSeq.getUnaligned();
1965 string rawRSequence = reverseSeq.getUnaligned();
1966 int success = pdiffs + 1; //guilty until proven innocent
1968 //can you find the forward barcode
1969 for(map<int,oligosPair>::iterator it=ipprimers.begin();it!=ipprimers.end();it++){
1970 string foligo = it->second.forward;
1971 string roligo = it->second.reverse;
1973 if(rawFSequence.length() < foligo.length()){ //let's just assume that the barcodes are the same length
1974 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
1977 if(rawRSequence.length() < roligo.length()){ //let's just assume that the barcodes are the same length
1978 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
1982 if((compareDNASeq(foligo, rawFSequence.substr(0,foligo.length()))) && (compareDNASeq(roligo, rawRSequence.substr((rawRSequence.length()-roligo.length()),roligo.length())))) {
1984 forwardSeq.setUnaligned(rawFSequence.substr(foligo.length()));
1985 reverseSeq.setUnaligned(rawRSequence.substr(roligo.length()));
1991 //if you found the barcode or if you don't want to allow for diffs
1992 if ((pdiffs == 0) || (success == 0)) { return success; }
1993 else { //try aligning and see if you can find it
1994 Alignment* alignment;
1995 if (ifprimers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFPrimerLength+pdiffs+1)); }
1996 else{ alignment = NULL; }
1998 //can you find the barcode
2001 vector< vector<int> > minFGroup;
2002 vector<int> minFPos;
2004 //the pair can have duplicates, but we only want to search for the unique forward and reverses, example
2013 only want to look for forward = Sarah, John, Anna, Pat, Gail
2014 reverse = Westcott, Schloss, Brown, Moore
2015 but if best match forward = 4, and reverse = 1, we want to count as a valid match because forward 1 and forward 4 are the same. so both barcodes map to same group.
2017 //cout << endl << forwardSeq.getName() << endl;
2018 for(map<string, vector<int> >::iterator it=ifprimers.begin();it!=ifprimers.end();it++){
2019 string oligo = it->first;
2021 if(rawFSequence.length() < maxFPrimerLength){ //let's just assume that the barcodes are the same length
2022 success = pdiffs + 10;
2025 //cout << "before = " << oligo << '\t' << rawFSequence.substr(0,oligo.length()+pdiffs) << endl;
2026 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2027 alignment->alignPrimer(oligo, rawFSequence.substr(0,oligo.length()+pdiffs));
2028 oligo = alignment->getSeqAAln();
2029 string temp = alignment->getSeqBAln();
2031 int alnLength = oligo.length();
2033 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
2034 oligo = oligo.substr(0,alnLength);
2035 temp = temp.substr(0,alnLength);
2036 int numDiff = countDiffs(oligo, temp);
2038 if (m->debug) { m->mothurOut("[DEBUG]: forward " + forwardSeq.getName() + " aligned fragment=" + temp + ", primer=" + oligo + ", numDiffs=" + toString(numDiff) + ".\n"); }
2040 if (alnLength == 0) { numDiff = pdiffs + 100; }
2041 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
2043 if(numDiff < minDiff){
2047 minFGroup.push_back(it->second);
2048 int tempminFPos = 0;
2050 for(int i=0;i<alnLength;i++){
2055 minFPos.push_back(tempminFPos);
2056 }else if(numDiff == minDiff){
2057 minFGroup.push_back(it->second);
2058 int tempminFPos = 0;
2059 for(int i=0;i<alnLength;i++){
2064 minFPos.push_back(tempminFPos);
2068 //cout << minDiff << '\t' << minCount << '\t' << endl;
2069 if(minDiff > pdiffs) { success = minDiff; } //no good matches
2071 //check for reverse match
2072 if (alignment != NULL) { delete alignment; }
2074 if (irbarcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxRPrimerLength+pdiffs+1)); }
2075 else{ alignment = NULL; }
2077 //can you find the barcode
2080 vector< vector<int> > minRGroup;
2081 vector<int> minRPos;
2083 for(map<string, vector<int> >::iterator it=irprimers.begin();it!=irprimers.end();it++){
2084 string oligo = it->first;
2085 //cout << "before = " << oligo << '\t' << rawRSequence.substr(0,oligo.length()+pdiffs) << endl;
2086 if(rawRSequence.length() < maxRPrimerLength){ //let's just assume that the barcodes are the same length
2087 success = pdiffs + 10;
2091 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2092 alignment->alignPrimer(oligo, rawRSequence.substr(0,oligo.length()+pdiffs));
2093 oligo = alignment->getSeqAAln();
2094 string temp = alignment->getSeqBAln();
2096 int alnLength = oligo.length();
2097 for(int i=oligo.length()-1;i>=0;i--){ if(oligo[i] != '-'){ alnLength = i+1; break; } }
2098 oligo = oligo.substr(0,alnLength);
2099 temp = temp.substr(0,alnLength);
2100 int numDiff = countDiffs(oligo, temp);
2102 if (m->debug) { m->mothurOut("[DEBUG]: reverse " + forwardSeq.getName() + " aligned fragment=" + temp + ", primer=" + oligo + ", numDiffs=" + toString(numDiff) + ".\n"); }
2104 if (alnLength == 0) { numDiff = pdiffs + 100; }
2106 //cout << "after = " << oligo << '\t' << temp << '\t' << numDiff << endl;
2107 if(numDiff < minDiff){
2111 minRGroup.push_back(it->second);
2112 int tempminRPos = 0;
2114 for(int i=0;i<alnLength;i++){
2119 minRPos.push_back(tempminRPos);
2120 }else if(numDiff == minDiff){
2121 int tempminRPos = 0;
2122 for(int i=0;i<alnLength;i++){
2127 minRPos.push_back(tempminRPos);
2128 minRGroup.push_back(it->second);
2133 if(minDiff > pdiffs) { success = minDiff; } //no good matches
2135 bool foundMatch = false;
2136 vector<int> matches;
2137 for (int i = 0; i < minFGroup.size(); i++) {
2138 for (int j = 0; j < minFGroup[i].size(); j++) {
2139 for (int k = 0; k < minRGroup.size(); k++) {
2140 if (m->inUsersGroups(minFGroup[i][j], minRGroup[k])) { matches.push_back(minFGroup[i][j]); k+= minRGroup.size(); }
2147 if (matches.size() == 1) {
2150 fStart = minFPos[0];
2151 rStart = minRPos[0];
2154 //we have an acceptable match for the forward and reverse, but do they match?
2156 forwardSeq.setUnaligned(rawFSequence.substr(fStart));
2157 reverseSeq.setUnaligned(rawRSequence.substr(rStart));
2159 }else { success = pdiffs + 100; }
2163 if (alignment != NULL) { delete alignment; }
2169 catch(exception& e) {
2170 m->errorOut(e, "TrimOligos", "stripIForward");
2175 //*******************************************************************/
2176 int TrimOligos::stripBarcode(Sequence& seq, int& group){
2179 if (paired) { int success = stripPairedBarcode(seq, group); return success; }
2181 string rawSequence = seq.getUnaligned();
2182 int success = bdiffs + 1; //guilty until proven innocent
2184 //can you find the barcode
2185 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
2186 string oligo = it->first;
2187 if(rawSequence.length() < oligo.length()){ //let's just assume that the barcodes are the same length
2188 success = bdiffs + 10; //if the sequence is shorter than the barcode then bail out
2192 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
2194 seq.setUnaligned(rawSequence.substr(oligo.length()));
2201 //if you found the barcode or if you don't want to allow for diffs
2202 if ((bdiffs == 0) || (success == 0)) { return success; }
2204 else { //try aligning and see if you can find it
2205 Alignment* alignment;
2206 if (barcodes.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFBarcodeLength+bdiffs+1)); }
2207 else{ alignment = NULL; }
2209 //can you find the barcode
2215 for(map<string,int>::iterator it=barcodes.begin();it!=barcodes.end();it++){
2216 string oligo = it->first;
2217 // int length = oligo.length();
2219 if(rawSequence.length() < maxFBarcodeLength){ //let's just assume that the barcodes are the same length
2220 success = bdiffs + 10;
2224 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2225 alignment->alignPrimer(oligo, rawSequence.substr(0,oligo.length()+bdiffs));
2226 oligo = alignment->getSeqAAln();
2227 string temp = alignment->getSeqBAln();
2229 int alnLength = oligo.length();
2231 for(int i=oligo.length()-1;i>=0;i--){
2232 if(oligo[i] != '-'){ alnLength = i+1; break; }
2234 oligo = oligo.substr(0,alnLength);
2235 temp = temp.substr(0,alnLength);
2237 int numDiff = countDiffs(oligo, temp);
2239 if (m->debug) { m->mothurOut("[DEBUG]: " + seq.getName() + " aligned fragment =" + temp + ", barcode =" + oligo + ", numDiffs = " + toString(numDiff) + "\n"); }
2241 if(numDiff < minDiff){
2244 minGroup = it->second;
2246 for(int i=0;i<alnLength;i++){
2252 else if(numDiff == minDiff){
2258 if(minDiff > bdiffs) { success = minDiff; } //no good matches
2259 else if(minCount > 1) { success = bdiffs + 100; } //can't tell the difference between multiple barcodes
2260 else{ //use the best match
2262 seq.setUnaligned(rawSequence.substr(minPos));
2266 if (alignment != NULL) { delete alignment; }
2273 catch(exception& e) {
2274 m->errorOut(e, "TrimOligos", "stripBarcode");
2280 /********************************************************************/
2281 int TrimOligos::stripForward(Sequence& seq, int& group){
2284 if (paired) { int success = stripPairedPrimers(seq, group); return success; }
2286 string rawSequence = seq.getUnaligned();
2287 int success = pdiffs + 1; //guilty until proven innocent
2289 //can you find the primer
2290 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
2291 string oligo = it->first;
2292 if(rawSequence.length() < oligo.length()){ //let's just assume that the primers are the same length
2293 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
2297 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
2299 seq.setUnaligned(rawSequence.substr(oligo.length()));
2305 //if you found the barcode or if you don't want to allow for diffs
2306 if ((pdiffs == 0) || (success == 0)) { return success; }
2308 else { //try aligning and see if you can find it
2309 Alignment* alignment;
2310 if (primers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFPrimerLength+pdiffs+1)); }
2311 else{ alignment = NULL; }
2313 //can you find the barcode
2319 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
2320 string oligo = it->first;
2321 // int length = oligo.length();
2323 if(rawSequence.length() < maxFPrimerLength){
2324 success = pdiffs + 100;
2328 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2329 alignment->alignPrimer(oligo, rawSequence.substr(0,oligo.length()+pdiffs));
2330 oligo = alignment->getSeqAAln();
2331 string temp = alignment->getSeqBAln();
2333 int alnLength = oligo.length();
2335 for(int i=oligo.length()-1;i>=0;i--){
2336 if(oligo[i] != '-'){ alnLength = i+1; break; }
2338 oligo = oligo.substr(0,alnLength);
2339 temp = temp.substr(0,alnLength);
2341 int numDiff = countDiffs(oligo, temp);
2343 if(numDiff < minDiff){
2346 minGroup = it->second;
2348 for(int i=0;i<alnLength;i++){
2354 else if(numDiff == minDiff){
2360 if(minDiff > pdiffs) { success = minDiff; } //no good matches
2361 else if(minCount > 1) { success = pdiffs + 10; } //can't tell the difference between multiple primers
2362 else{ //use the best match
2364 seq.setUnaligned(rawSequence.substr(minPos));
2368 if (alignment != NULL) { delete alignment; }
2375 catch(exception& e) {
2376 m->errorOut(e, "TrimOligos", "stripForward");
2380 //*******************************************************************/
2381 int TrimOligos::stripForward(Sequence& seq, QualityScores& qual, int& group, bool keepForward){
2384 if (paired) { int success = stripPairedPrimers(seq, qual, group, keepForward); return success; }
2386 string rawSequence = seq.getUnaligned();
2387 int success = pdiffs + 1; //guilty until proven innocent
2389 //can you find the primer
2390 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
2391 string oligo = it->first;
2392 if(rawSequence.length() < oligo.length()){ //let's just assume that the primers are the same length
2393 success = pdiffs + 10; //if the sequence is shorter than the barcode then bail out
2397 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
2399 if (!keepForward) { seq.setUnaligned(rawSequence.substr(oligo.length())); }
2400 if(qual.getName() != ""){
2401 if (!keepForward) { qual.trimQScores(oligo.length(), -1); }
2408 //if you found the barcode or if you don't want to allow for diffs
2409 if ((pdiffs == 0) || (success == 0)) { return success; }
2411 else { //try aligning and see if you can find it
2412 Alignment* alignment;
2413 if (primers.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxFPrimerLength+pdiffs+1)); }
2414 else{ alignment = NULL; }
2416 //can you find the barcode
2422 for(map<string,int>::iterator it=primers.begin();it!=primers.end();it++){
2423 string oligo = it->first;
2424 // int length = oligo.length();
2426 if(rawSequence.length() < maxFPrimerLength){
2427 success = pdiffs + 100;
2431 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2432 alignment->alignPrimer(oligo, rawSequence.substr(0,oligo.length()+pdiffs));
2433 oligo = alignment->getSeqAAln();
2434 string temp = alignment->getSeqBAln();
2436 int alnLength = oligo.length();
2438 for(int i=oligo.length()-1;i>=0;i--){
2439 if(oligo[i] != '-'){ alnLength = i+1; break; }
2441 oligo = oligo.substr(0,alnLength);
2442 temp = temp.substr(0,alnLength);
2444 int numDiff = countDiffs(oligo, temp);
2446 if(numDiff < minDiff){
2449 minGroup = it->second;
2451 for(int i=0;i<alnLength;i++){
2457 else if(numDiff == minDiff){
2463 if(minDiff > pdiffs) { success = minDiff; } //no good matches
2464 else if(minCount > 1) { success = pdiffs + 10; } //can't tell the difference between multiple primers
2465 else{ //use the best match
2467 if (!keepForward) { seq.setUnaligned(rawSequence.substr(minPos)); }
2468 if(qual.getName() != ""){
2469 if (!keepForward) { qual.trimQScores(minPos, -1); }
2474 if (alignment != NULL) { delete alignment; }
2481 catch(exception& e) {
2482 m->errorOut(e, "TrimOligos", "stripForward");
2486 //******************************************************************/
2487 bool TrimOligos::stripReverse(Sequence& seq, QualityScores& qual){
2489 string rawSequence = seq.getUnaligned();
2490 bool success = 0; //guilty until proven innocent
2492 for(int i=0;i<revPrimer.size();i++){
2493 string oligo = revPrimer[i];
2495 if(rawSequence.length() < oligo.length()){
2500 if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
2501 seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
2502 if(qual.getName() != ""){
2503 qual.trimQScores(-1, rawSequence.length()-oligo.length());
2512 catch(exception& e) {
2513 m->errorOut(e, "TrimOligos", "stripReverse");
2517 //******************************************************************/
2518 bool TrimOligos::stripReverse(Sequence& seq){
2521 string rawSequence = seq.getUnaligned();
2522 bool success = 0; //guilty until proven innocent
2524 for(int i=0;i<revPrimer.size();i++){
2525 string oligo = revPrimer[i];
2527 if(rawSequence.length() < oligo.length()){
2532 if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
2533 seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
2542 catch(exception& e) {
2543 m->errorOut(e, "TrimOligos", "stripReverse");
2547 //******************************************************************/
2548 bool TrimOligos::stripLinker(Sequence& seq, QualityScores& qual){
2550 string rawSequence = seq.getUnaligned();
2551 bool success = ldiffs + 1; //guilty until proven innocent
2553 for(int i=0;i<linker.size();i++){
2554 string oligo = linker[i];
2556 if(rawSequence.length() < oligo.length()){
2557 success = ldiffs + 10;
2561 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
2562 seq.setUnaligned(rawSequence.substr(oligo.length()));
2563 if(qual.getName() != ""){
2564 qual.trimQScores(oligo.length(), -1);
2571 //if you found the linker or if you don't want to allow for diffs
2572 if ((ldiffs == 0) || (success == 0)) { return success; }
2574 else { //try aligning and see if you can find it
2575 Alignment* alignment;
2576 if (linker.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLinkerLength+ldiffs+1)); }
2577 else{ alignment = NULL; }
2579 //can you find the barcode
2584 for(int i = 0; i < linker.size(); i++){
2585 string oligo = linker[i];
2586 // int length = oligo.length();
2588 if(rawSequence.length() < maxLinkerLength){ //let's just assume that the barcodes are the same length
2589 success = ldiffs + 10;
2593 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2594 alignment->alignPrimer(oligo, rawSequence.substr(0,oligo.length()+ldiffs));
2595 oligo = alignment->getSeqAAln();
2596 string temp = alignment->getSeqBAln();
2598 int alnLength = oligo.length();
2600 for(int i=oligo.length()-1;i>=0;i--){
2601 if(oligo[i] != '-'){ alnLength = i+1; break; }
2603 oligo = oligo.substr(0,alnLength);
2604 temp = temp.substr(0,alnLength);
2606 int numDiff = countDiffs(oligo, temp);
2608 if(numDiff < minDiff){
2612 for(int i=0;i<alnLength;i++){
2618 else if(numDiff == minDiff){
2624 if(minDiff > ldiffs) { success = minDiff; } //no good matches
2625 else if(minCount > 1) { success = ldiffs + 100; } //can't tell the difference between multiple barcodes
2626 else{ //use the best match
2627 seq.setUnaligned(rawSequence.substr(minPos));
2629 if(qual.getName() != ""){
2630 qual.trimQScores(minPos, -1);
2635 if (alignment != NULL) { delete alignment; }
2643 catch(exception& e) {
2644 m->errorOut(e, "TrimOligos", "stripLinker");
2648 //******************************************************************/
2649 bool TrimOligos::stripLinker(Sequence& seq){
2652 string rawSequence = seq.getUnaligned();
2653 bool success = ldiffs +1; //guilty until proven innocent
2655 for(int i=0;i<linker.size();i++){
2656 string oligo = linker[i];
2658 if(rawSequence.length() < oligo.length()){
2659 success = ldiffs +10;
2663 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
2664 seq.setUnaligned(rawSequence.substr(oligo.length()));
2670 //if you found the linker or if you don't want to allow for diffs
2671 if ((ldiffs == 0) || (success == 0)) { return success; }
2673 else { //try aligning and see if you can find it
2674 Alignment* alignment;
2675 if (linker.size() > 0) {alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLinkerLength+ldiffs+1)); }
2676 else{ alignment = NULL; }
2678 //can you find the barcode
2683 for(int i = 0; i < linker.size(); i++){
2684 string oligo = linker[i];
2685 // int length = oligo.length();
2687 if(rawSequence.length() < maxLinkerLength){ //let's just assume that the barcodes are the same length
2688 success = ldiffs + 10;
2692 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2693 alignment->alignPrimer(oligo, rawSequence.substr(0,oligo.length()+ldiffs));
2694 oligo = alignment->getSeqAAln();
2695 string temp = alignment->getSeqBAln();
2697 int alnLength = oligo.length();
2699 for(int i=oligo.length()-1;i>=0;i--){
2700 if(oligo[i] != '-'){ alnLength = i+1; break; }
2702 oligo = oligo.substr(0,alnLength);
2703 temp = temp.substr(0,alnLength);
2705 int numDiff = countDiffs(oligo, temp);
2707 if(numDiff < minDiff){
2711 for(int i=0;i<alnLength;i++){
2717 else if(numDiff == minDiff){
2723 if(minDiff > ldiffs) { success = minDiff; } //no good matches
2724 else if(minCount > 1) { success = ldiffs + 100; } //can't tell the difference between multiple barcodes
2725 else{ //use the best match
2726 seq.setUnaligned(rawSequence.substr(minPos));
2730 if (alignment != NULL) { delete alignment; }
2737 catch(exception& e) {
2738 m->errorOut(e, "TrimOligos", "stripLinker");
2743 //******************************************************************/
2744 bool TrimOligos::stripSpacer(Sequence& seq, QualityScores& qual){
2746 string rawSequence = seq.getUnaligned();
2747 bool success = sdiffs+1; //guilty until proven innocent
2749 for(int i=0;i<spacer.size();i++){
2750 string oligo = spacer[i];
2752 if(rawSequence.length() < oligo.length()){
2753 success = sdiffs+10;
2757 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
2758 seq.setUnaligned(rawSequence.substr(oligo.length()));
2759 if(qual.getName() != ""){
2760 qual.trimQScores(oligo.length(), -1);
2767 //if you found the spacer or if you don't want to allow for diffs
2768 if ((sdiffs == 0) || (success == 0)) { return success; }
2770 else { //try aligning and see if you can find it
2771 Alignment* alignment;
2772 if (spacer.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxSpacerLength+sdiffs+1)); }
2773 else{ alignment = NULL; }
2775 //can you find the barcode
2780 for(int i = 0; i < spacer.size(); i++){
2781 string oligo = spacer[i];
2782 // int length = oligo.length();
2784 if(rawSequence.length() < maxSpacerLength){ //let's just assume that the barcodes are the same length
2785 success = sdiffs + 10;
2789 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2790 alignment->alignPrimer(oligo, rawSequence.substr(0,oligo.length()+sdiffs));
2791 oligo = alignment->getSeqAAln();
2792 string temp = alignment->getSeqBAln();
2794 int alnLength = oligo.length();
2796 for(int i=oligo.length()-1;i>=0;i--){
2797 if(oligo[i] != '-'){ alnLength = i+1; break; }
2799 oligo = oligo.substr(0,alnLength);
2800 temp = temp.substr(0,alnLength);
2802 int numDiff = countDiffs(oligo, temp);
2804 if(numDiff < minDiff){
2808 for(int i=0;i<alnLength;i++){
2814 else if(numDiff == minDiff){
2820 if(minDiff > sdiffs) { success = minDiff; } //no good matches
2821 else if(minCount > 1) { success = sdiffs + 100; } //can't tell the difference between multiple barcodes
2822 else{ //use the best match
2823 seq.setUnaligned(rawSequence.substr(minPos));
2825 if(qual.getName() != ""){
2826 qual.trimQScores(minPos, -1);
2831 if (alignment != NULL) { delete alignment; }
2839 catch(exception& e) {
2840 m->errorOut(e, "TrimOligos", "stripSpacer");
2844 //******************************************************************/
2845 bool TrimOligos::stripSpacer(Sequence& seq){
2848 string rawSequence = seq.getUnaligned();
2849 bool success = sdiffs+1; //guilty until proven innocent
2851 for(int i=0;i<spacer.size();i++){
2852 string oligo = spacer[i];
2854 if(rawSequence.length() < oligo.length()){
2855 success = sdiffs+10;
2859 if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
2860 seq.setUnaligned(rawSequence.substr(oligo.length()));
2866 //if you found the spacer or if you don't want to allow for diffs
2867 if ((sdiffs == 0) || (success == 0)) { return success; }
2869 else { //try aligning and see if you can find it
2870 Alignment* alignment;
2871 if (spacer.size() > 0) { alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxSpacerLength+sdiffs+1)); }
2872 else{ alignment = NULL; }
2874 //can you find the barcode
2879 for(int i = 0; i < spacer.size(); i++){
2880 string oligo = spacer[i];
2881 // int length = oligo.length();
2883 if(rawSequence.length() < maxSpacerLength){ //let's just assume that the barcodes are the same length
2884 success = sdiffs + 10;
2888 //use needleman to align first barcode.length()+numdiffs of sequence to each barcode
2889 alignment->alignPrimer(oligo, rawSequence.substr(0,oligo.length()+sdiffs));
2890 oligo = alignment->getSeqAAln();
2891 string temp = alignment->getSeqBAln();
2893 int alnLength = oligo.length();
2895 for(int i=oligo.length()-1;i>=0;i--){
2896 if(oligo[i] != '-'){ alnLength = i+1; break; }
2898 oligo = oligo.substr(0,alnLength);
2899 temp = temp.substr(0,alnLength);
2901 int numDiff = countDiffs(oligo, temp);
2903 if(numDiff < minDiff){
2907 for(int i=0;i<alnLength;i++){
2913 else if(numDiff == minDiff){
2919 if(minDiff > sdiffs) { success = minDiff; } //no good matches
2920 else if(minCount > 1) { success = sdiffs + 100; } //can't tell the difference between multiple barcodes
2921 else{ //use the best match
2922 seq.setUnaligned(rawSequence.substr(minPos));
2926 if (alignment != NULL) { delete alignment; }
2933 catch(exception& e) {
2934 m->errorOut(e, "TrimOligos", "stripSpacer");
2939 //******************************************************************/
2940 bool TrimOligos::compareDNASeq(string oligo, string seq){
2943 int length = oligo.length();
2945 for(int i=0;i<length;i++){
2947 if(oligo[i] != seq[i]){
2948 if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C') { success = 0; }
2949 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N')) { success = 0; }
2950 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G')) { success = 0; }
2951 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T')) { success = 0; }
2952 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A')) { success = 0; }
2953 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G')) { success = 0; }
2954 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A')) { success = 0; }
2955 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G')) { success = 0; }
2956 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G')) { success = 0; }
2957 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G')) { success = 0; }
2958 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C')) { success = 0; }
2959 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G')) { success = 0; }
2961 if(success == 0) { break; }
2970 catch(exception& e) {
2971 m->errorOut(e, "TrimOligos", "compareDNASeq");
2976 //********************************************************************/
2977 int TrimOligos::countDiffs(string oligo, string seq){
2980 int length = oligo.length();
2983 for(int i=0;i<length;i++){
2985 if(oligo[i] != seq[i]){
2986 if(oligo[i] == 'A' || oligo[i] == 'T' || oligo[i] == 'G' || oligo[i] == 'C' || oligo[i] == '-' || oligo[i] == '.') { countDiffs++; }
2987 else if((oligo[i] == 'N' || oligo[i] == 'I') && (seq[i] == 'N')) { countDiffs++; }
2988 else if(oligo[i] == 'R' && (seq[i] != 'A' && seq[i] != 'G')) { countDiffs++; }
2989 else if(oligo[i] == 'Y' && (seq[i] != 'C' && seq[i] != 'T')) { countDiffs++; }
2990 else if(oligo[i] == 'M' && (seq[i] != 'C' && seq[i] != 'A')) { countDiffs++; }
2991 else if(oligo[i] == 'K' && (seq[i] != 'T' && seq[i] != 'G')) { countDiffs++; }
2992 else if(oligo[i] == 'W' && (seq[i] != 'T' && seq[i] != 'A')) { countDiffs++; }
2993 else if(oligo[i] == 'S' && (seq[i] != 'C' && seq[i] != 'G')) { countDiffs++; }
2994 else if(oligo[i] == 'B' && (seq[i] != 'C' && seq[i] != 'T' && seq[i] != 'G')) { countDiffs++; }
2995 else if(oligo[i] == 'D' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'G')) { countDiffs++; }
2996 else if(oligo[i] == 'H' && (seq[i] != 'A' && seq[i] != 'T' && seq[i] != 'C')) { countDiffs++; }
2997 else if(oligo[i] == 'V' && (seq[i] != 'A' && seq[i] != 'C' && seq[i] != 'G')) { countDiffs++; }
3004 catch(exception& e) {
3005 m->errorOut(e, "TrimOligos", "countDiffs");
3009 //********************************************************************/
3010 string TrimOligos::reverseOligo(string oligo){
3012 string reverse = "";
3014 for(int i=oligo.length()-1;i>=0;i--){
3016 if(oligo[i] == 'A') { reverse += 'T'; }
3017 else if(oligo[i] == 'T'){ reverse += 'A'; }
3018 else if(oligo[i] == 'U'){ reverse += 'A'; }
3020 else if(oligo[i] == 'G'){ reverse += 'C'; }
3021 else if(oligo[i] == 'C'){ reverse += 'G'; }
3023 else if(oligo[i] == 'R'){ reverse += 'Y'; }
3024 else if(oligo[i] == 'Y'){ reverse += 'R'; }
3026 else if(oligo[i] == 'M'){ reverse += 'K'; }
3027 else if(oligo[i] == 'K'){ reverse += 'M'; }
3029 else if(oligo[i] == 'W'){ reverse += 'W'; }
3030 else if(oligo[i] == 'S'){ reverse += 'S'; }
3032 else if(oligo[i] == 'B'){ reverse += 'V'; }
3033 else if(oligo[i] == 'V'){ reverse += 'B'; }
3035 else if(oligo[i] == 'D'){ reverse += 'H'; }
3036 else if(oligo[i] == 'H'){ reverse += 'D'; }
3038 else if(oligo[i] == '-'){ reverse += '-'; }
3039 else { reverse += 'N'; }
3045 catch(exception& e) {
3046 m->errorOut(e, "TrimOligos", "reverseOligo");
3051 /********************************************************************/