diff alignCustomAmplicon/utils.cpp @ 5:0aaf65fbb48a draft default tip

Uploaded

--- a/alignCustomAmplicon/utils.cpp	Wed Jan 09 00:45:53 2013 -0500
+++ b/alignCustomAmplicon/utils.cpp	Wed Mar 20 00:22:08 2013 -0400
@@ -6,7 +6,7 @@
 #define min(a,b,c) a < b ? ( a < c ? a : c ) : ( b < c ? b : c ) ;
 using namespace std;
 
-void printMatrix (double* m, int rows, int cols)
+void printMatrix (float* m, int rows, int cols)
 {
 	for(int i = 0; i < rows; i++)
 	{
@@ -44,27 +44,27 @@
 	return d[size1-1][size2-1]; // return edit distance
 }
 
-double nw_align2 (string &seq1, string &seq2, double matchScore, double missmatch_penalty, double gapPenOpen, double gapPenExt, bool noFrontGapPenalty, bool noEndGapPenalty, bool debug)
+float nw_align2 (string &seq1, string &seq2, float matchScore, float missmatch_penalty, float gapPenOpen, float gapPenExt, bool noFrontGapPenalty, bool noEndGapPenalty, bool debug)
 {
 	//Do NW alignment of 2 strings...
 	int size1 = seq1.length();
 	int size2 = seq2.length();
 		
-	double scores[size1+1][size2+1];		
+	float scores[size1+1][size2+1];		
 	char dir[size1+1][size2+1];
 	char gaps[size1+1][size2+1];
 	
 	for(int i = 0; i <= size1; i++) 
 	{
-		scores[i][0] = (double)(i-1)*gapPenExt + gapPenOpen;
+		scores[i][0] = (float)(i-1)*gapPenExt + gapPenOpen;
 		dir[i][0] = 'U';
-		gaps[i][0] = 'O';
+		gaps[i][0] = 'D';
 	}
 	for(int j = 0; j <= size2; j++)
 	{
-		scores[0][j] = (double)(j-1)*gapPenExt + gapPenOpen;
+		scores[0][j] = (float)(j-1)*gapPenExt + gapPenOpen;
 		dir[0][j] = 'L';
-		gaps[0][j] = 'O';
+		gaps[0][j] = 'I';
 	}
 	
 	scores[0][0] = 0;
@@ -76,9 +76,9 @@
 		for(int j = 0; j <= size2; j++) scores[0][j] = 0;
 	}
 			
-	double match;
-	double del;
-	double insert;
+	float match;
+	float del;
+	float insert;
 	
 	for(int i = 1; i <= size1; i++)
 		for(int j = 1; j <= size2; j++)
@@ -89,14 +89,24 @@
 			}
 			else
 			{
-				match = scores[i-1][j-1] + (double)(seq1[i-1] == seq2[j-1]? matchScore : missmatch_penalty);  //1.9f clustalw
+				match = scores[i-1][j-1] + (float)(seq1[i-1] == seq2[j-1]? matchScore : missmatch_penalty);  //1.9f clustalw
 			}
-			del = scores[i-1][j] + (gaps[i-1][j] == 'N'? gapPenOpen:gapPenExt);	
-			insert = scores[i][j-1] + (gaps[i][j-1] == 'N'? gapPenOpen:gapPenExt);
+			del = scores[i-1][j] + (gaps[i-1][j] != 'D'? gapPenOpen:gapPenExt);	
+			insert = scores[i][j-1] + (gaps[i][j-1] != 'I'? gapPenOpen:gapPenExt);
+			
 			//debug
-			//printf("I:%d J:%d match = %f del %f insert:%f\n",i,j,match,del,insert);
+			//if (match == del || match == insert)
+			//{
+				//printf("I:%d J:%d match = %f del %f insert %f %c %c\n",i,j,match,del,insert,seq1[i-1],seq2[j-1]);
+			//}
 			
-			if (match>del && match > insert) 
+			if (match > del && match > insert) 
+			{
+				scores[i][j] = match;
+				dir[i][j] = 'D';
+				gaps[i][j] = 'N';
+			}
+			else if(match >= del && match >= insert && seq1[i-1] != seq2[j-1])
 			{
 				scores[i][j] = match;
 				dir[i][j] = 'D';
@@ -106,13 +116,13 @@
 			{
 				scores[i][j] = del;
 				dir[i][j] = 'U';
-				gaps[i][j] = 'O';
+				gaps[i][j] = 'D';
 			}
 			else 
 			{
 				scores[i][j] = insert;
 				dir[i][j] = 'L';
-				gaps[i][j] = 'O';
+				gaps[i][j] = 'I';
 			}
 		}
 		
@@ -129,7 +139,7 @@
 	int cont2 = size2;
 	int max1 = size1;
 	int max2 = size2;
-	double num = scores[size1][size2];;
+	float num = scores[size1][size2];;
  	
 	if (noEndGapPenalty) 
 	{
@@ -240,7 +250,7 @@
 
 }
 
-vector <string> nw_alignAlignments (vector<string> a1, vector<string> a2, double matchScore, double missmatch_penalty, double gapPenOpen, double gapPenExt, bool noFrontGapPenalty, bool noEndGapPenalty, bool debug)
+vector <string> nw_alignAlignments (vector<string> a1, vector<string> a2, float matchScore, float missmatch_penalty, float gapPenOpen, float gapPenExt, bool noFrontGapPenalty, bool noEndGapPenalty, bool debug)
 {
 	//Do NW alignment of 2 strings...
 	int size1 = a1[0].length();
@@ -260,21 +270,21 @@
 
 
 
-	double scores[size1+1][size2+1];		
+	float scores[size1+1][size2+1];		
 	char dir[size1+1][size2+1];
 	char gaps[size1+1][size2+1];
 	
 	for(int i = 0; i <= size1; i++) 
 	{
-		scores[i][0] = (double)(i-1)*gapPenExt + gapPenOpen;
+		scores[i][0] = (float)(i-1)*gapPenExt + gapPenOpen;
 		dir[i][0] = 'U';
-		gaps[i][0] = 'O';
+		gaps[i][0] = 'D';
 	}
 	for(int j = 0; j <= size2; j++)
 	{
-		scores[0][j] = (double)(j-1)*gapPenExt + gapPenOpen;
+		scores[0][j] = (float)(j-1)*gapPenExt + gapPenOpen;
 		dir[0][j] = 'L';
-		gaps[0][j] = 'O';
+		gaps[0][j] = 'I';
 	}
 	
 	scores[0][0] = 0;
@@ -286,15 +296,17 @@
 		for(int j = 0; j <= size2; j++) scores[0][j] = 0;
 	}
 			
-	double match;
-	double del;
-	double insert;
+	float match;
+	float del;
+	float insert;
+	bool missmatch;
 	
 	for(int i = 1; i <= size1; i++)
 		for(int j = 1; j <= size2; j++)
 		{
 
 			match = del = insert = 0.0f;
+			missmatch = true;
 			for(int z=0;z<v1;++z)
 			{
 				for(int w=0;w<v2;++w)
@@ -305,10 +317,14 @@
 					}
 					else
 					{
-						match += scores[i-1][j-1] + (double)(a1[z][i-1] == a2[w][j-1]? matchScore:missmatch_penalty);
+						match += scores[i-1][j-1] + (float)(a1[z][i-1] == a2[w][j-1]? matchScore:missmatch_penalty);
+						if(a1[z][i-1] == a2[w][j-1])
+						{
+							missmatch = false;
+						}
 					}
-					del += scores[i-1][j] + (gaps[i-1][j] == 'N'? gapPenOpen:gapPenExt);	
-					insert += scores[i][j-1] + (gaps[i][j-1] == 'N'? gapPenOpen:gapPenExt);
+					del += scores[i-1][j] + (gaps[i-1][j] != 'D'? gapPenOpen:gapPenExt);	
+					insert += scores[i][j-1] + (gaps[i][j-1] != 'I'? gapPenOpen:gapPenExt);
 				}			
 			}
 			
@@ -322,17 +338,23 @@
 				dir[i][j] = 'D';
 				gaps[i][j] = 'N';
 			}
+			else if(match >= del && match >= insert && missmatch)
+			{
+				scores[i][j] = match;
+				dir[i][j] = 'D';
+				gaps[i][j] = 'N';
+			}
 			else if(del > insert)
 			{
 				scores[i][j] = del;
 				dir[i][j] = 'U';
-				gaps[i][j] = 'O';
+				gaps[i][j] = 'D';
 			}
 			else 
 			{
 				scores[i][j] = insert;
 				dir[i][j] = 'L';
-				gaps[i][j] = 'O';
+				gaps[i][j] = 'I';
 			}
 		}
 
@@ -346,7 +368,7 @@
 	int cont2 = size2;
 	int max1 = size1;
 	int max2 = size2;
-	double num = scores[size1][size2];
+	float num = scores[size1][size2];
  	
 	if (noEndGapPenalty) 
 	{