X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=kmer.cpp;h=50574f4a5ee0709daf2ebb67b8f0e49971e20ef0;hp=2c2783384c3849c2941509547d2a67b83dee6f9f;hb=cf9987b67aa49777a4c91c2d21f96e58bf17aa82;hpb=20a2d0350a737a434c89f303662d64a8eeea7b05

diff --git a/kmer.cpp b/kmer.cpp
index 2c27833..50574f4 100644
--- a/kmer.cpp
+++ b/kmer.cpp
@@ -7,47 +7,85 @@
  *
  */
 
-using namespace std;
-
-#include <string>
-#include <vector>
-
 #include "kmer.hpp"
 
 /**************************************************************************************************/
 
-Kmer::Kmer(int size) : kmerSize(size) {
-	
-	int power4s[9] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536 };
-	maxKmer = power4s[kmerSize]+1;// (int)pow(4.,k)+1;
+Kmer::Kmer(int size) : kmerSize(size) {	//	The constructor sets the size of the kmer
 	
+	int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 };
+																		//	No reason to waste the time of recalculating
+	maxKmer = power4s[kmerSize]+1;// (int)pow(4.,k)+1;					//	powers of 4 everytime through.  We need an
+																		//	extra kmer if we get a non-ATGCU base
 }
 
 /**************************************************************************************************/
 
-string Kmer::getKmerString(string sequence){
-	int length = sequence.length();
-	int nKmers = length - kmerSize + 1;
+string Kmer::getKmerString(string sequence){	//	Calculate kmer for each position in the sequence, count the freq
+	int length = sequence.length();				//	of each kmer, and convert it to an ascii character with base '!'.
+	int nKmers = length - kmerSize + 1;			//	Export the string of characters as a string
 	vector<int> counts(maxKmer, 0);
 	
-	for(int i=0;i<nKmers;i++){
-		int kmerNumber = getKmerNumber(sequence, i);
+	for(int i=0;i<nKmers;i++){					//	Go though sequence and get the number between 0 and maxKmer for that
+		int kmerNumber = getKmerNumber(sequence, i);//	kmer.  Increase the count for the kmer in the counts vector
 		counts[kmerNumber]++;
 	}
 	
-	string kmerString = "";
-	for(int i=0;i<maxKmer;i++){
-		kmerString += getASCII(counts[i]);
+	string kmerString = "";						
+	for(int i=0;i<maxKmer;i++){					//	Scan through the vector of counts and convert each element of the 
+		kmerString += getASCII(counts[i]);		//	vector to an ascii character that is equal to or greater than '!'
 	}		
 	
 	return kmerString;	
 }
+/**************************************************************************************************/
+
+vector< map<int, int> > Kmer::getKmerCounts(string sequence){	//	Calculate kmer for each position in the sequence, save info in a map
+	int length = sequence.length();				//	so you know at each spot in the sequence what kmers were found
+	int nKmers = length - kmerSize + 1;			//	
+	vector< map<int, int> > counts; counts.resize(nKmers);  // a map kmer counts for each spot
+	 map<int, int>::iterator it;
+	
+	for(int i=0;i<nKmers;i++){					//	Go though sequence and get the number between 0 and maxKmer for that
+		if (i == 0) {
+			int kmerNumber = getKmerNumber(sequence, i);//	kmer. 
+			counts[i][kmerNumber] = 1;			// add this kmer if not already there
+		}else { 
+			//your count is everything that came before and whatever you find now
+			counts[i] = counts[i-1];
+			int kmerNumber = getKmerNumber(sequence, i);//	kmer.  
+			
+			it = counts[i].find(kmerNumber);
+			if (it!= counts[i].end()) {   counts[i][kmerNumber]++;  }  //increment number of times you have seen this kmer
+			else {   counts[i][kmerNumber] = 1;   }  // add this kmer since not already there
+		}
+	}
+	
+	return counts;	
+}
+	
 	
 /**************************************************************************************************/
 
 int Kmer::getKmerNumber(string sequence, int index){
-	int power4s[9] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536 };
+	
+//	Here we convert a kmer to a number between 0 and maxKmer.  For example, AAAA would equal 0 and TTTT would equal 255.
+//	If there's an N in the kmer, it is set to 256 (if we are looking at 4mers).  The largest we can look at are 8mers,
+//	this could be easily increased.
+//	
+//	Example:   ATGCAT (kSize = 6)
+//		  i:   012345
+//
+//		Score	= 0*4^(6-0-1) + 3*4^(6-1-1) + 2*4^(6-2-1) + 1*4^(6-3-1) + 0*4^(6-4-1) + 3*4^(6-5-1)
+//				= 0*4^5	+	3*4^4	+	2*4^3	+	1*4^2	+	0*4^1	+	3*4^0
+//				= 0 + 3*256 + 2*64 + 1*16 + 0*4 + 3*1
+//				= 0 + 768 + 128 + 16 + 0 + 3
+//				= 915
+	
+	int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 };
+	
 	int kmer = 0;
+	
 	for(int i=0;i<kmerSize;i++){
 		if(toupper(sequence[i+index]) == 'A')		{	kmer += (0 * power4s[kmerSize-i-1]);	}
 		else if(toupper(sequence[i+index]) == 'C')	{	kmer += (1 * power4s[kmerSize-i-1]);	}
@@ -62,19 +100,30 @@ int Kmer::getKmerNumber(string sequence, int index){
 /**************************************************************************************************/
 	
 string Kmer::getKmerBases(int kmerNumber){
-	int power4s[9] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536 };
+	
+//	Here we convert the kmer number into the kmer in terms of bases.
+//
+//	Example:	Score = 915 (for a 6-mer)
+//				Base6 = (915 / 4^0) % 4 = 915 % 4 = 3 => T	[T]
+//				Base5 = (915 / 4^1) % 4 = 228 % 4 = 0 => A	[AT]
+//				Base4 = (915 / 4^2) % 4 = 57 % 4 = 1 => C	[CAT]
+//				Base3 = (915 / 4^3) % 4 = 14 % 4 = 2 => G	[GCAT]
+//				Base2 = (915 / 4^4) % 4 = 3 % 4 = 3 => T	[TGCAT]
+//				Base1 = (915 / 4^5) % 4 = 0 % 4 = 0 => A	[ATGCAT] -> this checks out with the previous method
+	
+	int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 };
 	
 	string kmer = "";
 	
-	if(kmerNumber == power4s[kmerSize]){//pow(4.,7)){	
-		for(int i=0;i<kmerSize;i++){
+	if(kmerNumber == power4s[kmerSize]){//pow(4.,7)){	//	if the kmer number is the same as the maxKmer then it must
+		for(int i=0;i<kmerSize;i++){					//	have had an N in it and so we'll just call it N x kmerSize
 			kmer += 'N';
 		}
 	}
 	else{
 		for(int i=0;i<kmerSize;i++){
-			int nt = (int)(kmerNumber / (float)power4s[i]) % 4;
-			if(nt == 0)		{	kmer = 'A' + kmer;	}
+			int nt = (int)(kmerNumber / (float)power4s[i]) % 4;		//	the '%' operator returns the remainder 
+			if(nt == 0)		{	kmer = 'A' + kmer;	}				//	from int-based division ]
 			else if(nt == 1){	kmer = 'C' + kmer;	}
 			else if(nt == 2){	kmer = 'G' + kmer;	}
 			else if(nt == 3){	kmer = 'T' + kmer;	}
@@ -82,13 +131,34 @@ string Kmer::getKmerBases(int kmerNumber){
 	}
 	return kmer;
 }
-
 /**************************************************************************************************/
 
-char Kmer::getASCII(int number)		{	return (char)(33+number);			}
+int Kmer::getReverseKmerNumber(int kmerNumber){
+		
+	string kmerString = getKmerBases(kmerNumber);
+	
+	//get Reverse
+	string reverse = "";
+	for(int i=kmerString.length()-1;i>=0;i--){
+		if(kmerString[i] == 'A')		{	reverse += 'T';	}
+		else if(kmerString[i] == 'T'){	reverse += 'A';	}
+		else if(kmerString[i] == 'G'){	reverse += 'C';	}
+		else if(kmerString[i] == 'C'){	reverse += 'G';	}
+		else						{	reverse += 'N';	}
+	}
+	
+	int reverseNumber = getKmerNumber(reverse, 0);
+	
+	return reverseNumber;
+	
+}
+
+/**************************************************************************************************/
 
+char Kmer::getASCII(int number)		{	return (char)(33+number);			}	// '!' is the first printable char and
+																				// has the int value of 33
 /**************************************************************************************************/
 
-int Kmer::getNumber(char character)	{	return ((int)(character-'!'));		}
+int Kmer::getNumber(char character)	{	return ((int)(character-'!'));		}	// '!' has the value of 33
 
 /**************************************************************************************************/