X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=needlemanoverlap.cpp;h=1c2ef223539c65a95ca1dac66a823d226b6108b8;hb=a0f1fca79d2ddfa7ad36b4485039c68b5704fe8d;hp=8f5d4649013e003e11cec9110586783c488028d5;hpb=526a868606faa50caf86e7399f7554c0335b39e5;p=mothur.git

diff --git a/needlemanoverlap.cpp b/needlemanoverlap.cpp
index 8f5d464..1c2ef22 100644
--- a/needlemanoverlap.cpp
+++ b/needlemanoverlap.cpp
@@ -17,8 +17,6 @@
  *
  */
 
-using namespace std;
-
 #include "alignmentcell.hpp"
 #include "alignment.hpp"
 #include "overlap.hpp"
@@ -26,20 +24,25 @@ using namespace std;
 
 /**************************************************************************************************/
 
-NeedlemanOverlap::NeedlemanOverlap(float gO, float m, float mm, int r) ://	note that we don't have a gap extend
-gap(gO), match(m), mismatch(mm), Alignment(r) {							//	the gap openning penalty is assessed for
-																		//	every gapped position
-	for(int i=1;i<nCols;i++){
-		alignment[0][i].prevCell = 'l';					//	initialize first row by pointing all poiters to the left
-		alignment[0][i].cValue = 0;						//	and the score to zero
-	}
+NeedlemanOverlap::NeedlemanOverlap(float gO, float f, float mm, int r) ://	note that we don't have a gap extend
+gap(gO), match(f), mismatch(mm), Alignment(r) {							//	the gap openning penalty is assessed for
+	try {																	//	every gapped position
+		for(int i=1;i<nCols;i++){
+			alignment[0][i].prevCell = 'l';					//	initialize first row by pointing all poiters to the left
+			alignment[0][i].cValue = 0;						//	and the score to zero
+		}
+		
+		for(int i=1;i<nRows;i++){
+			alignment[i][0].prevCell = 'u';					//	initialize first column by pointing all poiters upwards
+			alignment[i][0].cValue = 0;						//	and the score to zero
+		}
 	
-	for(int i=1;i<nRows;i++){
-		alignment[i][0].prevCell = 'u';					//	initialize first column by pointing all poiters upwards
-		alignment[i][0].cValue = 0;						//	and the score to zero
+	}
+	catch(exception& e) {
+		m->errorOut(e, "NeedlemanOverlap", "NeedlemanOverlap");
+		exit(1);
 	}
 }
-
 /**************************************************************************************************/
 
 NeedlemanOverlap::~NeedlemanOverlap(){	/*	do nothing	*/	}
@@ -47,46 +50,151 @@ NeedlemanOverlap::~NeedlemanOverlap(){	/*	do nothing	*/	}
 /**************************************************************************************************/
 
 void NeedlemanOverlap::align(string A, string B){
+	try {
 	
-	seqA = ' ' + A;	lA = seqA.length();		//	algorithm requires a dummy space at the beginning of each string
-	seqB = ' ' + B;	lB = seqB.length();		//	algorithm requires a dummy space at the beginning of each string
+		seqA = ' ' + A;	lA = seqA.length();		//	algorithm requires a dummy space at the beginning of each string
+		seqB = ' ' + B;	lB = seqB.length();		//	algorithm requires a dummy space at the beginning of each string
+
+		if (lA > nRows) { m->mothurOut("One of your candidate sequences is longer than you longest template sequence. Your longest template sequence is " + toString(nRows) + ". Your candidate is " + toString(lA) + "."); m->mothurOutEndLine();  }
+		
+		for(int i=1;i<lB;i++){					//	This code was largely translated from Perl code provided in Ex 3.1 
+		
+			for(int j=1;j<lA;j++){				//	of the O'Reilly BLAST book.  I found that the example output had a
 	
-	for(int i=1;i<lB;i++){					//	This code was largely translated from Perl code provided in Ex 3.1 
-		for(int j=1;j<lA;j++){				//	of the O'Reilly BLAST book.  I found that the example output had a
-											//	number of errors
-			float diagonal;
-			if(seqB[i] == seqA[j])	{	diagonal = alignment[i-1][j-1].cValue + match;		}
-			else					{	diagonal = alignment[i-1][j-1].cValue + mismatch;	}
-			
-			float up	= alignment[i-1][j].cValue + gap;
-			float left	= alignment[i][j-1].cValue + gap;
+				//	number of errors
+				float diagonal;
+				if(seqB[i] == seqA[j])	{	diagonal = alignment[i-1][j-1].cValue + match;		}
+				else					{	diagonal = alignment[i-1][j-1].cValue + mismatch;	}
 			
-			if(diagonal >= up){
-				if(diagonal >= left){
-					alignment[i][j].cValue = diagonal;
-					alignment[i][j].prevCell = 'd';
+				float up	= alignment[i-1][j].cValue + gap;
+				float left	= alignment[i][j-1].cValue + gap;
+				
+				if(diagonal >= up){
+					if(diagonal >= left){
+						alignment[i][j].cValue = diagonal;
+						alignment[i][j].prevCell = 'd';
+					}
+					else{
+						alignment[i][j].cValue = left;
+						alignment[i][j].prevCell = 'l';
+					}
 				}
 				else{
-					alignment[i][j].cValue = left;
-					alignment[i][j].prevCell = 'l';
+					if(up >= left){
+						alignment[i][j].cValue = up;
+						alignment[i][j].prevCell = 'u';
+					}
+					else{
+						alignment[i][j].cValue = left;
+						alignment[i][j].prevCell = 'l';
+					}
 				}
 			}
-			else{
-				if(up >= left){
-					alignment[i][j].cValue = up;
-					alignment[i][j].prevCell = 'u';
+		}
+
+		Overlap over;						
+		over.setOverlap(alignment, lA, lB, 0);		//	Fix gaps at the beginning and end of the sequences
+		traceBack();								//	Traceback the alignment to populate seqAaln and seqBaln
+	
+	}
+	catch(exception& e) {
+		m->errorOut(e, "NeedlemanOverlap", "align");
+		exit(1);
+	}
+
+}
+/**************************************************************************************************/
+
+void NeedlemanOverlap::alignPrimer(string A, string B){
+	try {
+        
+		seqA = ' ' + A;	lA = seqA.length();		//	algorithm requires a dummy space at the beginning of each string
+		seqB = ' ' + B;	lB = seqB.length();		//	algorithm requires a dummy space at the beginning of each string
+        
+		if (lA > nRows) { m->mothurOut("One of your candidate sequences is longer than you longest template sequence. Your longest template sequence is " + toString(nRows) + ". Your candidate is " + toString(lA) + "."); m->mothurOutEndLine();  }
+		
+		for(int i=1;i<lB;i++){					//	This code was largely translated from Perl code provided in Ex 3.1
+            
+			for(int j=1;j<lA;j++){				//	of the O'Reilly BLAST book.  I found that the example output had a
+                
+				//	number of errors
+				float diagonal;
+				if(isEquivalent(seqB[i],seqA[j]))	{	diagonal = alignment[i-1][j-1].cValue + match;		}
+				else                                {	diagonal = alignment[i-1][j-1].cValue + mismatch;	}
+                
+				float up	= alignment[i-1][j].cValue + gap;
+				float left	= alignment[i][j-1].cValue + gap;
+				
+				if(diagonal >= up){
+					if(diagonal >= left){
+						alignment[i][j].cValue = diagonal;
+						alignment[i][j].prevCell = 'd';
+					}
+					else{
+						alignment[i][j].cValue = left;
+						alignment[i][j].prevCell = 'l';
+					}
 				}
 				else{
-					alignment[i][j].cValue = left;
-					alignment[i][j].prevCell = 'l';
+					if(up >= left){
+						alignment[i][j].cValue = up;
+						alignment[i][j].prevCell = 'u';
+					}
+					else{
+						alignment[i][j].cValue = left;
+						alignment[i][j].prevCell = 'l';
+					}
 				}
 			}
 		}
+        
+		Overlap over;
+		over.setOverlap(alignment, lA, lB, 0);		//	Fix gaps at the beginning and end of the sequences
+		traceBack();								//	Traceback the alignment to populate seqAaln and seqBaln
+        
+	}
+	catch(exception& e) {
+		m->errorOut(e, "NeedlemanOverlap", "alignPrimer");
+		exit(1);
 	}
-	Overlap over;						
-	over.setOverlap(alignment, lA, lB, 0);		//	Fix gaps at the beginning and end of the sequences
-	traceBack();								//	Traceback the alignment to populate seqAaln and seqBaln
+    
 }
+//********************************************************************/
+bool NeedlemanOverlap::isEquivalent(char oligo, char seq){
+	try {
+		
+        bool same = true;
+					
+        oligo = toupper(oligo);
+        seq = toupper(seq);
+        
+        if(oligo != seq){
+            if(oligo == 'A' && (seq != 'A' && seq != 'M' && seq != 'R' && seq != 'W' && seq != 'D' && seq != 'H' && seq != 'V'))       {	same = false;	}
+            else if(oligo == 'C' && (seq != 'C' && seq != 'Y' && seq != 'M' && seq != 'S' && seq != 'B' && seq != 'H' && seq != 'V'))       {	same = false;	}
+            else if(oligo == 'G' && (seq != 'G' && seq != 'R' && seq != 'K' && seq != 'S' && seq != 'B' && seq != 'D' && seq != 'V'))       {	same = false;	}
+            else if(oligo == 'T' && (seq != 'T' && seq != 'Y' && seq != 'K' && seq != 'W' && seq != 'B' && seq != 'D' && seq != 'H'))       {	same = false;	}
+            else if((oligo == '.' || oligo == '-'))           {	same = false;	}
+            else if((oligo == 'N' || oligo == 'I') && (seq == 'N'))                         {	same = false;	}
+            else if(oligo == 'R' && (seq != 'A' && seq != 'G'))                        {	same = false;	}
+            else if(oligo == 'Y' && (seq != 'C' && seq != 'T'))                        {	same = false;	}
+            else if(oligo == 'M' && (seq != 'C' && seq != 'A'))                        {	same = false;	}
+            else if(oligo == 'K' && (seq != 'T' && seq != 'G'))                        {	same = false;	}
+            else if(oligo == 'W' && (seq != 'T' && seq != 'A'))                        {	same = false;	}
+            else if(oligo == 'S' && (seq != 'C' && seq != 'G'))                        {	same = false;	}
+            else if(oligo == 'B' && (seq != 'C' && seq != 'T' && seq != 'G'))       {	same = false;	}
+            else if(oligo == 'D' && (seq != 'A' && seq != 'T' && seq != 'G'))       {	same = false;	}
+            else if(oligo == 'H' && (seq != 'A' && seq != 'T' && seq != 'C'))       {	same = false;	}
+            else if(oligo == 'V' && (seq != 'A' && seq != 'C' && seq != 'G'))       {	same = false;	}
+        }
 
+		
+		
+		return same;
+	}
+	catch(exception& e) {
+		m->errorOut(e, "TrimOligos", "countDiffs");
+		exit(1);
+	}
+}
 /**************************************************************************************************/