RabbitTClust v.2.2.0, support incremental clustering by --append option

Author: XiaomingXu1995

README.md

@@ -1,12 +1,12 @@
 ![RabbitTClust](rabbittclust.png)
 
-# `RabbitTClust v.2.1.0`
+# `RabbitTClust v.2.2.0`
 RabbitTClust is a fast and memory-efficient genome clustering tool based on sketch-based distance estimations.
 It enables processing of large-scale datasets by combining dimensionality reduction techniques with streaming and parallelization on modern multi-core platforms.
 RabbitTClust supports classical single-linkage hierarchical (clust-mst) and greedy incremental clustering (clust-greedy) algorithms for different scenarios. 
 
 ## Installation
-RabbitTClust version 2.1.0 can only support 64-bit Linux Systems.
+`RabbitTClust v.2.2.0` can only support 64-bit Linux Systems.
 
 The detailed update information for this version, as well as the version history, can be found in the [`version_history`](version_history/history.md) document.
 
@@ -15,45 +15,12 @@ The detailed update information for this version, as well as the version history
 * c++14
 * [zlib](https://zlib.net/)
 
-### Compile and install automatically
+### Compile and install
 ```bash
 git clone --recursive https://github.com/RabbitBio/RabbitTClust.git
 cd RabbitTClust
 ./install.sh
 ```
-
-### Compile and install manually 
-```bash
-git clone --recursive https://github.com/RabbitBio/RabbitTClust.git
-cd RabbitTClust
-
-#make rabbitSketch library
-cd RabbitSketch &&
-mkdir -p build && cd build &&
-cmake -DCXXAPI=ON -DCMAKE_INSTALL_PREFIX=. .. &&
-make -j8 && make install &&
-cd ../../ &&
-
-#make rabbitFX library
-cd RabbitFX && 
-mkdir -p build && cd build &&
-cmake -DCMAKE_INSTALL_PREFIX=. .. &&
-make -j8 && make install && 
-cd ../../ &&
-
-#compile the clust-greedy
-mkdir -p build && cd build &&
-cmake -DUSE_RABBITFX=ON -DUSE_GREEDY=ON .. && 
-make -j8 && make install &&
-cd ../ &&
-
-#compile the clust-mst
-cd build &&
-cmake -DUSE_RABBITFX=ON -DUSE_GREEDY=OFF .. &&
-make -j8 && make install &&
-cd ../ 
-```
-
 ## Usage
 ```bash
 # clust-mst, minimum-spanning-tree-based module for RabbitTClust
@@ -62,65 +29,75 @@ Options:
   -h,--help                   Print this help message and exit
   -t,--threads INT            set the thread number, default all CPUs of the platform
   -m,--min-length UINT        set the filter minimum length (minLen), genome length less than minLen will be ignore, default 10,000
-  -c,--containment INT        use AAF distance with containment coefficient, set the containCompress, the sketch size is in proportion with 1/containCompress
-  -k,--kmer-size INT          set the kmer size
+  -c,--containment INT        use AAF distance with containment coefficient, set the containCompress, the sketch size is in proportion with 1/containCompress  -k,--kmer-size INT          set the kmer size
   -s,--sketch-size INT        set the sketch size for Jaccard Index and Mash distance, default 1000
-  -l,--inputlist              input is genome list, one genome per line
+  -l,--list                   input is genome list, one genome per line
   -e,--no-save                not save the intermediate files, such as sketches or MST
   -d,--threshold FLOAT        set the distance threshold for clustering
   -F,--function TEXT          set the sketch function, such as MinHash, KSSD, default MinHash
   -o,--output TEXT REQUIRED   set the output name of cluster result
-  -i,--input TEXT             set the input file
+  -i,--input TEXT Excludes: --append
+                              set the input file, single FASTA genome file (without -l option) or genome list file (with -l option)
   --presketched TEXT          clustering by the pre-generated sketch files rather than genomes
   --premsted TEXT             clustering by the pre-generated mst files rather than genomes for clust-mst
+  --append TEXT Excludes: --input
+                              append genome file or file list with the pre-generated sketch or MST files
 
 # clust-greedy, greedy incremental clustering module for RabbitTClust
 Usage: ./clust-greedy [OPTIONS]
 Options:
   -h,--help                   Print this help message and exit
   -t,--threads INT            set the thread number, default all CPUs of the platform
   -m,--min-length UINT        set the filter minimum length (minLen), genome length less than minLen will be ignore, default 10,000
-  -c,--containment INT        use AAF distance with containment coefficient, set the containCompress, the sketch size is in proportion with 1/containCompress
-  -k,--kmer-size INT          set the kmer size
+  -c,--containment INT        use AAF distance with containment coefficient, set the containCompress, the sketch size is in proportion with 1/containCompress  -k,--kmer-size INT          set the kmer size
   -s,--sketch-size INT        set the sketch size for Jaccard Index and Mash distance, default 1000
-  -l,--inputlist              input is genome list, one genome per line
+  -l,--list                   input is genome list, one genome per line
   -e,--no-save                not save the intermediate files, such as sketches or MST
   -d,--threshold FLOAT        set the distance threshold for clustering
   -F,--function TEXT          set the sketch function, such as MinHash, KSSD, default MinHash
   -o,--output TEXT REQUIRED   set the output name of cluster result
-  -i,--input TEXT             set the input file
+  -i,--input TEXT Excludes: --append
+                              set the input file, single FASTA genome file (without -l option) or genome list file (with -l option)
   --presketched TEXT          clustering by the pre-generated sketch files rather than genomes
+  --append TEXT Excludes: --input
+                              append genome file or file list with the pre-generated sketch or MST files
 ```
 
 ## Example:
 ```bash
-#input is a file list, one genome path per line:
+# input is a file list, one genome path per line:
 ./clust-mst -l -i bact_refseq.list -o bact_refseq.mst.clust
 ./clust-greedy -l -i bact_genbank.list -o bact_genbank.greedy.clust
 
-#input is a single genome file in FASTA format, one genome as a sequence:
+# input is a single genome file in FASTA format, one genome as a sequence:
 ./clust-mst -i bacteria.fna -o bacteria.mst.clust
 ./clust-greedy -i bacteria.fna -o bacteria.greedy.clust
 
-#the sketch size (reciprocal of sampling proportion), kmer size, and distance threshold can be specified by -s (-c), -k, and -d options.
+# the sketch size (reciprocal of sampling proportion), kmer size, and distance threshold can be specified by -s (-c), -k, and -d options.
 ./clust-mst -l -k 21 -s 1000 -d 0.05 -i bact_refseq.list -o bact_refseq.mst.clust
 ./clust-greedy -l -k 21 -c 1000 -d 0.05 -i bact_genbank.list -o bact_genbank.greedy.clust
 
 
-#for redundancy detection with clust-greedy, input is a genome file list:
-#use -d to specify the distance threshold corresponding to various degrees of redundancy.
+# for redundancy detection with clust-greedy, input is a genome file list:
+# use -d to specify the distance threshold corresponding to various degrees of redundancy.
 ./clust-greedy -d 0.001 -l -i bacteria.list -o bacteria.out
 
-#for last running of clust-mst, it generated a folder name in year_month_day_hour-minute-second format, such as 2023_05_06_08-49-15.
-#this folder contains the sketch, mst files.
-#for generator cluster from exist MST with a distance threshold of 0.045:
-./clust-mst -d 0.045 --premsted 2023_05_06_08-49-15 -o bacteria.mst.d.045.clust
-#for generator cluster from exist sketches files of clust-mst with a distance threshold of 0.045:
-./clust-mst -d 0.045 --presketched 2023_05_06_08-49-15 -o bacteria.mst.d.045.clust
+# v.2.1.0 or later
+# for last running of clust-mst, it generated a folder name in year_month_day_hour-minute-second format, such as 2023_05_06_08-49-15.
+# this folder contains the sketch, mst files.
+# for generator cluster from exist MST with a distance threshold of 0.045:
+./clust-mst -d 0.045 --premsted 2023_05_06_08-49-15/ -o bact_refseq.mst.d.045.clust
+# for generator cluster from exist sketches files of clust-mst with a distance threshold of 0.045:
+./clust-mst -d 0.045 --presketched 2023_05_06_08-49-15/ -o bact_refseq.mst.d.045.clust
 
-#for generator cluster from exist sketches of clust-greedy with a distance threshold of 0.001:
+# for generator cluster from exist sketches of clust-greedy with a distance threshold of 0.001:
 # folder 2023_05_06_08-49-15 contains the sketch files.
-./clust-greedy -d 0.001 --presketched 2023_05_06_08-49-15 -o bact_genbank.greedy.d.001.clust
+./clust-greedy -d 0.001 --presketched 2023_05_06_09-37-23/ -o bact_genbank.greedy.d.001.clust
+
+# v.2.2.0 or later
+# for generator cluster from exist part sketches (presketch_A_dir) and append genome set (genome_B.list) to incrementally clustering 
+./clust-mst --presketched 2023_05_06_08-49-15/ -l --append genome_B.list -o append_refseq.mst.clust
+./clust-mst --presketched 2023_05_06_09-37-23/ -l --append genome_B.list -o append_genbank.greedy.clust
 ```
 ## Output
 The output file is in a CD-HIT output format and is slightly different when running with varying input options (*-l* and *-i*).  

src/MST.cpp

@@ -171,11 +171,12 @@ vector<int> getNoiseNode(vector<PairInt> densePairArr, int alpha){
 	return noiseArr;
 }
 
-vector<EdgeInfo> modifyMST(vector<SketchInfo>& sketches, int sketch_func_id, int threads, int** &denseArr, int denseSpan, uint64_t* &aniArr, string prefixName, double threshold){
+
+
+vector<EdgeInfo> modifyMST(vector<SketchInfo>& sketches, int start_index, int sketch_func_id, int threads, int** &denseArr, int denseSpan, uint64_t* &aniArr){
 	//int denseSpan = 10;
 	double step = 1.0 / denseSpan;
 
-	
 	//double step = threshold / denseSpan;
 	//cerr << "the threshold is: " << threshold << endl;
 	//cerr << "the step is: " << step << endl;
@@ -211,13 +212,14 @@ vector<EdgeInfo> modifyMST(vector<SketchInfo>& sketches, int sketch_func_id, int
 	//int N = sketches.size();
 	uint64_t totalCompNum = (uint64_t)N * (uint64_t)(N-1)/2;
 	uint64_t percentNum = totalCompNum / 100;
-	//cerr << "the percentNum is: " << percentNum << endl;
+	cerr << "---the percentNum is: " << percentNum << endl;
+	cerr << "---the start_index is: " << start_index << endl;
 	uint64_t percentId = 0;
 	#pragma omp parallel for num_threads(threads) schedule (dynamic)
 	for(id = 0; id < sketches.size() - tailNum; id+=subSize){
 		int thread_id = omp_get_thread_num();
 		for(int i = id; i < id+subSize; i++){
-			for(int j = i+1; j < sketches.size(); j++){
+			for(int j = max(i+1, start_index); j < sketches.size(); j++){
 				double tmpDist;
 				if(sketch_func_id == 0)
 				{

src/MST.h

@@ -44,7 +44,9 @@ std::vector<EdgeInfo> kruskalAlgorithm(std::vector<EdgeInfo>graph, int vertices)
 
 vector<EdgeInfo> generateMST(vector<SketchInfo>& sketches, string sketchFunc, int threads);
 
-vector<EdgeInfo> modifyMST(vector<SketchInfo>& sketches, int sketch_func_id, int threads, int** &denseArr, int denseSpan, uint64_t* &aniArr, string prefixName, double threshold);
+vector<EdgeInfo> append_MST(vector<SketchInfo>& pre_sketches, vector<SketchInfo>& append_sketches, int sketch_func_id, int threads, int ** &denseArr, int denseSpan, uint64_t* &aniArr);
+
+vector<EdgeInfo> modifyMST(vector<SketchInfo>& sketches, int start_index, int sketch_func_id, int threads, int** &denseArr, int denseSpan, uint64_t* &aniArr);
 
 std::vector<EdgeInfo> generateForest(std::vector <EdgeInfo> mst, double threshhold);
 

src/MST_IO.cpp

@@ -26,15 +26,30 @@ void loadDense(int** &denseArr, string folderPath, int& denseSpan, int& genome_n
 	cerr << "-----read the dense file from: " << file_dense << endl;
 }
 
+void loadANI(string folderPath, uint64_t* &aniArr, int sketch_func_id){
+	if(sketch_func_id != 0 && sketch_func_id != 1){
+		cerr << "ERROR: saveANI(), save ANI can only support MinHash and KSSD functions" << endl;
+		return;
+	}
+	string file_ani = folderPath + '/' + "mst.ani";
+	FILE* fp_ani = fopen(file_ani.c_str(), "r");
+	if(!fp_ani){
+		cerr << "ERROR: saveANI(), cannot open file: " << file_ani << endl;
+		exit(1);
+	}
+	aniArr = new uint64_t[101];
+	fread(aniArr, sizeof(uint64_t), 101, fp_ani);
+	fclose(fp_ani);
+	cerr << "-----read the ani file from: " << file_ani << endl;
+}
 
-bool loadMST(string folderPath, vector<SketchInfo>& sketches, vector<EdgeInfo>& mst)
+void loadMST(string folderPath, vector<EdgeInfo>& mst)
 {
-	bool sketch_by_file = load_genome_info(folderPath, "mst", sketches);
 	//load the mst edge 
 	string file_mst = folderPath + '/' + "edge.mst";
 	FILE* fp_mst = fopen(file_mst.c_str(), "r");
 	if(!fp_mst){
-		cerr << "ERROR: saveMST(), cannot open the file: " <<  file_mst << endl;
+		cerr << "ERROR: loadMST(), cannot open the file: " <<  file_mst << endl;
 		exit(1);
 	}
 	size_t mst_size;
@@ -51,7 +66,6 @@ bool loadMST(string folderPath, vector<SketchInfo>& sketches, vector<EdgeInfo>&
 	}
 	fclose(fp_mst);
 	cerr << "-----read the mst file from " << file_mst << endl;
-	return sketch_by_file;
 }
 
 void printResult(vector<vector<int>>& cluster, vector<SketchInfo>& sketches, bool sketchByFile, string outputFile)
@@ -88,7 +102,7 @@ void printResult(vector<vector<int>>& cluster, vector<SketchInfo>& sketches, boo
 
 }
 
-void saveMST(vector<SketchInfo> sketches, vector<EdgeInfo> mst, string folderPath, bool sketchByFile){
+void saveMST(vector<SketchInfo>& sketches, vector<EdgeInfo>& mst, string folderPath, bool sketchByFile){
 	save_genome_info(sketches, folderPath, "mst", sketchByFile);
 	string file_mst = folderPath + '/' + "edge.mst";
 	FILE* fp_mst = fopen(file_mst.c_str(), "w+");
@@ -137,7 +151,7 @@ void saveANI(string folderPath, uint64_t* aniArr, int sketch_func_id){
 	}
 	fwrite(aniArr, sizeof(uint64_t), 101, fp_ani);
 	fclose(fp_ani);
-	cerr << "-----save the ani file into: " << folderPath << endl;
+	cerr << "-----save the ani file into: " << file_ani << endl;
 }
 
 

src/MST_IO.h

@@ -11,22 +11,14 @@ struct ClusterInfo{
 	int id;
 	uint64_t length;
 };
-//void loadDense(int** &denseArr, string inputFile, int denseSpan, vector<SketchInfo> sketches);
-
-//bool loadMSTs(string inputFile, string inputFile1, vector<SketchInfo>& sketches, vector<EdgeInfo>& mst);
 
 void printResult(std::vector<std::vector<int>>& clusterOrigin, std::vector<SketchInfo>& sketches, bool sketchByFile, string outputFile);
 
-bool loadMST(string folderPath, vector<SketchInfo>& sketches, vector<EdgeInfo>& mst);
+void loadMST(string folderPath, vector<EdgeInfo>& mst);
 void loadDense(int** &denseArr, string folderPath, int& denseSpan, int& genome_number);
-//void saveMST(string folderPath, string inputFile, string sketchFunc, bool isContainment, int containCompress, vector<SketchInfo> sketches, vector<EdgeInfo> mst, bool sketchByFile, int sketchSize, int kmerSize);
-//
-//void saveDense(string folderPath, string prefixName, int** denseArr, int denseSpan, vector<SketchInfo> sketches);
-//
-//void saveANI(string folderPath, string prefixName, uint64_t* aniArr, string sketchFunc);
-
+void loadANI(string folderPath, uint64_t* &aniArr, int sketch_func_id);
 
-void saveMST(vector<SketchInfo> sketches, vector<EdgeInfo> mst, string folderPath, bool sketchByFile);
+void saveMST(vector<SketchInfo>& sketches, vector<EdgeInfo>& mst, string folderPath, bool sketchByFile);
 void saveDense(string folderPath, int** denseArr, int denseSpan, int genome_number);
 void saveANI(string folderPath, uint64_t* aniArr, int sketch_func_id);
 

src/SketchInfo.h

@@ -35,6 +35,8 @@ struct SketchInfo{
 };
 
 
+bool cmpGenomeSize(SketchInfo s1, SketchInfo s2);
+bool cmpSeqSize(SketchInfo s1, SketchInfo s2);
 
 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, int minLen, string sketchFunc, bool isContainment, int containCompress, vector<SketchInfo>& sketches, int threads);

src/Sketch_IO.cpp

@@ -9,12 +9,36 @@ bool cmpIndex(SketchInfo s1, SketchInfo s2){
 	return s1.id < s2.id;
 }
 
+void read_sketch_parameters(string folder_path, int& sketch_func_id, int& kmer_size, bool& is_containment, int& contain_compress, int& sketch_size, int& half_k, int& half_subk, int& drlevel){
+	string hash_file = folder_path + '/' + "hash.sketch";
+	FILE * fp_hash = fopen(hash_file.c_str(), "r");
+	if(!fp_hash){
+		cerr << "ERROR: read_sketch_parameters(), cannot open the file: " << hash_file << endl;
+		exit(1);
+	}
+	fread(&sketch_func_id, sizeof(int), 1, fp_hash);
+	if(sketch_func_id == 0){
+		fread(&kmer_size, sizeof(int), 1, fp_hash);
+		fread(&is_containment, sizeof(bool), 1, fp_hash);
+		if(is_containment)
+			fread(&contain_compress, sizeof(int), 1, fp_hash);
+		else
+			fread(&sketch_size, sizeof(int), 1, fp_hash);
+	}
+	else if(sketch_func_id == 1){
+		fread(&half_k, sizeof(int), 1, fp_hash);
+		fread(&half_subk, sizeof(int), 1, fp_hash);
+		fread(&drlevel, sizeof(int), 1, fp_hash);
+	}
+	fclose(fp_hash); 
+}
+
 void save_genome_info(vector<SketchInfo>& sketches, string folderPath, string type, bool sketchByFile){
 	assert(type == "sketch" || type == "mst");
 	string info_file = folderPath + '/' + "info." + type;
 	FILE * fp_info = fopen(info_file.c_str(), "w+");
 	if(!fp_info){
-		cerr << "ERROR: saveSketches(), cannot open the file: " << info_file << endl;
+		cerr << "ERROR: save_genome_info(), cannot open the file: " << info_file << endl;
 		exit(1);
 	}
 	fwrite(&sketchByFile, sizeof(bool), 1, fp_info);

src/Sketch_IO.h

@@ -3,6 +3,7 @@
 #include "SketchInfo.h"
 #include "common.hpp"
 
+void read_sketch_parameters(string folder_path, int& sketch_func_id, int& kmer_size, bool& is_containment, int& contain_compress, int& sketch_size, int& half_k, int& half_subk, int& drlevel);
 void save_genome_info(vector<SketchInfo>& sketches, string folderPath, string type, bool sketchByFile);
 void saveSketches(vector<SketchInfo>& sketches, string folderPath, bool sketchByFile, string sketchFunc, bool isContainment, int containCompress, int sketchSize, int kmerSize);
 

src/main.cpp

@@ -95,8 +95,10 @@ int main(int argc, char * argv[]){
 #ifndef GREEDY_CLUST
 	auto option_premsted = app.add_option("--premsted", folder_path, "clustering by the pre-generated mst files rather than genomes for clust-mst");
 #endif
+	auto option_append = app.add_option("--append", inputFile, "append genome file or file list with the pre-generated sketch or MST files");
 
 	option_output->required();
+	option_append->excludes(option_input);
 
 	CLI11_PARSE(app, argc, argv);
 
@@ -132,22 +134,41 @@ int main(int argc, char * argv[]){
 
 #ifndef GREEDY_CLUST
 //======clust-mst=========================================================================
-	if(*option_premsted){
+	if(*option_premsted && !*option_append){
 		clust_from_mst(folder_path, outputFile, threshold, threads);
 		return 0;
 	}
+	if(*option_append && !*option_presketched && !*option_premsted){
+		cerr << "ERROR option --append, option --presketched or --premsted needed" << endl;
+		return 1;
+	}
+	if(*option_append && (*option_premsted || *option_presketched)){
+		append_clust_mst(folder_path, inputFile, outputFile, sketchByFile, minLen, noSave, threshold, threads);
+		return 0;
+	}
 //======clust-mst=========================================================================
+#else
+//======clust-greedy======================================================================
+	if(*option_append && !*option_presketched){
+		cerr << "ERROR option --append, option --presketched needed" << endl;
+		return 1;
+	}
+	if(*option_append && *option_presketched){
+		append_clust_greedy(folder_path, inputFile, outputFile, sketchByFile, minLen, noSave, threshold, threads);
+		return 0;
+	}
+//======clust-greedy======================================================================
 #endif
 
-	if(*option_presketched){
+	if(*option_presketched && !*option_append){
 		clust_from_sketches(folder_path, outputFile, threshold, threads);
 		return 0;
 	}
 
 	if(!tune_parameters(sketchByFile, isSetKmer, inputFile, threads, minLen, isContainment, isJaccard, kmerSize, threshold, containCompress, sketchSize)){
 		return 1;
 	}
-
+	
 	clust_from_genomes(inputFile, outputFile, sketchByFile, kmerSize, sketchSize, threshold,sketchFunc, isContainment, containCompress, minLen, folder_path, noSave, threads);
 
 	return 0;