X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=alignment.cpp;h=9ab17004bc146bc97168373995c98b932dcb9c3f;hp=95e9dc2b115888a9cd38353aeb67353c1fb71806;hb=050a3ff02473a3d4c0980964e1a9ebe52e55d6b8;hpb=bffbd9ad0d837bc9523d95e7b35c34cfd2631046

diff --git a/alignment.cpp b/alignment.cpp
index 95e9dc2..9ab1700 100644
--- a/alignment.cpp
+++ b/alignment.cpp
@@ -15,77 +15,137 @@
 
 /**************************************************************************************************/
 
-Alignment::Alignment() {	/*	do nothing	*/	}
+Alignment::Alignment() {	m = MothurOut::getInstance(); /*	do nothing	*/	}
 
 /**************************************************************************************************/
 
 Alignment::Alignment(int A) : nCols(A), nRows(A) {
-	
-	alignment.resize(nRows);			//	For the Gotoh and Needleman-Wunsch we initialize the dynamic programming
-	for(int i=0;i<nRows;i++){			//	matrix by initializing a matrix that is A x A.  By default we will set A
-		alignment[i].resize(nCols);		//	at 2000 for 16S rRNA gene sequences
-	}	
-	
+	try {
+ 
+		m = MothurOut::getInstance();
+		alignment.resize(nRows);			//	For the Gotoh and Needleman-Wunsch we initialize the dynamic programming
+		for(int i=0;i<nRows;i++){			//	matrix by initializing a matrix that is A x A.  By default we will set A
+			alignment[i].resize(nCols);		//	at 2000 for 16S rRNA gene sequences
+		}	
+	}
+	catch(exception& e) {
+		m->errorOut(e, "Alignment", "Alignment");
+		exit(1);
+	}
+}
+/**************************************************************************************************/
+void Alignment::resize(int A) {
+	try {
+		nCols = A;
+		nRows = A;
+
+		alignment.resize(nRows);			
+		for(int i=0;i<nRows;i++){			
+			alignment[i].resize(nCols);		
+		}	
+	}
+	catch(exception& e) {
+		m->errorOut(e, "Alignment", "resize");
+		exit(1);
+	}
 }
-
 /**************************************************************************************************/
 
 void Alignment::traceBack(){			//	This traceback routine is used by the dynamic programming algorithms
-										//	to fill the values of seqAaln and seqBaln
-	seqAaln = "";
-	seqBaln = "";
-	int row = lB-1;
-	int column = lA-1;
-//	seqAstart = 1;
-//	seqAend = column;
-	
-	AlignmentCell currentCell = alignment[row][column];	//	Start the traceback from the bottom-right corner of the
-														//	matrix
-
-	if(currentCell.prevCell == 'x'){	seqAaln = seqBaln = "NOALIGNMENT";		}//If there's an 'x' in the bottom-
-	else{												//	right corner bail out because it means nothing got aligned
-		while(currentCell.prevCell != 'x'){				//	while the previous cell isn't an 'x', keep going...
-			
-			if(currentCell.prevCell == 'u'){			//	if the pointer to the previous cell is 'u', go up in the
-				seqAaln = '-' + seqAaln;				//	matrix.  this indicates that we need to insert a gap in
-				seqBaln = seqB[row] + seqBaln;			//	seqA and a base in seqB
-				currentCell = alignment[--row][column];
-			}
-			else if(currentCell.prevCell == 'l'){		//	if the pointer to the previous cell is 'l', go to the left
-				seqBaln = '-' + seqBaln;				//	in the matrix.  this indicates that we need to insert a gap
-				seqAaln = seqA[column] + seqAaln;		//	in seqB and a base in seqA
-				currentCell = alignment[row][--column];
-			}
-			else{
-				seqAaln = seqA[column] + seqAaln;		//	otherwise we need to go diagonally up and to the left,
-				seqBaln = seqB[row] + seqBaln;			//	here we add a base to both alignments
-				currentCell = alignment[--row][--column];
+	try {	
+		BBaseMap.clear();
+        ABaseMap.clear(); //	to fill the values of seqAaln and seqBaln
+		seqAaln = "";
+		seqBaln = "";
+		int row = lB-1;
+		int column = lA-1;
+		//	seqAstart = 1;
+		//	seqAend = column;
+		
+		AlignmentCell currentCell = alignment[row][column];	//	Start the traceback from the bottom-right corner of the
+		//	matrix
+		
+		if(currentCell.prevCell == 'x'){	seqAaln = seqBaln = "NOALIGNMENT";		}//If there's an 'x' in the bottom-
+		else{	//	right corner bail out because it means nothing got aligned
+            int count = 0;
+			while(currentCell.prevCell != 'x'){				//	while the previous cell isn't an 'x', keep going...
+				
+				if(currentCell.prevCell == 'u'){			//	if the pointer to the previous cell is 'u', go up in the
+					seqAaln = '-' + seqAaln;				//	matrix.  this indicates that we need to insert a gap in
+					seqBaln = seqB[row] + seqBaln;			//	seqA and a base in seqB
+                    BBaseMap[row] = count;
+					currentCell = alignment[--row][column];
+				}
+				else if(currentCell.prevCell == 'l'){		//	if the pointer to the previous cell is 'l', go to the left
+					seqBaln = '-' + seqBaln;				//	in the matrix.  this indicates that we need to insert a gap
+					seqAaln = seqA[column] + seqAaln;		//	in seqB and a base in seqA
+                    ABaseMap[column] = count;
+					currentCell = alignment[row][--column];
+				}
+				else{
+					seqAaln = seqA[column] + seqAaln;		//	otherwise we need to go diagonally up and to the left,
+					seqBaln = seqB[row] + seqBaln;			//	here we add a base to both alignments
+                    BBaseMap[row] = count;
+                    ABaseMap[column] = count;
+					currentCell = alignment[--row][--column];
+				}
+                count++;
 			}
 		}
+		
+       
+        pairwiseLength = seqAaln.length();
+		seqAstart = 1;	seqAend = 0;
+		seqBstart = 1;	seqBend = 0;
+        //flip maps since we now know the total length
+        map<int, int> newAMap;
+        for (map<int, int>::iterator it = ABaseMap.begin(); it != ABaseMap.end(); it++) {
+            int spot = it->second;
+            newAMap[pairwiseLength-spot-1] = it->first-1;
+        }
+        ABaseMap = newAMap;
+        map<int, int> newBMap;
+        for (map<int, int>::iterator it = BBaseMap.begin(); it != BBaseMap.end(); it++) {
+            int spot = it->second;
+            newBMap[pairwiseLength-spot-1] = it->first-1;
+        }
+		BBaseMap = newBMap;
+        
+		for(int i=0;i<seqAaln.length();i++){
+			if(seqAaln[i] != '-' && seqBaln[i] == '-')		{	seqAstart++;	}
+			else if(seqAaln[i] == '-' && seqBaln[i] != '-')	{	seqBstart++;	}
+			else											{	break;			}
+		}
+		
+		pairwiseLength -= (seqAstart + seqBstart - 2);
+		
+		for(int i=seqAaln.length()-1; i>=0;i--){
+			if(seqAaln[i] != '-' && seqBaln[i] == '-')		{	seqAend++;		}
+			else if(seqAaln[i] == '-' && seqBaln[i] != '-')	{	seqBend++;		}
+			else											{	break;			}
+		}
+		pairwiseLength -= (seqAend + seqBend);
+		
+		seqAend = seqA.length() - seqAend - 1;
+		seqBend = seqB.length() - seqBend - 1;
 	}
-	
-	pairwiseLength = seqAaln.length();
-	seqAstart = 1;	seqAend = 0;
-	seqBstart = 1;	seqBend = 0;
-	
-	for(int i=0;i<seqAaln.length();i++){
-		if(seqAaln[i] != '-' && seqBaln[i] == '-')		{	seqAstart++;	}
-		else if(seqAaln[i] == '-' && seqBaln[i] != '-')	{	seqBstart++;	}
-		else											{	break;			}
-	}
-	
-	pairwiseLength -= (seqAstart + seqBstart - 2);
-	
-	for(int i=seqAaln.length()-1; i>=0;i--){
-		if(seqAaln[i] != '-' && seqBaln[i] == '-')		{	seqAend++;		}
-		else if(seqAaln[i] == '-' && seqBaln[i] != '-')	{	seqBend++;		}
-		else											{	break;			}
+	catch(exception& e) {
+		m->errorOut(e, "Alignment", "traceBack");
+		exit(1);
 	}
-	pairwiseLength -= (seqAend + seqBend);
-
-	seqAend = seqA.length() - seqAend - 1;
-	seqBend = seqB.length() - seqBend - 1;
+}
+/**************************************************************************************************/
 
+Alignment::~Alignment(){
+	try {
+		for (int i = 0; i < alignment.size(); i++) {
+			for (int j = (alignment[i].size()-1); j >= 0; j--) {  alignment[i].pop_back();  }
+		}
+	}
+	catch(exception& e) {
+		m->errorOut(e, "Alignment", "~Alignment");
+		exit(1);
+	}
 }
 
 /**************************************************************************************************/
@@ -117,7 +177,16 @@ int Alignment::getCandidateEndPos(){
 int Alignment::getTemplateStartPos(){
 	return seqBstart;									//	this is called to report the quality of the alignment
 }
+/**************************************************************************************************/
 
+map<int, int> Alignment::getSeqAAlnBaseMap(){
+	return ABaseMap;									
+}
+/**************************************************************************************************/
+
+map<int, int> Alignment::getSeqBAlnBaseMap(){
+	return BBaseMap;									
+}
 /**************************************************************************************************/
 
 int Alignment::getTemplateEndPos(){