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

diff --git a/alignment.cpp b/alignment.cpp
index 25b2760..9ab1700 100644
--- a/alignment.cpp
+++ b/alignment.cpp
@@ -15,28 +15,46 @@
 
 /**************************************************************************************************/
 
-Alignment::Alignment() {	/*	do nothing	*/	}
+Alignment::Alignment() {	m = MothurOut::getInstance(); /*	do nothing	*/	}
 
 /**************************************************************************************************/
 
 Alignment::Alignment(int A) : nCols(A), nRows(A) {
 	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) {
-		errorOut(e, "Alignment", "Alignment");
+		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
 	try {	
-									//	to fill the values of seqAaln and seqBaln
+		BBaseMap.clear();
+        ABaseMap.clear(); //	to fill the values of seqAaln and seqBaln
 		seqAaln = "";
 		seqBaln = "";
 		int row = lB-1;
@@ -48,31 +66,51 @@ void Alignment::traceBack(){			//	This traceback routine is used by the dynamic
 		//	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
+		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();
+       
+        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++;	}
@@ -92,7 +130,7 @@ void Alignment::traceBack(){			//	This traceback routine is used by the dynamic
 		seqBend = seqB.length() - seqBend - 1;
 	}
 	catch(exception& e) {
-		errorOut(e, "Alignment", "traceBack");
+		m->errorOut(e, "Alignment", "traceBack");
 		exit(1);
 	}
 }
@@ -105,7 +143,7 @@ Alignment::~Alignment(){
 		}
 	}
 	catch(exception& e) {
-		errorOut(e, "Alignment", "~Alignment");
+		m->errorOut(e, "Alignment", "~Alignment");
 		exit(1);
 	}
 }
@@ -139,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(){