Merge pull request #145 from jaebeom-kim/master

jaebeom-kim · web-flow · commit 56836b6d28bb · 2025-05-29T10:13:36.000+09:00
- Fix a bug of printing header line multiple times
- Improve command grouping
- add --validate-db option
- use sse4.1 and ips4o for windows os.
diff --git a/README.md b/README.md
@@ -163,6 +163,7 @@ metabuli classify --seq-mode 3 read.fna dbdir outdir jobid
 
   * Important parameters:
    --validate-input : Validate query file format (0 by default)
+   --validate-db : Validate DB files (0 by default)
    --threads : The number of threads used (all by default)
    --max-ram : The maximum RAM usage. (128 GiB by default)
    --min-score : The minimum score to be classified 
@@ -313,6 +314,7 @@ metabuli build --gtdb 1 <DBDIR> <FASTA_LIST> <GTDB_TAXDUMP/taxid.map> --taxonomy
    --accession-level : Set 1 to creat a DB for accession level classification (0 by default).
    --cds-info : List of absolute paths to CDS files.
    --validate-input : Validate FASTA file format (0 by default)
+   --validate-db : Validate created DB files (0 by default)
   
 ```
 This will generate **diffIdx**, **info**, **split**, and **taxID_list** and some other files. You can delete `*_diffIdx` and `*_info` files.
@@ -340,6 +342,7 @@ metabuli updateDB --gtdb 1 <NEW DBDIR> <FASTA_LIST> <GTDB_TAXDUMP/taxid.map> <OL
   --accession-level: Set 1 to add new sequences for accession level classification (0 by default).
   --cds-info: List of absolute paths to CDS files.
   --validate-input : Validate FASTA file format (0 by default)
+  --validate-db : Validate created DB files (0 by default)
 ```
 
 #### \<Add sequences of new taxa>
@@ -429,6 +432,7 @@ metabuli build <DBDIR> <FASTA_LIST> <accession2taxid> --taxonomy-path <TAXDUMP>
   --accession-level: Set 1 to creat a DB for accession level classification (0 by default).
   --cds-info: List of absolute paths to CDS files.
   --validate-input : Validate FASTA file format (0 by default)
+  --validate-db : Validate created DB files (0 by default)
 ```
 This will generate **diffIdx**, **info**, **split**, and **taxID_list** and some other files. You can delete `*_diffIdx` and `*_info` files and `DATE-TIME` folder (e.g., `2025-1-24-10-32`) if generated.
 
@@ -465,6 +469,7 @@ metabuli updateDB <NEW DBDIR> <FASTA_LIST> <accession2taxid> <OLD DBDIR> [option
   --make-library : Make species library for faster execution (1 by default).
   --new-taxa : List of new taxa to be added.
   --validate-input : Validate FASTA file format (0 by default)
+  --validate-db : Validate created DB files (0 by default)
 ```
 
 #### \<Add sequences of new taxa> - Please refer [this section](#add-sequences-of-new-taxa).
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
@@ -214,7 +214,7 @@ jobs:
       #     %BUILD_SOURCESDIRECTORY%\build\metabuli\metabuli.bat version
       #   displayName: Setup Busybox
       - script: |
-          %CYGWIN_ROOT%\bin\bash.exe -cl "${BUILD_SOURCESDIRECTORY}/util/Metabuli-regression/run_regression.sh $(cygpath ${BUILD_SOURCESDIRECTORY}/build/build_sse2/src/metabuli.exe) $(cygpath ${BUILD_SOURCESDIRECTORY}/examples)"
+          %CYGWIN_ROOT%\bin\bash.exe -cl "${BUILD_SOURCESDIRECTORY}/util/Metabuli-regression/run_regression.sh $(cygpath ${BUILD_SOURCESDIRECTORY}/build/build_sse41/src/metabuli.exe) $(cygpath ${BUILD_SOURCESDIRECTORY}/examples)"
         displayName: Run Regression Suite
         condition: eq(variables['regression'], 1)
 
diff --git a/lib/fastq_utils/fastq.cpp b/lib/fastq_utils/fastq.cpp
@@ -62,7 +62,7 @@ gzFile fastq_open(const char* filename,const char *mode);
 static void fastq_close(gzFile fd);
 
 void fastq_print_version() {
-  fprintf(stderr,"fastq_utils %s\n",VERSION);
+  fprintf(stderr,"fastq_utils %s\n",FASTQ_UTIL_VERSION);
 }
 
 
diff --git a/lib/fastq_utils/fastq.h b/lib/fastq_utils/fastq.h
@@ -20,7 +20,7 @@
 #
 # =========================================================
 */
-#define VERSION "0.25.2"
+#define FASTQ_UTIL_VERSION "0.25.2"
 
 #define DEFAULT  0
 #define CASAVA18 1
diff --git a/lib/mmseqs b/lib/mmseqs
@@ -1 +1 @@
-Subproject commit 995f37678532fcaccd65020912b2493e9011eb91
+Subproject commit 62f8d3e17b429825695273e347775871339cda75
diff --git a/src/LocalCommandDeclarations.h b/src/LocalCommandDeclarations.h
@@ -26,5 +26,6 @@ extern int accession2taxid(int argc, const char **argv, const Command& command);
 extern int editNames(int argc, const char **argv, const Command& command);
 extern int createnewtaxalist(int argc, const char **argv, const Command& command);
 extern int classifiedRefiner(int argc, const char **argv, const Command& command);
+extern int validateDatabase(int argc, const char **argv, const Command& command);
 
 #endif //ADCLASSIFIER2_LOCALCOMMANDDECLARATIONS_H
diff --git a/src/MetabuliBase.cpp b/src/MetabuliBase.cpp
@@ -20,16 +20,16 @@ std::vector<Command> metabuliCommands = {
           {"sequenceDB", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::sequenceDb },
           {"tmpDir", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::directory }}},
         {"build", build, &localPar.build, COMMAND_DATABASE_CREATION,
-                "Build database based on the list of FASTA files.",
+                "Build database based on the list of FASTA files",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
-                "<DB dir> <FASTA list> <Accesssion2taxid>",
+                "<database directory> <FASTA list> <accesssion2taxid>",
                 CITATION_SPACEPHARER,
                 {{"Directory where the DB will be generated", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::empty},
                  {"A list of FASTA files", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile},
                  {"Mapping file (accession to tax ID)", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile}}},
         {"updateDB", updateDB, &localPar.updateDB, COMMAND_DATABASE_CREATION,
-                "Add new sequences to the database",
+                "Add new sequences to an existing database",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
                 " <new database directory> <FASTA list> <accesssion2taxid> <old database directory>",
@@ -59,15 +59,15 @@ std::vector<Command> metabuliCommands = {
                 "<info>",
                 CITATION_SPACEPHARER,
                 {{"info", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile}}},        
-        {"database-report", databaseReport, &localPar.databaseReport, COMMAND_DB,
-                "It generates a report of taxa in a database.",
+        {"database-report", databaseReport, &localPar.databaseReport, COMMAND_DATABASE_CREATION,
+                "Generate a taxonomy report of a database",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
-                "<i: DBDIR> ",
+                "<i: database directory> ",
                 CITATION_SPACEPHARER,
                 {{"Directory where the DB will be generated", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::directory}}},
         {"classify", classify, &localPar.classify, COMMAND_MAIN,
-                "Assigning taxonomy label to query reads",
+                "Assign taxonomic labels to query reads",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
                 "<i:query file(s)> <i:database directory> <o:output directory> <job ID> ",
@@ -76,8 +76,8 @@ std::vector<Command> metabuliCommands = {
                  {"DB dir", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::directory},
                  {"out dir", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::directory},
                  {"job ID", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::flatfile}}},
-        {"extract", extract, &localPar.extract, COMMAND_MAIN,
-                "It extracts reads classified into a certain taxon. It should be used after classification.",
+        {"extract", extract, &localPar.extract, COMMAND_FORMAT_CONVERSION,
+                "Extract reads classified to a certain taxon (Used after classification)",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
                 "<i:query file(s)> <i:read-by-read result> <i:database directory>",
@@ -114,7 +114,7 @@ std::vector<Command> metabuliCommands = {
         //              {"NCBI style accession2taxid file. It should be consistent to tax dump files.", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile},
         //              {"DB directory", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::directory}}},
         {"make-report", binning2report, &localPar.binning2report, COMMAND_FORMAT_CONVERSION,
-                "It generates Kraken style report file from binning results",
+                "Generate a Kraken-style taxonomy report using read-by-read classifications",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
                 "<i:Binning Result> <o:OUT DIR> <o:JOB ID> <i: TAXONOMY DIR> ",
@@ -183,13 +183,20 @@ std::vector<Command> metabuliCommands = {
                  {"accession2taxid", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile},
                  {"output prefix", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::flatfile}}},
         {"classifiedRefiner", classifiedRefiner, &localPar.classifiedRefiner, COMMAND_FORMAT_CONVERSION,
-                "Generates refined classification file from classified file",
+                "Refine read-by-read classification file",
                 nullptr,
                 "Siyoung Choi <remy0502@snu.ac.kr>",
                 "<i: classified file> <i: taxonomy dump>",
                 CITATION_SPACEPHARER,
                 {{"classified file", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile},
-                {"taxonomy dump", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::directory}}}
+                {"taxonomy dump", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::directory}}},
+        {"validatedb", validateDatabase, &localPar.validateDatabase, COMMAND_DATABASE_CREATION,
+                "Validate a database",
+                nullptr,
+                "Jaebeom Kim <jbeom0731@gmail.com>",
+                "<i: database directory>",
+                CITATION_SPACEPHARER,
+                {{"database directory", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::directory}}}
 };
 
 std::vector<KmerThreshold> externalThreshold = {};
diff --git a/src/commons/Classifier.cpp b/src/commons/Classifier.cpp
@@ -54,7 +54,7 @@ void Classifier::startClassify(const LocalParameters &par) {
 
         // Get splits for remaining sequences
         if (tries == 1) {
-                cout << "Indexing query file ...";
+                cout << "Indexing query file ..." << std::flush;
         }
         queryIndexer->setBytesPerKmer(matchPerKmer);
         queryIndexer->indexQueryFile(processedReadCnt);
@@ -124,15 +124,12 @@ void Classifier::startClassify(const LocalParameters &par) {
 
                 numOfTatalQueryKmerCnt += queryKmerBuffer.startIndexOfReserve;
             } else { // search was incomplete
-                // Increase matchPerKmer and try again 
+                // increase matchPerKmer and try again 
                 matchPerKmer += 4;
-                // delete kseq1;
-                // delete kseq2;
                 cout << "--match-per-kmer was increased to " << matchPerKmer << " and searching again..." << endl;
-                // cout << "The search was incomplete. Increasing --match-per-kmer to " << matchPerKmer << " and trying again..." << endl;
                 break;
             }
-         }
+        }
          
         delete kseq1;
         if (par.seqMode == 2) {
diff --git a/src/commons/KmerMatcher.cpp b/src/commons/KmerMatcher.cpp
@@ -26,7 +26,7 @@ KmerMatcher::~KmerMatcher() {
 }
 
 void KmerMatcher::loadTaxIdList(const LocalParameters & par) {
-    cout << "Loading the list for taxonomy IDs ... ";
+    cout << "Loading the list for taxonomy IDs ... " << std::flush;
     if (par.contamList != "") {
         vector<string> contams = Util::split(par.contamList, ",");
         for (auto &contam : contams) {
diff --git a/src/commons/LocalParameters.cpp b/src/commons/LocalParameters.cpp
@@ -39,6 +39,13 @@ LocalParameters::LocalParameters() :
                         typeid(int),
                         (void *) &skipRedundancy,
                         "[0-1]"),
+        VALIDATE_DB(VALIDATE_DB_ID,
+                     "--validate-db",
+                     "Validate the database",
+                     "Validate the database. It checks if all required files are present and if the k-mer count is consistent.",
+                     typeid(int),
+                     (void *) &validateDb,
+                     "[0-1]"),
         SEQ_MODE(SEQ_MODE_ID,
                  "--seq-mode",
                  "Sequencing type",
@@ -514,6 +521,7 @@ LocalParameters::LocalParameters() :
     build.push_back(&MAKE_LIBRARY);
     build.push_back(&GTDB);
     build.push_back(&VALIDATE_INPUT);
+    build.push_back(&VALIDATE_DB);
 
     // updateDB
     updateDB.push_back(&PARAM_THREADS);
@@ -529,6 +537,7 @@ LocalParameters::LocalParameters() :
     updateDB.push_back(&MAKE_LIBRARY);
     updateDB.push_back(&GTDB);
     updateDB.push_back(&VALIDATE_INPUT);
+    updateDB.push_back(&VALIDATE_DB);
 
     //classify
     classify.push_back(&PARAM_THREADS);
@@ -549,6 +558,7 @@ LocalParameters::LocalParameters() :
     classify.push_back(&SKIP_REDUNDANCY);
     classify.push_back(&PRINT_LINEAGE);
     classify.push_back(&VALIDATE_INPUT);
+    classify.push_back(&VALIDATE_DB);
 
     // extract
     extract.push_back(&TAXONOMY_PATH);
diff --git a/src/commons/LocalParameters.h b/src/commons/LocalParameters.h
@@ -43,6 +43,7 @@ class LocalParameters : public Parameters {
     std::vector<MMseqsParameter*> editNames;
     std::vector<MMseqsParameter*> createnewtaxalist;
     std::vector<MMseqsParameter*> classifiedRefiner;
+    std::vector<MMseqsParameter*> validateDatabase;
     
 
     // Superkingdom taxonomy id
@@ -52,6 +53,7 @@ class LocalParameters : public Parameters {
 
     // DB and classify
     PARAMETER(SKIP_REDUNDANCY)
+    PARAMETER(VALIDATE_DB)
 
     // Classify
     PARAMETER(SEQ_MODE)
@@ -136,6 +138,7 @@ class LocalParameters : public Parameters {
 
     // DB and classify
     int skipRedundancy;
+    int validateDb;
 
     // Classify
     int seqMode;
diff --git a/src/commons/Match.h b/src/commons/Match.h
@@ -26,7 +26,7 @@ struct Match { // 24 byte
 
     void printMatch() const {
         std::cout << qInfo.sequenceID << " " << qInfo.pos << " " << qInfo.frame << " "
-        << targetId << " " << speciesId << " " << rightEndHamming << " " << (int)hamming << " " << getScore() << std::endl;
+        << targetId << " " << speciesId << " " << rightEndHamming << " " << (int)hamming << " " << getScore() << "\n";
     }
 
     float getScore(float score = 0.0f, int cnt = 0) const { 
diff --git a/src/commons/QueryIndexer.cpp b/src/commons/QueryIndexer.cpp
@@ -22,8 +22,8 @@ QueryIndexer::QueryIndexer(const LocalParameters & par) {
 }
 
 void QueryIndexer::setAvailableRam() {
-    availableRam = (((size_t) maxRam * 1024 * 1024 * 1024) * 85) / 100
-               - ((size_t) 128 * 1024 * 1024 * threads);
+    availableRam = ((size_t) maxRam * (size_t) 1024 * 1024 * 1024)
+                         - ((size_t) 128 * 1024 * 1024 * (size_t) threads);
     // std::cout << "availableRam: " << availableRam << "\n";
 }
 
@@ -47,7 +47,7 @@ void QueryIndexer::indexQueryFile(size_t processedQueryNum) {
         while (kseq->ReadEntry()) {
             seqPos++;
             if (unlikely(kseq->entry.sequence.l == 0 || kseq->entry.name.l == 0)) {
-                std::cout << seqPos << "th entry has no sequence or name" << std::endl;
+                std::cout << seqPos << "th entry has no sequence or name." << std::endl;
                 exit(1);
             }
             readNum_1++;
@@ -67,12 +67,6 @@ void QueryIndexer::indexQueryFile(size_t processedQueryNum) {
             }
         }
         querySplits.emplace_back(start, readNum_1, kmerCnt, seqCnt);
-        // Print elements
-        // for (auto & querySplit : querySplits) {
-        //     std::cout << "start: " << querySplit.start << "\t";
-        //     std::cout << "end: " << querySplit.end << "\t";
-        //     std::cout << "kmerCnt: " << querySplit.kmerCnt << "\n";
-        // }
         delete kseq;
     } else {
         KSeqWrapper* kseq_1 = KSeqFactory(queryPath_1.c_str());
@@ -98,7 +92,7 @@ void QueryIndexer::indexQueryFile(size_t processedQueryNum) {
                 seqPos1++;
                 if (unlikely(kseq_1->entry.sequence.l == 0 || kseq_1->entry.name.l == 0)) {
                     std::cout << "In file " << queryPath_1 << ", " << std::endl << "\t";
-                    std::cout << seqPos1 << "th entry has no sequence or name" << std::endl;
+                    std::cout << seqPos1 << "th entry has no sequence or name." << std::endl;
                     exit(1);
                 }
                 readNum_1++;
@@ -113,7 +107,7 @@ void QueryIndexer::indexQueryFile(size_t processedQueryNum) {
                 seqPos2++;
                 if (unlikely(kseq_2->entry.sequence.l == 0 || kseq_2->entry.name.l == 0)) {
                     std::cout << "In file " << queryPath_2 << ", " << std::endl << "\t";
-                    std::cout << seqPos2 << "th entry has no sequence or name" << std::endl;
+                    std::cout << seqPos2 << "th entry has no sequence or name." << std::endl;
                     exit(1);
                 }
                 readNum_2++;
diff --git a/src/commons/Reporter.cpp b/src/commons/Reporter.cpp
@@ -26,11 +26,14 @@ void Reporter::openReadClassificationFile() {
 }
 
 void Reporter::writeReadClassification(const vector<Query> & queryList, bool classifiedOnly) {
-    readClassificationFile << "#is_classified\tname\ttaxID\tquery_length\tscore\trank";
-    if (par.printLineage) {
-        readClassificationFile << "\tlineage";
+    if (isFirstTime) {
+        readClassificationFile << "#query_id\tname\ttaxID\tquery_length\tscore\trank";
+        if (par.printLineage) {
+            readClassificationFile << "\tlineage";
+        }
+        readClassificationFile << "\ttaxID:match_count\n";
+        isFirstTime = false;
     }
-    readClassificationFile << "\ttaxID:match_count\n";
     for (size_t i = 0; i < queryList.size(); i++) {
         if (classifiedOnly && !queryList[i].isClassified) {
             continue;
diff --git a/src/commons/Reporter.h b/src/commons/Reporter.h
@@ -22,6 +22,8 @@ class Reporter {
     string readClassificationFileName;
     ofstream readClassificationFile;
 
+    bool isFirstTime = true;
+
 public:
     Reporter(const LocalParameters &par, TaxonomyWrapper *taxonomy, const std::string &customReportFileName = "");
     // Write report
diff --git a/src/metabuli.cpp b/src/metabuli.cpp
@@ -3,7 +3,7 @@
 
 const char* binary_name = "metabuli";
 const char* tool_name = "metabuli";
-const char* tool_introduction = "Metabuli is a taxonomical classifier that jointly analyzes amino acid and DNA sequences.";
+const char* tool_introduction = "Metabuli is a metagenomic taxonomic classifier that jointly analyzes amino acid and DNA sequences.";
 const char* main_author = "Jaebeom Kim <jbeom0731@gmail.com>";
 const char* show_extended_help = "1";
 const char* show_bash_info = nullptr;
@@ -15,7 +15,19 @@ void (*validatorUpdate)(void) = 0;
 
 extern std::vector<Command> baseCommands;
 extern std::vector<Command> metabuliCommands;
+extern std::vector<Categories> categories;
 void init() {
+        for (size_t i = 0; i < categories.size(); ++i) {
+            if (categories[i].mode & COMMAND_MAIN) {
+                categories[i].title = "Classification workflows";
+            }
+            if (categories[i].mode & COMMAND_DATABASE_CREATION) {
+                categories[i].title = "Database workflows";
+            }
+            if (categories[i].mode & COMMAND_FORMAT_CONVERSION) {
+                categories[i].title = "Downstream workflows";
+            }
+        }
     registerCommands(&metabuliCommands);
 }
 void (*initCommands)(void) = init;
diff --git a/src/util/CMakeLists.txt b/src/util/CMakeLists.txt
@@ -21,4 +21,6 @@ set(util_source_files
         util/createnewtaxalist.cpp
         util/createnewtaxalist.h
         util/classifiedRefiner.cpp
+        util/validateDatabase.cpp
+        util/validateDatabase.h
         PARENT_SCOPE)
diff --git a/src/util/validateDatabase.cpp b/src/util/validateDatabase.cpp
diff --git a/src/util/validateDatabase.h b/src/util/validateDatabase.h
diff --git a/src/workflow/build.cpp b/src/workflow/build.cpp
diff --git a/src/workflow/classify.cpp b/src/workflow/classify.cpp
diff --git a/src/workflow/updateDB.cpp b/src/workflow/updateDB.cpp
diff --git a/util/build_windows.sh b/util/build_windows.sh

Original file line number	Diff line number	Diff line change
`@@ -62,7 +62,7 @@ gzFile fastq_open(const char* filename,const char *mode);`
`62`	`62`	`static void fastq_close(gzFile fd);`
`63`	`63`
`64`	`64`	`void fastq_print_version() {`
`65`		`- fprintf(stderr,"fastq_utils %s\n",VERSION);`
	`65`	`+ fprintf(stderr,"fastq_utils %s\n",FASTQ_UTIL_VERSION);`
`66`	`66`	`}`
`67`	`67`
`68`	`68`
Original file line number	Diff line number	Diff line change
`@@ -20,7 +20,7 @@`
`20`	`20`	`#`
`21`	`21`	`# =========================================================`
`22`	`22`	`*/`
`23`		`-#define VERSION "0.25.2"`
	`23`	`+#define FASTQ_UTIL_VERSION "0.25.2"`
`24`	`24`
`25`	`25`	`#define DEFAULT 0`
`26`	`26`	`#define CASAVA18 1`
Original file line number	Diff line number	Diff line change
`@@ -26,7 +26,7 @@ KmerMatcher::~KmerMatcher() {`
`26`	`26`	`}`
`27`	`27`
`28`	`28`	`void KmerMatcher::loadTaxIdList(const LocalParameters & par) {`
`29`		`- cout << "Loading the list for taxonomy IDs ... ";`
	`29`	`+ cout << "Loading the list for taxonomy IDs ... " << std::flush;`
`30`	`30`	`if (par.contamList != "") {`
`31`	`31`	`vector<string> contams = Util::split(par.contamList, ",");`
`32`	`32`	`for (auto &contam : contams) {`