X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=trimoligos.cpp;h=8c523ce68d8a042bb702900d468def2aadad157b;hb=2009a1a1f47e7467094d844e7c07ab8ddf7bb447;hp=245cfdcac8819d076782b85795300c0dfb034462;hpb=d6c0a11d1cecfac18b323285e7ffadb7f58e848f;p=mothur.git

diff --git a/trimoligos.cpp b/trimoligos.cpp
index 245cfdc..8c523ce 100644
--- a/trimoligos.cpp
+++ b/trimoligos.cpp
@@ -154,7 +154,7 @@ int TrimOligos::stripBarcode(Sequence& seq, QualityScores& qual, int& group){
 			else{													//use the best match
 				group = minGroup;
 				seq.setUnaligned(rawSequence.substr(minPos));
-				
+    
 				if(qual.getName() != ""){
 					qual.trimQScores(minPos, -1);
 				}
@@ -576,25 +576,105 @@ bool TrimOligos::stripReverse(Sequence& seq){
 bool TrimOligos::stripLinker(Sequence& seq, QualityScores& qual){
 	try {
 		string rawSequence = seq.getUnaligned();
-		bool success = 0;	//guilty until proven innocent
+		bool success = ldiffs + 1;	//guilty until proven innocent
 		
 		for(int i=0;i<linker.size();i++){
 			string oligo = linker[i];
 			
 			if(rawSequence.length() < oligo.length()){
-				success = 0;
+				success = ldiffs + 10;
 				break;
 			}
 			
-			if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
-				seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
+			if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
+				seq.setUnaligned(rawSequence.substr(oligo.length()));
 				if(qual.getName() != ""){
-					qual.trimQScores(-1, rawSequence.length()-oligo.length());
+					qual.trimQScores(oligo.length(), -1);
 				}
-				success = 1;
+				success = 0;
 				break;
 			}
-		}	
+		}
+        
+        //if you found the linker or if you don't want to allow for diffs
+		if ((ldiffs == 0) || (success == 0)) { return success;  }
+		
+		else { //try aligning and see if you can find it
+			
+			int maxLength = 0;
+			
+			Alignment* alignment;
+			if (linker.size() > 0) {
+				for(int i = 0; i < linker.size(); i++){
+					if(linker[i].length() > maxLength){
+						maxLength = linker[i].length();
+					}
+				}
+				alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+ldiffs+1));  
+				
+			}else{ alignment = NULL; } 
+			
+			//can you find the barcode
+			int minDiff = 1e6;
+			int minCount = 1;
+			int minPos = 0;
+			
+			for(int i = 0; i < linker.size(); i++){
+				string oligo = linker[i];
+				//				int length = oligo.length();
+				
+				if(rawSequence.length() < maxLength){	//let's just assume that the barcodes are the same length
+					success = ldiffs + 10;
+					break;
+				}
+				
+				//use needleman to align first barcode.length()+numdiffs of sequence to each barcode
+				alignment->align(oligo, rawSequence.substr(0,oligo.length()+ldiffs));
+				oligo = alignment->getSeqAAln();
+				string temp = alignment->getSeqBAln();
+				
+				int alnLength = oligo.length();
+				
+				for(int i=oligo.length()-1;i>=0;i--){
+					if(oligo[i] != '-'){	alnLength = i+1;	break;	}
+				}
+				oligo = oligo.substr(0,alnLength);
+				temp = temp.substr(0,alnLength);
+				
+				int numDiff = countDiffs(oligo, temp);
+				
+				if(numDiff < minDiff){
+					minDiff = numDiff;
+					minCount = 1;
+					minPos = 0;
+					for(int i=0;i<alnLength;i++){
+						if(temp[i] != '-'){
+							minPos++;
+						}
+					}
+				}
+				else if(numDiff == minDiff){
+					minCount++;
+				}
+				
+			}
+			
+			if(minDiff > ldiffs)	{	success = minDiff;		}	//no good matches
+			else if(minCount > 1)	{	success = ldiffs + 100;	}	//can't tell the difference between multiple barcodes
+			else{													//use the best match
+				seq.setUnaligned(rawSequence.substr(minPos));
+				
+				if(qual.getName() != ""){
+					qual.trimQScores(minPos, -1);
+				}
+				success = minDiff;
+			}
+			
+			if (alignment != NULL) {  delete alignment;  }
+			
+		}
+
+       
 		return success;
 		
 	}
@@ -608,23 +688,97 @@ bool TrimOligos::stripLinker(Sequence& seq){
 	try {
 		
 		string rawSequence = seq.getUnaligned();
-		bool success = 0;	//guilty until proven innocent
+		bool success = ldiffs +1;	//guilty until proven innocent
 		
 		for(int i=0;i<linker.size();i++){
 			string oligo = linker[i];
 			
 			if(rawSequence.length() < oligo.length()){
-				success = 0;
+				success = ldiffs +10;
 				break;
 			}
 			
-			if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
-				seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
-				success = 1;
+			if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
+				seq.setUnaligned(rawSequence.substr(oligo.length()));
+				success = 0;
 				break;
 			}
 		}	
 		
+        //if you found the linker or if you don't want to allow for diffs
+		if ((ldiffs == 0) || (success == 0)) { return success;  }
+		
+		else { //try aligning and see if you can find it
+			
+			int maxLength = 0;
+			
+			Alignment* alignment;
+			if (linker.size() > 0) {
+				for(int i = 0; i < linker.size(); i++){
+					if(linker[i].length() > maxLength){
+						maxLength = linker[i].length();
+					}
+				}
+				alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+ldiffs+1));  
+				
+			}else{ alignment = NULL; } 
+			
+			//can you find the barcode
+			int minDiff = 1e6;
+			int minCount = 1;
+			int minPos = 0;
+			
+			for(int i = 0; i < linker.size(); i++){
+				string oligo = linker[i];
+				//				int length = oligo.length();
+				
+				if(rawSequence.length() < maxLength){	//let's just assume that the barcodes are the same length
+					success = ldiffs + 10;
+					break;
+				}
+				
+				//use needleman to align first barcode.length()+numdiffs of sequence to each barcode
+				alignment->align(oligo, rawSequence.substr(0,oligo.length()+ldiffs));
+				oligo = alignment->getSeqAAln();
+				string temp = alignment->getSeqBAln();
+				
+				int alnLength = oligo.length();
+				
+				for(int i=oligo.length()-1;i>=0;i--){
+					if(oligo[i] != '-'){	alnLength = i+1;	break;	}
+				}
+				oligo = oligo.substr(0,alnLength);
+				temp = temp.substr(0,alnLength);
+				
+				int numDiff = countDiffs(oligo, temp);
+				
+				if(numDiff < minDiff){
+					minDiff = numDiff;
+					minCount = 1;
+					minPos = 0;
+					for(int i=0;i<alnLength;i++){
+						if(temp[i] != '-'){
+							minPos++;
+						}
+					}
+				}
+				else if(numDiff == minDiff){
+					minCount++;
+				}
+				
+			}
+			
+			if(minDiff > ldiffs)	{	success = minDiff;		}	//no good matches
+			else if(minCount > 1)	{	success = ldiffs + 100;	}	//can't tell the difference between multiple barcodes
+			else{													//use the best match
+				seq.setUnaligned(rawSequence.substr(minPos));
+				success = minDiff;
+			}
+			
+			if (alignment != NULL) {  delete alignment;  }
+			
+		}
+
 		return success;
 		
 	}
@@ -638,25 +792,105 @@ bool TrimOligos::stripLinker(Sequence& seq){
 bool TrimOligos::stripSpacer(Sequence& seq, QualityScores& qual){
 	try {
 		string rawSequence = seq.getUnaligned();
-		bool success = 0;	//guilty until proven innocent
+		bool success = sdiffs+1;	//guilty until proven innocent
 		
 		for(int i=0;i<spacer.size();i++){
 			string oligo = spacer[i];
 			
 			if(rawSequence.length() < oligo.length()){
-				success = 0;
+				success = sdiffs+10;
 				break;
 			}
 			
-			if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
-				seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
+			if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
+				seq.setUnaligned(rawSequence.substr(oligo.length()));
 				if(qual.getName() != ""){
-					qual.trimQScores(-1, rawSequence.length()-oligo.length());
+					qual.trimQScores(oligo.length(), -1);
 				}
-				success = 1;
+				success = 0;
 				break;
 			}
-		}	
+		}
+        
+        //if you found the spacer or if you don't want to allow for diffs
+		if ((sdiffs == 0) || (success == 0)) { return success;  }
+		
+		else { //try aligning and see if you can find it
+			
+			int maxLength = 0;
+			
+			Alignment* alignment;
+			if (spacer.size() > 0) {
+				for(int i = 0; i < spacer.size(); i++){
+					if(spacer[i].length() > maxLength){
+						maxLength = spacer[i].length();
+					}
+				}
+				alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+sdiffs+1));  
+				
+			}else{ alignment = NULL; } 
+			
+			//can you find the barcode
+			int minDiff = 1e6;
+			int minCount = 1;
+			int minPos = 0;
+			
+			for(int i = 0; i < spacer.size(); i++){
+				string oligo = spacer[i];
+				//				int length = oligo.length();
+				
+				if(rawSequence.length() < maxLength){	//let's just assume that the barcodes are the same length
+					success = sdiffs + 10;
+					break;
+				}
+				
+				//use needleman to align first barcode.length()+numdiffs of sequence to each barcode
+				alignment->align(oligo, rawSequence.substr(0,oligo.length()+sdiffs));
+				oligo = alignment->getSeqAAln();
+				string temp = alignment->getSeqBAln();
+				
+				int alnLength = oligo.length();
+				
+				for(int i=oligo.length()-1;i>=0;i--){
+					if(oligo[i] != '-'){	alnLength = i+1;	break;	}
+				}
+				oligo = oligo.substr(0,alnLength);
+				temp = temp.substr(0,alnLength);
+				
+				int numDiff = countDiffs(oligo, temp);
+				
+				if(numDiff < minDiff){
+					minDiff = numDiff;
+					minCount = 1;
+					minPos = 0;
+					for(int i=0;i<alnLength;i++){
+						if(temp[i] != '-'){
+							minPos++;
+						}
+					}
+				}
+				else if(numDiff == minDiff){
+					minCount++;
+				}
+				
+			}
+			
+			if(minDiff > sdiffs)	{	success = minDiff;		}	//no good matches
+			else if(minCount > 1)	{	success = sdiffs + 100;	}	//can't tell the difference between multiple barcodes
+			else{													//use the best match
+				seq.setUnaligned(rawSequence.substr(minPos));
+				
+				if(qual.getName() != ""){
+					qual.trimQScores(minPos, -1);
+				}
+				success = minDiff;
+			}
+			
+			if (alignment != NULL) {  delete alignment;  }
+			
+		}
+        
+
 		return success;
 		
 	}
@@ -670,23 +904,97 @@ bool TrimOligos::stripSpacer(Sequence& seq){
 	try {
 		
 		string rawSequence = seq.getUnaligned();
-		bool success = 0;	//guilty until proven innocent
+		bool success = sdiffs+1;	//guilty until proven innocent
 		
 		for(int i=0;i<spacer.size();i++){
 			string oligo = spacer[i];
 			
 			if(rawSequence.length() < oligo.length()){
-				success = 0;
+				success = sdiffs+10;
 				break;
 			}
 			
-			if(compareDNASeq(oligo, rawSequence.substr(rawSequence.length()-oligo.length(),oligo.length()))){
-				seq.setUnaligned(rawSequence.substr(0,rawSequence.length()-oligo.length()));
-				success = 1;
+			if(compareDNASeq(oligo, rawSequence.substr(0,oligo.length()))){
+				seq.setUnaligned(rawSequence.substr(oligo.length()));
+				success = 0;
 				break;
 			}
 		}	
 		
+        //if you found the spacer or if you don't want to allow for diffs
+		if ((sdiffs == 0) || (success == 0)) { return success;  }
+		
+		else { //try aligning and see if you can find it
+			
+			int maxLength = 0;
+			
+			Alignment* alignment;
+			if (spacer.size() > 0) {
+				for(int i = 0; i < spacer.size(); i++){
+					if(spacer[i].length() > maxLength){
+						maxLength = spacer[i].length();
+					}
+				}
+				alignment = new NeedlemanOverlap(-1.0, 1.0, -1.0, (maxLength+sdiffs+1));  
+				
+			}else{ alignment = NULL; } 
+			
+			//can you find the barcode
+			int minDiff = 1e6;
+			int minCount = 1;
+			int minPos = 0;
+			
+			for(int i = 0; i < spacer.size(); i++){
+				string oligo = spacer[i];
+				//				int length = oligo.length();
+				
+				if(rawSequence.length() < maxLength){	//let's just assume that the barcodes are the same length
+					success = sdiffs + 10;
+					break;
+				}
+				
+				//use needleman to align first barcode.length()+numdiffs of sequence to each barcode
+				alignment->align(oligo, rawSequence.substr(0,oligo.length()+sdiffs));
+				oligo = alignment->getSeqAAln();
+				string temp = alignment->getSeqBAln();
+				
+				int alnLength = oligo.length();
+				
+				for(int i=oligo.length()-1;i>=0;i--){
+					if(oligo[i] != '-'){	alnLength = i+1;	break;	}
+				}
+				oligo = oligo.substr(0,alnLength);
+				temp = temp.substr(0,alnLength);
+				
+				int numDiff = countDiffs(oligo, temp);
+				
+				if(numDiff < minDiff){
+					minDiff = numDiff;
+					minCount = 1;
+					minPos = 0;
+					for(int i=0;i<alnLength;i++){
+						if(temp[i] != '-'){
+							minPos++;
+						}
+					}
+				}
+				else if(numDiff == minDiff){
+					minCount++;
+				}
+				
+			}
+			
+			if(minDiff > sdiffs)	{	success = minDiff;		}	//no good matches
+			else if(minCount > 1)	{	success = sdiffs + 100;	}	//can't tell the difference between multiple barcodes
+			else{													//use the best match
+				seq.setUnaligned(rawSequence.substr(minPos));
+				success = minDiff;
+			}
+			
+			if (alignment != NULL) {  delete alignment;  }
+			
+		}
+
 		return success;
 		
 	}