+//***************************************************************************************************************
+bool ClusterFragmentsCommand::isFragment(string seq1, string seq2){
+ try {
+ bool fragment = false;
+
+ //exact match
+ int pos = seq1.find(seq2);
+ if (pos != string::npos) { return true; }
+ //no match, no diffs wanted
+ else if ((diffs == 0) && (percent == 0)) { return false; }
+ else { //try aligning and see if you can find it
+
+ //find number of acceptable differences for this sequence fragment
+ int totalDiffs;
+ if (diffs == 0) { //you didnt set diffs you want a percentage
+ totalDiffs = floor((seq2.length() * (percent / 100.0)));
+ }else if (percent == 0) { //you didn't set percent you want diffs
+ totalDiffs = diffs;
+ }else if ((percent != 0) && (diffs != 0)) { //you want both, set total diffs to smaller of 2
+ totalDiffs = diffs;
+ int percentDiff = floor((seq2.length() * (percent / 100.0)));
+ if (percentDiff < totalDiffs) { totalDiffs = percentDiff; }
+ }
+
+ Alignment* alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (seq1.length()+totalDiffs+1));
+
+ //use needleman to align
+ alignment->align(seq2, seq1);
+ string tempSeq2 = alignment->getSeqAAln();
+ string temp = alignment->getSeqBAln();
+
+ delete alignment;
+
+ //chop gap ends
+ int startPos = 0;
+ int endPos = tempSeq2.length()-1;
+ for (int i = 0; i < tempSeq2.length(); i++) { if (isalpha(tempSeq2[i])) { startPos = i; break; } }
+ for (int i = tempSeq2.length()-1; i >= 0; i--) { if (isalpha(tempSeq2[i])) { endPos = i; break; } }
+
+ //count number of diffs
+ int numDiffs = 0;
+ for (int i = startPos; i <= endPos; i++) {
+ if (tempSeq2[i] != temp[i]) { numDiffs++; }
+ }
+
+ if (numDiffs <= totalDiffs) { fragment = true; }
+
+ }
+
+ return fragment;
+
+ }
+ catch(exception& e) {
+ m->errorOut(e, "ClusterFragmentsCommand", "isFragment");
+ exit(1);
+ }
+}