X-Git-Url: https://git.donarmstrong.com/?p=mothur.git;a=blobdiff_plain;f=kmerdb.cpp;h=9a5d23500be955b459413beee989601f6c0ec78f;hp=0ec822f763feb3869a3afe63fecd78a71dec3e82;hb=b206f634aae1b4ce13978d203247fb64757d5482;hpb=510b1cfc25cd79391d6973ca20c5ec25fb1bb3b2

diff --git a/kmerdb.cpp b/kmerdb.cpp
index 0ec822f..9a5d235 100644
--- a/kmerdb.cpp
+++ b/kmerdb.cpp
@@ -29,124 +29,227 @@
 
 /**************************************************************************************************/
 
-KmerDB::KmerDB(string fastaFileName, int kSize) : Database(fastaFileName), kmerSize(kSize) {
-
-	string kmerDBName = fastaFileName.substr(0,fastaFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
-	ifstream kmerFileTest(kmerDBName.c_str());
-	
-	int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 };
+KmerDB::KmerDB(string fastaFileName, int kSize) : Database(), kmerSize(kSize) {
+	try { 
 	
-	maxKmer = power4s[kmerSize];
-	kmerLocations.resize(maxKmer+1);
-	
-	if(!kmerFileTest){		//	if we can open the kmer db file, then read it in...
-		mothurOut("Generating the " + kmerDBName + " database...\t");	cout.flush();
-		generateKmerDB(kmerDBName);	
-	}
-	else{					//	...otherwise generate it.
-		mothurOut("Reading in the " + kmerDBName + " database...\t");	cout.flush();
-		readKmerDB(kmerDBName, kmerFileTest);
+		kmerDBName = fastaFileName.substr(0,fastaFileName.find_last_of(".")+1) + char('0'+ kmerSize) + "mer";
+		
+		int power4s[14] = { 1, 4, 16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576, 4194304, 16777216, 67108864 };
+		count = 0;
+		
+		maxKmer = power4s[kmerSize];
+		kmerLocations.resize(maxKmer+1);
+		
 	}
-	mothurOut("DONE."); mothurOutEndLine();	mothurOutEndLine(); cout.flush();
+	catch(exception& e) {
+		m->errorOut(e, "KmerDB", "KmerDB");
+		exit(1);
+	}	
 
 }
 /**************************************************************************************************/
+KmerDB::KmerDB() : Database() {}
+/**************************************************************************************************/
 
-KmerDB::~KmerDB(){														
-	
-		//for (int i = 0; i < templateSequences.size(); i++) {  delete templateSequences[i]; }
-		// templateSequences.clear(); 
-}
+KmerDB::~KmerDB(){}
 
 /**************************************************************************************************/
 
-Sequence KmerDB::findClosestSequence(Sequence* candidateSeq){
-	
-	Kmer kmer(kmerSize);
+vector<int> KmerDB::findClosestSequences(Sequence* candidateSeq, int num){
+	try {
+		if (num > numSeqs) { m->mothurOut("[WARNING]: you requested " + toString(num) + " closest sequences, but the template only contains " + toString(numSeqs) + ", adjusting."); m->mothurOutEndLine(); num = numSeqs; }
+		
+		vector<int> topMatches;
+		Kmer kmer(kmerSize);
+		searchScore = 0;
+		Scores.clear();
+		
+		vector<int> matches(numSeqs, 0);						//	a record of the sequences with shared kmers
+		vector<int> timesKmerFound(kmerLocations.size()+1, 0);	//	a record of the kmers that we have already found
+		
+		int numKmers = candidateSeq->getNumBases() - kmerSize + 1;	
 	
-	searchScore = 0;
-	int maxSequence = 0;
-
-	vector<int> matches(numSeqs, 0);						//	a record of the sequences with shared kmers
-	vector<int> timesKmerFound(kmerLocations.size()+1, 0);	//	a record of the kmers that we have already found
-
-	int numKmers = candidateSeq->getNumBases() - kmerSize + 1;	
-
-	for(int i=0;i<numKmers;i++){
-		int kmerNumber = kmer.getKmerNumber(candidateSeq->getUnaligned(), i);		//	go through the query sequence and get a kmer number
-		if(timesKmerFound[kmerNumber] == 0){				//	if we haven't seen it before...
-			for(int j=0;j<kmerLocations[kmerNumber].size();j++){//increase the count for each sequence that also has
-				matches[kmerLocations[kmerNumber][j]]++;	//	that kmer
+		for(int i=0;i<numKmers;i++){
+			int kmerNumber = kmer.getKmerNumber(candidateSeq->getUnaligned(), i);		//	go through the query sequence and get a kmer number
+			if(timesKmerFound[kmerNumber] == 0){				//	if we haven't seen it before...
+				for(int j=0;j<kmerLocations[kmerNumber].size();j++){//increase the count for each sequence that also has
+					matches[kmerLocations[kmerNumber][j]]++;	//	that kmer
+				}
 			}
+			timesKmerFound[kmerNumber] = 1;						//	ok, we've seen the kmer now
 		}
-		timesKmerFound[kmerNumber] = 1;						//	ok, we've seen the kmer now
-	}
-
-	for(int i=0;i<numSeqs;i++){								//	find the db sequence that shares the most kmers with
-		if(matches[i] > searchScore){					//	the query sequence
-			searchScore = matches[i];
-			maxSequence = i;
+		
+		if (num != 1) {
+			vector<seqMatch> seqMatches; seqMatches.resize(numSeqs);
+			for(int i=0;i<numSeqs;i++){		
+				seqMatches[i].seq = i;
+				seqMatches[i].match = matches[i];
+			}
+			
+			//sorts putting largest matches first
+			sort(seqMatches.begin(), seqMatches.end(), compareSeqMatches);
+			
+			searchScore = seqMatches[0].match;
+			searchScore = 100 * searchScore / (float) numKmers;		//	return the Sequence object corresponding to the db
+		
+			//save top matches
+			for (int i = 0; i < num; i++) {
+				topMatches.push_back(seqMatches[i].seq);
+				float thisScore = 100 * seqMatches[i].match / (float) numKmers;
+				Scores.push_back(thisScore);
+			}
+		}else{
+			int bestIndex = 0;
+			int bestMatch = -1;
+			for(int i=0;i<numSeqs;i++){	
+				
+				if (matches[i] > bestMatch) {
+					bestIndex = i;
+					bestMatch = matches[i];
+				}
+			}
+			
+			searchScore = bestMatch;
+			searchScore = 100 * searchScore / (float) numKmers;		//	return the Sequence object corresponding to the db
+			topMatches.push_back(bestIndex);
+			Scores.push_back(searchScore);
 		}
+		return topMatches;		
 	}
-
-	searchScore = 100 * searchScore / (float) numKmers;		//	return the Sequence object corresponding to the db
-	return templateSequences[maxSequence];					//	sequence with the most shared kmers.
+	catch(exception& e) {
+		m->errorOut(e, "KmerDB", "findClosestSequences");
+		exit(1);
+	}	
 }
 
 /**************************************************************************************************/
 
-void KmerDB::generateKmerDB(string kmerDBName){
-	
-	Kmer kmer(kmerSize);
+void KmerDB::generateDB(){
+	try {
+		
+		ofstream kmerFile;										//	once we have the kmerLocations folder print it out
+		m->openOutputFile(kmerDBName, kmerFile);					//	to a file
+		
+		//output version
+		kmerFile << "#" << m->getVersion() << endl;
+		
+		for(int i=0;i<maxKmer;i++){								//	step through all of the possible kmer numbers
+			kmerFile << i << ' ' << kmerLocations[i].size();	//	print the kmer number and the number of sequences with
+			for(int j=0;j<kmerLocations[i].size();j++){			//	that kmer.  then print out the indices of the sequences
+				kmerFile << ' ' << kmerLocations[i][j];			//	with that kmer.
+			}
+			kmerFile << endl;
+		}
+		kmerFile.close();
+		
+	}
+	catch(exception& e) {
+		m->errorOut(e, "KmerDB", "generateDB");
+		exit(1);
+	}	
 	
-	for(int i=0;i<numSeqs;i++){								//	for all of the template sequences...
-
-		string seq = templateSequences[i].getUnaligned();	//	...take the unaligned sequence...
-		int numKmers = seq.length() - kmerSize + 1;
+}
+/**************************************************************************************************/
+void KmerDB::addSequence(Sequence seq) {
+	try {
+		Kmer kmer(kmerSize);
 		
+		string unaligned = seq.getUnaligned();	//	...take the unaligned sequence...
+		int numKmers = unaligned.length() - kmerSize + 1;
+			
 		vector<int> seenBefore(maxKmer+1,0);
 		for(int j=0;j<numKmers;j++){						//	...step though the sequence and get each kmer...
-			int kmerNumber = kmer.getKmerNumber(seq, j);
+			int kmerNumber = kmer.getKmerNumber(unaligned, j);
 			if(seenBefore[kmerNumber] == 0){
-				kmerLocations[kmerNumber].push_back(i);		//	...insert the sequence index into kmerLocations for
+				kmerLocations[kmerNumber].push_back(count);		//	...insert the sequence index into kmerLocations for
 			}												//	the appropriate kmer number
 			seenBefore[kmerNumber] = 1;
 		}													
-	}
-	
-	ofstream kmerFile;										//	once we have the kmerLocations folder print it out
-	openOutputFile(kmerDBName, kmerFile);					//	to a file
 	
-	for(int i=0;i<maxKmer;i++){								//	step through all of the possible kmer numbers
-		kmerFile << i << ' ' << kmerLocations[i].size();	//	print the kmer number and the number of sequences with
-		for(int j=0;j<kmerLocations[i].size();j++){			//	that kmer.  then print out the indices of the sequences
-			kmerFile << ' ' << kmerLocations[i][j];			//	with that kmer.
-		}
-		kmerFile << endl;
+		count++;
 	}
-	kmerFile.close();
-	
+	catch(exception& e) {
+		m->errorOut(e, "KmerDB", "addSequence");
+		exit(1);
+	}	
 }
-
 /**************************************************************************************************/
 
-void KmerDB::readKmerDB(string kmerDBName, ifstream& kmerDBFile){
+void KmerDB::readKmerDB(ifstream& kmerDBFile){
+	try {
+					
+		kmerDBFile.seekg(0);									//	start at the beginning of the file
+		
+		//read version
+		string line = m->getline(kmerDBFile); m->gobble(kmerDBFile);
+		
+		string seqName;
+		int seqNumber;
+
+		for(int i=0;i<maxKmer;i++){
+			int numValues = 0;	
+			kmerDBFile >> seqName >> numValues;
+			
+			for(int j=0;j<numValues;j++){						//	for each kmer number get the...
+				kmerDBFile >> seqNumber;						//		1. number of sequences with the kmer number
+				kmerLocations[i].push_back(seqNumber);			//		2. sequence indices
+			}
+		}
+		kmerDBFile.close();
+		
+	}
+	catch(exception& e) {
+		m->errorOut(e, "KmerDB", "readKmerDB");
+		exit(1);
+	}	
+}
 
-	kmerDBFile.seekg(0);									//	start at the beginning of the file
-	
-	string seqName;
-	int seqNumber;
+/**************************************************************************************************/
+int KmerDB::getCount(int kmer) {
+	try {
+		if (kmer < 0) { return 0; }  //if user gives negative number
+		else if (kmer > maxKmer) {	return 0;	}  //or a kmer that is bigger than maxkmer
+		else {	return kmerLocations[kmer].size();	}  // kmer is in vector range
+	}
+	catch(exception& e) {
+		m->errorOut(e, "KmerDB", "getCount");
+		exit(1);
+	}	
+}
+/**************************************************************************************************/
+int KmerDB::getReversed(int kmerNumber) {
+	try {
+		Kmer kmer(kmerSize);
+		
+		if (kmerNumber < 0) { return 0; }  //if user gives negative number
+		else if (kmerNumber > maxKmer) {	return 0;	}  //or a kmer that is bigger than maxkmer
+		else {	return kmer.getReverseKmerNumber(kmerNumber);	}  // kmer is in vector range
+	}
+	catch(exception& e) {
+		m->errorOut(e, "KmerDB", "getReversed");
+		exit(1);
+	}	
+}
+/**************************************************************************************************/
+vector<int> KmerDB::getSequencesWithKmer(int kmer) {
+	try {
+		
+		vector<int> seqs;
 	
-	for(int i=0;i<maxKmer;i++){
-		int numValues;	
-		kmerDBFile >> seqName >> numValues;
+		if (kmer < 0) { }  //if user gives negative number
+		else if (kmer > maxKmer) {	}  //or a kmer that is bigger than maxkmer
+		else {	seqs = kmerLocations[kmer];	}
 		
-		for(int j=0;j<numValues;j++){						//	for each kmer number get the...
-			kmerDBFile >> seqNumber;						//		1. number of sequences with the kmer number
-			kmerLocations[i].push_back(seqNumber);			//		2. sequence indices
-		}
+		return seqs;
 	}
-	kmerDBFile.close();
+	catch(exception& e) {
+		m->errorOut(e, "KmerDB", "getSequencesWithKmer");
+		exit(1);
+	}	
 }
+/**************************************************************************************************/
+
 
 /**************************************************************************************************/
+
+