add the parameter -m to set the minimum genome length to filter out short genome

Author: XiaomingXu1995

src/SketchInfo.cpp

@@ -112,15 +112,15 @@ int producer_fasta_task(std::string file, rabbit::fa::FastaDataPool* fastaPool,
 	return 0;
 }
 
-void consumer_fasta_seqSize(rabbit::fa::FastaDataPool* fastaPool, FaChunkQueue &dq, uint64_t* maxSize, uint64_t* minSize, uint64_t* totalSize,	uint64_t* number, uint64_t* badNumber){
+void consumer_fasta_seqSize(rabbit::fa::FastaDataPool* fastaPool, FaChunkQueue &dq, uint64_t minLen, uint64_t* maxSize, uint64_t* minSize, uint64_t* totalSize,	uint64_t* number, uint64_t* badNumber){
 	rabbit::int64 id = 0;
 	rabbit::fa::FastaChunk *faChunk;
 	while(dq.Pop(id, faChunk)){
 		std::vector<Reference> data;
 		int ref_num = rabbit::fa::chunkListFormat(*faChunk, data);
 		for(Reference &r: data){
 			uint64_t length = (uint64_t)r.length;
-			if(length < 10000){
+			if(length < minLen){
 				*badNumber = *badNumber + 1;
 				continue;
 			}
@@ -209,7 +209,7 @@ void consumer_fasta_task(rabbit::fa::FastaDataPool* fastaPool, FaChunkQueue &dq,
 }
 #endif
 
-void calSize(bool sketchByFile, string inputFile, int threads, uint64_t &maxSize, uint64_t& minSize, uint64_t& averageSize){
+void calSize(bool sketchByFile, string inputFile, int threads, uint64_t minLen, uint64_t &maxSize, uint64_t& minSize, uint64_t& averageSize){
 	maxSize = 0;
 	minSize = 1 << 31;
 	averageSize = 0;
@@ -238,7 +238,7 @@ void calSize(bool sketchByFile, string inputFile, int threads, uint64_t &maxSize
 				stat(line.c_str(), &statbuf);
 				curSize = statbuf.st_size;
 			}
-			if(curSize < 10000){
+			if(curSize < minLen){
 				badNumber++;
 				continue;
 			}
@@ -270,7 +270,7 @@ void calSize(bool sketchByFile, string inputFile, int threads, uint64_t &maxSize
 		std::thread producer(producer_fasta_task, inputFile, fastaPool, std::ref(queue1));
 		std::thread **threadArr = new std::thread* [th];
 		for(int t = 0; t < th; t++){
-			threadArr[t] = new std::thread(std::bind(consumer_fasta_seqSize, fastaPool, std::ref(queue1), &maxSizeArr[t], &minSizeArr[t], &totalSizeArr[t], &numArr[t], &badNumArr[t]));
+			threadArr[t] = new std::thread(std::bind(consumer_fasta_seqSize, fastaPool, std::ref(queue1), minLen, &maxSizeArr[t], &minSizeArr[t], &totalSizeArr[t], &numArr[t], &badNumArr[t]));
 		}
 		producer.join();
 
@@ -294,7 +294,7 @@ void calSize(bool sketchByFile, string inputFile, int threads, uint64_t &maxSize
 			while(1){
 				int length = kseq_read(ks1);
 				if(length < 0) break;
-				if(length < 10000){
+				if(length < minLen){
 					badNumber++;
 					continue;
 				}
@@ -313,7 +313,7 @@ void calSize(bool sketchByFile, string inputFile, int threads, uint64_t &maxSize
 
 	if((double)badNumber / totalNumber >= 0.2)
 	{
-		fprintf(stderr, "Warning: there are %d poor quality (length < 10000) genome assemblies in the total %d genome assemblied.\n", badNumber, totalNumber);
+		fprintf(stderr, "Warning: there are %d poor quality (length < %ld) genome assemblies in the total %d genome assemblied.\n", badNumber, minLen, totalNumber);
 	}
 	cerr << "\t===the totalSize is: " << totalSize << endl;
 	cerr << "\t===the maxSize is: " << maxSize << endl;
@@ -520,7 +520,7 @@ bool sketchSequences(string inputFile, int kmerSize, int sketchSize, string sket
 	return true;
 }
 
-bool sketchFiles(string inputFile, int kmerSize, int sketchSize, string sketchFunc, bool isContainment, int containCompress, vector<SketchInfo>& sketches, int threads){
+bool sketchFiles(string inputFile, uint64_t minLen, int kmerSize, int sketchSize, string sketchFunc, bool isContainment, int containCompress, vector<SketchInfo>& sketches, int threads){
 	fprintf(stderr, "input fileList, sketch by file\n");
 	fstream fs(inputFile);
 	if(!fs){
@@ -691,7 +691,7 @@ bool sketchFiles(string inputFile, int kmerSize, int sketchSize, string sketchFu
 			tmpSketchInfo.fileName = fileList[i];
 			tmpSketchInfo.totalSeqLength = totalLength;
 			tmpSketchInfo.fileSeqs = curFileSeqs;
-			if(totalLength >= 10000)//filter the poor quality genome assemblies whose length less than 10k bp(fastANI paper)
+			if(totalLength >= minLen)//filter the poor quality genome assemblies whose length less than minLen(fastANI paper)
 				sketches.push_back(tmpSketchInfo);
 			if(i % 10000 == 0)	cerr << "finished sketching: " << i << " genomes" << endl;
 		}

src/SketchInfo.h

@@ -36,9 +36,9 @@ struct SketchInfo{
 
 
 
-void calSize(bool sketchByFile, string inputFile, int threads, uint64_t &maxSize, uint64_t& minSize, uint64_t& averageSize);
+void calSize(bool sketchByFile, string inputFile, int threads, uint64_t minLen, uint64_t &maxSize, uint64_t& minSize, uint64_t& averageSize);
 bool sketchSequences(string inputFile, int kmerSize, int sketchSize, string sketchFunc, bool isContainment, int containCompress, vector<SketchInfo>& sketches, int threads);
-bool sketchFiles(string inputFile, int kmerSize, int sketchSize, string sketchFunc, bool isContainment, int containCompress, vector<SketchInfo>& sketches, int threads);
+bool sketchFiles(string inputFile, uint64_t minLen, int kmerSize, int sketchSize, string sketchFunc, bool isContainment, int containCompress, vector<SketchInfo>& sketches, int threads);
 bool cmpSketch(SketchInfo s1, SketchInfo s2);
 
 

src/main.cpp

@@ -63,6 +63,7 @@ int main(int argc, char * argv[]){
 	string mstSketchFile = "sketch.info";
 	bool isSave = true;
 	bool isSetKmer = false;
+	uint64_t minLen = 0;
 	while(argIndex < argc){
 		switch(argv[argIndex][1]){
 			case 't':
@@ -84,6 +85,10 @@ int main(int argc, char * argv[]){
 			case 'e':
 				isSave = false;
 				break;
+			case 'm':
+				minLen = stoi(argv[++argIndex]);
+				fprintf(stderr, "set the filter minimum length: %ld\n", minLen);
+				break;
 			case 'k':
 				isSetKmer = true;
 				kmerSize = stoi(argv[++argIndex]);
@@ -237,7 +242,7 @@ int main(int argc, char * argv[]){
 	}//end useFile
 
 	uint64_t maxSize, minSize, averageSize;
-	calSize(sketchByFile, inputFile, threads, maxSize, minSize, averageSize);
+	calSize(sketchByFile, inputFile, threads, minLen, maxSize, minSize, averageSize);
 
 	if(isContainment && isJaccard){
 		cerr << "conflict distance measurement of Mash distance (fixed-sketch-size) and AAF distance (variable-sketch-size) " << endl;
@@ -334,7 +339,7 @@ int main(int argc, char * argv[]){
 
 		}//end sketch by sequence
 		else{
-			if(!sketchFiles(inputFile, kmerSize, sketchSize, sketchFunc, isContainment, containCompress, sketches, threads)){
+			if(!sketchFiles(inputFile, minLen, kmerSize, sketchSize, sketchFunc, isContainment, containCompress, sketches, threads)){
 				printUsage();
 				return 1;
 			}

src/parameter.h

@@ -47,6 +47,7 @@ inline void printUsage(){
 	fprintf(stdout, "usage: clust-greedy [-h] [-l] [-t] <int> [-d] <double> [-F] <string> [-i] <string> [-o] <string> \n");
 	fprintf(stdout, "usage: clust-greedy [-h] [-f] [-d] <double> [-i] <string> <string> [-o] <string>\n");
 	fprintf(stdout, "	-h\t\t: this help message\n");
+	fprintf(stdout, "	-m <int>\t: set the filter minimum genome length (minLen), genome with total length less the minLen will be ignore, for both clust-mst and clust-greedy\n");
 	fprintf(stdout, "	-k <int>\t: set kmer size, automatically calculate the kmer size without -k option, for both clust-mst and clust-greedy\n");
 	fprintf(stdout, "	-s <int>\t: set sketch size, default 1000, for both clust-mst and clust-greedy\n");
 	fprintf(stdout, "	-c <int>\t: set sampling ratio to compute viriable sketchSize, sketchSize = genomeSize/samplingRatio, only support with MinHash sketch function of clust-greedy\n");