Merge read-group branch into master (#165)

* first compiling version after merging read-grouping module

* now 1-based temporarily

* implement SyncmerScanner_dna2aa

* optimize read grouping

* cadd parameters for ceateingcommon k-mer list

* extract reference k-mers from six-frame translations

* extract reference k-mers from six-frame translations

* code clean

* KmerDbReader only reads k-mer values not IDs

* add parameters to read grouping module

* use only k-mer values during common k-mer filtering

* add yxml into lib

* WIP: making read grouping flexible to original classification format

* fix read index error in reporting updated results

* fix tiny-size db creation error

* regression test updated

* hide developing features

* fix errors in benchmark set generation

* fix issue 108

* update database moudle

Author: jaebeom-kim

data/metabulidatabases.sh

@@ -92,7 +92,7 @@ mkdir -p "${OUTDB}"
 case "${SELECTION}" in
     "GTDB")
         if notExists "${TMP_PATH}/gtdb.tar.gz"; then
-            downloadFile "https://metabuli.steineggerlab.workers.dev/gtdb.tar.gz" "${TMP_PATH}/gtdb.tar.gz"
+            downloadFile "https://metabuli.steineggerlab.workers.dev/gtdb+virus+human.tar.gz" "${TMP_PATH}/gtdb.tar.gz"
         fi
         tar zxvf "${TMP_PATH}/gtdb.tar.gz" -C "${OUTDB}"
        # push_back "${TMP_PATH}/gtdb"

src/CMakeLists.txt

@@ -3,6 +3,7 @@ include_directories(workflow)
 include_directories(util)
 include_directories(benchmark)
 include_directories(uniref)
+include_directories(read-group)
 include_directories(../lib/prodigal)
 include_directories(../lib/fasta_validator)
 include_directories(../lib/fastq_utils)
@@ -14,13 +15,15 @@ add_subdirectory(version)
 add_subdirectory(workflow)
 add_subdirectory(benchmark)
 add_subdirectory(uniref)
+add_subdirectory(read-group)
 
 add_executable(metabuli
         ${commons_source_files}
         ${util_source_files}
         ${workflow_source_files}
         ${benchmark_source_files}
         ${uniref_source_files}
+        ${read_group_source_files}
         metabuli.cpp
         LocalCommandDeclarations.h 
         util/filter_by_genus.cpp

src/MetabuliBase.cpp

@@ -28,7 +28,7 @@ std::vector<Command> metabuliCommands = {
                 {{"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}}},
-        {"create-uniref-db", create_unirefdb, &localPar.buildUnirefDb, COMMAND_DATABASE_CREATION,
+        {"create-uniref-db", create_unirefdb, &localPar.buildUnirefDb, COMMAND_EXPERT,
                 "Build a UniRef k-mer database",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
@@ -37,15 +37,15 @@ std::vector<Command> metabuliCommands = {
                 {{"Database directory", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::empty},
                  {"UniRef100 FASTA", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile},
                  {"UniRef tree DB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile}}},
-        {"create-uniref-tree", create_uniref_tree, &localPar.buildUnirefTree, COMMAND_DATABASE_CREATION,
+        {"create-uniref-tree", create_uniref_tree, &localPar.buildUnirefTree, COMMAND_EXPERT,
                 "Build a tree of UniRef100/90/50 clusters",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
                 "<database directory> <UniRef100 XML file>",
                 CITATION_SPACEPHARER,
                 {{"Database directory", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::empty},
                  {"UniRef100 XML file", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile}}},
-        {"create-common-kmer-list", create_common_kmer_list, &localPar.build, COMMAND_EXPERT,
+        {"create-common-kmer-list", create_common_kmer_list, &localPar.createCommonKmerList, COMMAND_EXPERT,
                 "Create a list of k-mers that appear in multiple species",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
@@ -118,7 +118,7 @@ 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}}},
-        {"assign_uniref", assign_uniref, &localPar.assignUniref, COMMAND_MAIN,
+        {"assign_uniref", assign_uniref, &localPar.assignUniref, COMMAND_EXPERT,
                 "Assign UniRef cluster labels to query protein sequences.",
                 nullptr,
                 "Jaebeom Kim <jbeom0731@gmail.com>",
@@ -274,6 +274,17 @@ std::vector<Command> metabuliCommands = {
                 {{"Assembly accessions", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile},
                  {"Taxonomy dump", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::directory},
                  {"Assembly accession to taxid mapping", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile}}},
+        // {"grouping", groupGeneration, &localPar.groupGeneration, COMMAND_MAIN,
+        //         "Assigning group representative taxonomy label to unclassified query reads",
+        //         nullptr,
+        //         "Luna Sung-eun Jang <jse9512@gmail.com>",
+        //         "<i:query file(s)> <i:common k-mer database> <i: taxonomy directory> <i: read-by-read result> <o: output directory>", // TODO: follow this
+        //         CITATION_SPACEPHARER,
+        //         {{"FASTA/Q", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA | DbType::VARIADIC, &DbValidator::flatfile},
+        //          {"common k-mer database", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::directory},
+        //          {"taxonomy directory", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::directory},
+        //          {"read-by-read result", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::flatfile},
+        //          {"output directory", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::directory}}},
         {"count-common-kmers", count_common_kmers, &localPar.expand_diffidx, COMMAND_EXPERT,
                 "Count k-mers common at each taxonomy level",
                 nullptr,
@@ -296,7 +307,7 @@ std::vector<DatabaseDownload> externalDownloads = {
         },
         {
                 "GTDB",
-                "GTDB 214.1 (Complete/Chromosome level only, CheckM completeness > 90 and contamination < 5) and a human genome (T2T-CHM13v2.0)",
+                "GTDB 220 (Complete/Chromosome level only, CheckM completeness > 90 and contamination < 5), a human genome (T2T-CHM13v2.0), RefSeq viruses",
                 "Donovan et al. GTDB: an ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy. Nucleic Acids Res. (2022)",
                 "https://gtdb.ecogenomic.org/",
                 true, LocalParameters::DBTYPE_INDEX_DB, metabulidatabases_sh, metabulidatabases_sh_len,

src/benchmark/makeBenchmarkSet.cpp

@@ -1,4 +1,5 @@
 #include "IndexCreator.h"
+#include "editNames.h"
 #include <iostream>
 #include <istream>
 #include <string>
@@ -43,6 +44,8 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
     string databaseAssemblyList = assemblyList + ".databaseAssembly";
     string totalExludedAssemblyList = assemblyList + ".totalExcludedAssembly";
 
+    // editNames(taxonomyPath + "/names.dmp", par.assacc2taxid);
+
     // Load taxonomy
     TaxonomyWrapper taxonomy(taxonomyPath + "/names.dmp",
                              taxonomyPath + "/nodes.dmp",
@@ -52,18 +55,21 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
 
     cout << "Making name2taxid map...";
     std::unordered_map<std::string, TaxID> name2InternalTaxId;
+    std::unordered_map<std::string, TaxID> tempMap;
     taxonomy.getName2InternalTaxid(name2InternalTaxId);
     for (auto &it : name2InternalTaxId) {
         if (it.first.find(".") == std::string::npos) {
             continue;
         }
         string accessionNoVersion = it.first.substr(0, it.first.find("."));
-        name2InternalTaxId[accessionNoVersion] = it.second;
+        tempMap[accessionNoVersion] = it.second;
+    }
+    for (auto &it : tempMap) {
+        name2InternalTaxId[it.first] = it.second;
     }
     cout << "done." << endl;
 
-    cout << "Making observedAcc2taxid map...";
-    std::unordered_map<std::string, TaxID> observedAcc2taxid;
+
     cout << "Loading assembly list...";
     vector<string> totalAssemblyAccessions;
     vector<Assembly> assemblies;
@@ -72,8 +78,11 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
         cerr << "Error: could not open assembly list file " << assemblyList << endl;
         return 1;
     }
+
     string assemblyAccession;
     TaxID taxid;
+    cout << "Making observedAcc2taxid map...";
+    std::unordered_map<std::string, TaxID> observedAcc2taxid;
     while (getline(assemblyListFile, assemblyAccession)) {
         string assAccNoVersion = assemblyAccession.substr(0, assemblyAccession.find("."));
 
@@ -88,7 +97,7 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
         } else if (name2InternalTaxId.find(assAccNoVersion) != name2InternalTaxId.end()) {
             taxid = name2InternalTaxId[assAccNoVersion];
         } else {
-            cerr << "Error: accession " << assemblyAccession << " not found in the taxonomy" << endl;
+            cerr << "Error: accession " << assemblyAccession << " " << assAccNoVersion << " not found in the taxonomy" << endl;
             return 1;
         }
         observedAcc2taxid[assemblyAccession] = taxid;
@@ -257,7 +266,7 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
     }
 
     // Select n genera with multiple species
-    n = int(genusWithMultipleSpecies.size() / 4);
+    n = int(genusWithMultipleSpecies.size() / 3);
     cout << "Excluding " << n << " species" << endl;
     vector<TaxID> selectedGenera;
     for (int i = 0; i < n; i++) {
@@ -314,7 +323,7 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
     }
 
     // 2. Select n species with multiple assemblies
-    n = int(speciesWithMultipleAssemblies.size()/2);
+    n = int(speciesWithMultipleAssemblies.size()/3);
     cout << "Excluding " << n << " assemblies" << endl;
     vector<TaxID> selectedSpecies;
     for (int i = 0; i < n; i++) {
@@ -409,17 +418,23 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
             const string & rank = taxonomy.getString(lcaTaxon->rankIdx);
             if (rank == "order") {
                 orderCount++;
-            } else if (taxonomy.findRankIndex(rank) == -1) {
-                cout << "Error: LCA rank index is -1." << endl;
-                return 1;
-            } else if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("order")) {
-                cout << "Error: " << excludedGenusAssemblies[i] << " is not a valid family exclusion. LCA is below order rank." << endl;
-                cout << excludedGenusAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
+                continue;
+            } 
+            
+            if (strcmp(taxonomy.getString(lcaTaxon->nameIdx), "root") == 0) {
+                continue;
+            }
+
+            if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("order")) {    
+                cout << "Error: " << excludedFamilyAssemblies[i] << " is not a valid family exclusion. LCA is below order rank." << endl;
+                cout << "Name: " << taxonomy.getString(lcaTaxon->nameIdx) << " TaxID: " << lcaTaxon->taxId << endl;
+                cout << "Rank: " << rank << " RankIdx: " << lcaTaxon->rankIdx << endl;
+                cout << excludedFamilyAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
                 return 1;
             }
         }
         if (orderCount == 0) {
-            cout << "Error: " << excludedGenusAssemblies[i] << " is not a valid exclusion. No order rank LCA." << endl;
+            cout << "Error: " << excludedFamilyAssemblies[i] << " is not a valid exclusion. No order rank LCA." << endl;
             return 1;
         }
     }
@@ -438,14 +453,28 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
             const string & rank = taxonomy.getString(lcaTaxon->rankIdx);
             if (rank == "family") {
                 familyCount++;
-            } else if (taxonomy.findRankIndex(rank) == -1) {
-                cout << "Error: LCA rank index is -1." << endl;
-                return 1;
-            } else if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("family")) {
-                cout << "Error: " << excludedGenusAssemblies[i] << " is not a valid exclusion. LCA is below Family rank." << endl;
-                cout << excludedGenusAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
-                return 1;
+                continue;
+            } 
+
+            if (strcmp(taxonomy.getString(lcaTaxon->nameIdx), "root") == 0) {
+                continue;
             }
+
+            if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("family")) {    
+                cout << "Error: " << excludedFamilyAssemblies[i] << " is not a valid genus exclusion. LCA is below family rank." << endl;
+                cout << "Name: " << taxonomy.getString(lcaTaxon->nameIdx) << " TaxID: " << lcaTaxon->taxId << endl;
+                cout << "Rank: " << rank << " RankIdx: " << lcaTaxon->rankIdx << endl;
+                cout << excludedFamilyAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
+                return 1;
+            }            
+            // else if (taxonomy.findRankIndex(rank) == -1) {
+            //     cout << "Error: LCA rank index is -1." << endl;
+            //     return 1;
+            // } else if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("family")) {
+            //     cout << "Error: " << excludedGenusAssemblies[i] << " is not a valid exclusion. LCA is below Family rank." << endl;
+            //     cout << excludedGenusAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
+            //     return 1;
+            // }
         }
         if (familyCount == 0) {
             cout << "Error: " << excludedGenusAssemblies[i] << " is not a valid exclusion. No Family rank LCA." << endl;
@@ -465,16 +494,34 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
             TaxID databaseTaxid = observedAcc2taxid[databaseAssemblies[j]];
             const TaxonNode * lcaTaxon = taxonomy.taxonNode(taxonomy.LCA(excludedTaxid, databaseTaxid));
             const string & rank = taxonomy.getString(lcaTaxon->rankIdx);
+
             if (rank == "genus") {
                 genusCount++;
-            } else if (taxonomy.findRankIndex(rank) == -1) {
-                cout << "Error: LCA rank index is -1." << endl;
-                return 1;
-            } else if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("genus")) {
-                cout << "Error: " << excludedSpeciesAssemblies[i] << " is not a valid exclusion. LCA is below Genus rank." << endl;
-                cout << excludedSpeciesAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
-                return 1;
+                continue;
+            } 
+
+            if (strcmp(taxonomy.getString(lcaTaxon->nameIdx), "root") == 0) {
+                continue;
             }
+
+            if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("genus")) {    
+                cout << "Error: " << excludedFamilyAssemblies[i] << " is not a valid species exclusion. LCA is below genus rank." << endl;
+                cout << "Name: " << taxonomy.getString(lcaTaxon->nameIdx) << " TaxID: " << lcaTaxon->taxId << endl;
+                cout << "Rank: " << rank << " RankIdx: " << lcaTaxon->rankIdx << endl;
+                cout << excludedFamilyAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
+                return 1;
+            }            
+
+            // if (rank == "genus") {
+            //     genusCount++;
+            // } else if (taxonomy.findRankIndex(rank) == -1) {
+            //     cout << "Error: LCA rank index is -1." << endl;
+            //     return 1;
+            // } else if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("genus")) {
+            //     cout << "Error: " << excludedSpeciesAssemblies[i] << " is not a valid exclusion. LCA is below Genus rank." << endl;
+            //     cout << excludedSpeciesAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
+            //     return 1;
+            // }
         }
         if (genusCount == 0) {
             cout << "Error: " << excludedSpeciesAssemblies[i] << " is not a valid exclusion. No Genus rank LCA." << endl;
@@ -496,11 +543,20 @@ int makeGtdbBenchmarkSet(const LocalParameters & par) {
             const string & rank = taxonomy.getString(lcaTaxon->rankIdx);
             if (rank == "species") {
                 speciesCount++;
-            } else if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("species")) {
-                cout << "Error: " << excludedSubspeciesAssemblies[i] << " is not a valid exclusion. LCA is below Species rank." << endl;
-                cout << excludedSubspeciesAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
-                return 1;
+                continue;
+            } 
+
+            if (strcmp(taxonomy.getString(lcaTaxon->nameIdx), "root") == 0) {
+                continue;
             }
+
+            if (taxonomy.findRankIndex(rank) < taxonomy.findRankIndex("species")) {    
+                cout << "Error: " << excludedFamilyAssemblies[i] << " is not a valid exclusion. LCA is below species rank." << endl;
+                cout << "Name: " << taxonomy.getString(lcaTaxon->nameIdx) << " TaxID: " << lcaTaxon->taxId << endl;
+                cout << "Rank: " << rank << " RankIdx: " << lcaTaxon->rankIdx << endl;
+                cout << excludedFamilyAssemblies[i] << " " << databaseAssemblies[j] << " " << rank << endl;
+                return 1;
+            } 
         }
         if (speciesCount == 0) {
             cout << "Error: " << excludedSubspeciesAssemblies[i] << " is not a valid exclusion. No Species rank LCA." << endl;