Using SSHash submodule for graph representation

.gitignore

@@ -1,5 +1,6 @@
 *.fa
 *.fai
+.vscode
 !metagraph/tests/data/*.fa
 !metagraph/tests/data/*.fai
 metagraph/tests/data/*dump_test*

.gitmodules

@@ -68,3 +68,6 @@
 [submodule "metagraph/external-libraries/htslib"]
 	path = metagraph/external-libraries/htslib
 	url = https://github.com/samtools/htslib
+[submodule "metagraph/external-libraries/sshash"]
+	path = metagraph/external-libraries/sshash
+	url = https://github.com/ratschlab/sshash

metagraph/CMakeLists.txt

@@ -277,6 +277,10 @@ include_directories(
   PUBLIC SYSTEM
   external-libraries/asio/asio/include
   external-libraries/simde-no-tests
+  external-libraries/sshash
+  external-libraries/sshash/include
+  external-libraries/sshash/external/pthash
+  external-libraries/sshash/external/pthash/external/essentials/include
 )
 
 add_library(ips4o INTERFACE)
@@ -319,6 +323,7 @@ IF(APPLE)
     set(THREADS_PREFER_PTHREAD_FLAG ON)
 ENDIF()
 add_subdirectory(external-libraries/spdlog)
+add_subdirectory(external-libraries/sshash SYSTEM)
 add_subdirectory(external-libraries/DYNAMIC)
 add_subdirectory(external-libraries/zlib)
 target_compile_options(zlib
@@ -409,6 +414,7 @@ set(METALIBS ${METALIBS}
   mersenne_twister
   sdust
   spdlog::spdlog
+  sshash_static
   XXSDS_DYNAMIC
   ips4o
   caches

metagraph/src/cli/build.cpp

@@ -7,6 +7,7 @@
 #include "graph/representation/hash/dbg_hash_ordered.hpp"
 #include "graph/representation/hash/dbg_hash_string.hpp"
 #include "graph/representation/hash/dbg_hash_fast.hpp"
+#include "graph/representation/hash/dbg_sshash.hpp"
 #include "graph/representation/bitmap/dbg_bitmap.hpp"
 #include "graph/representation/bitmap/dbg_bitmap_construct.hpp"
 #include "graph/representation/succinct/dbg_succinct.hpp"
@@ -249,6 +250,13 @@ int build_graph(Config *config) {
                                                              get_verbose()));
         }
 
+    } else if (config->graph_type == Config::GraphType::SSHASH && !config->dynamic) {
+        graph.reset(new DBGSSHash(files.at(0), config->k, config->graph_mode, config->num_chars));
+        if (files.size() > 1) {
+            logger->error("DBGSSHash does not support multiple input files.");
+            exit(1);
+        }
+
     } else {
         //slower method
         switch (config->graph_type) {
@@ -284,6 +292,11 @@ int build_graph(Config *config) {
                               " in dynamic regime is not supported");
                 exit(1);
 
+            case Config::GraphType::SSHASH:
+                logger->error("SSHash-graph construction"
+                              " in dynamic regime is not supported");
+                exit(1);
+
             case Config::GraphType::INVALID:
                 assert(false);
         }

metagraph/src/cli/config/config.cpp

@@ -243,6 +243,8 @@ Config::Config(int argc, char *argv[]) {
             alignment_post_chain = true;
         } else if (!strcmp(argv[i], "--align-no-seed-complexity-filter")) {
             alignment_seed_complexity_filter = false;
+        } else if (!strcmp(argv[i], "--num-chars")) {
+            num_chars = atoll(get_value(i++));
         } else if (!strcmp(argv[i], "--max-hull-depth")) {
             max_hull_depth = atoll(get_value(i++));
         } else if (!strcmp(argv[i], "--batch-align")) {
@@ -815,6 +817,8 @@ Config::GraphType Config::string_to_graphtype(const std::string &string) {
     } else if (string == "bitmap") {
         return GraphType::BITMAP;
 
+    } else if (string == "sshash") {
+        return GraphType::SSHASH;
     } else {
         std::cerr << "Error: unknown graph representation" << std::endl;
         exit(1);
@@ -956,12 +960,13 @@ if (advanced) {
             fprintf(stderr, "\t   --max-count-q [INT] \tmax k-mer abundance quantile (max-count is used by default) [1.0]\n");
             fprintf(stderr, "\t   --reference [STR] \tbasename of reference sequence (for parsing VCF files) []\n");
             fprintf(stderr, "\n");
-            fprintf(stderr, "\t   --graph [STR] \tgraph representation: succinct / bitmap / hash / hashstr / hashfast [succinct]\n");
+            fprintf(stderr, "\t   --graph [STR] \tgraph representation: succinct / bitmap / hash / hashstr / hashfast [succinct] / sshash\n");
             fprintf(stderr, "\t   --state [STR] \tstate of succinct graph: small / dynamic / stat / fast [stat]\n");
             fprintf(stderr, "\t   --inplace \t\tconstruct succinct graph in-place and serialize without loading to RAM [off]\n");
             fprintf(stderr, "\t   --count-kmers \tcount k-mers and build weighted graph [off]\n");
             fprintf(stderr, "\t   --count-width \tnumber of bits used to represent k-mer abundance [8]\n");
             fprintf(stderr, "\t   --index-ranges [INT]\tindex all node ranges in BOSS for suffixes of given length [%zu]\n", kDefaultIndexSuffixLen);
+            fprintf(stderr, "\t   --num-chars [INT]\tif the number of characters is known beforehand, enter it here [0]\n");
             fprintf(stderr, "\t-k --kmer-length [INT] \tlength of the k-mer to use [3]\n");
 #if ! _PROTEIN_GRAPH
             fprintf(stderr, "\t   --mode \t\tk-mer indexing mode: basic / canonical / primary [basic]\n");

metagraph/src/cli/config/config.hpp

@@ -104,6 +104,7 @@ class Config {
     unsigned long long int num_rows_subsampled = 1'000'000;
     unsigned long long int num_singleton_kmers = 0;
     unsigned long long int max_hull_depth = -1;  // the default is a function of input
+    unsigned long long int num_chars = 0;
 
     uint8_t count_width = 8;
 
@@ -223,6 +224,7 @@ class Config {
         HASH_PACKED,
         HASH_STR,
         HASH_FAST,
+        SSHASH,
         BITMAP,
     };
 

metagraph/src/cli/load/load_graph.cpp

@@ -7,6 +7,7 @@
 #include "graph/representation/hash/dbg_hash_ordered.hpp"
 #include "graph/representation/hash/dbg_hash_string.hpp"
 #include "graph/representation/hash/dbg_hash_fast.hpp"
+#include "graph/representation/hash/dbg_sshash.hpp"
 #include "graph/representation/bitmap/dbg_bitmap.hpp"
 #include "graph/representation/succinct/dbg_succinct.hpp"
 #include "cli/config/config.hpp"
@@ -36,6 +37,9 @@ Config::GraphType parse_graph_type(const std::string &filename) {
     } else if (utils::ends_with(filename, graph::DBGBitmap::kExtension)) {
         return Config::GraphType::BITMAP;
 
+    } else if (utils::ends_with(filename, graph::DBGSSHash::kExtension)) {
+        return Config::GraphType::SSHASH;
+
     } else {
         return Config::GraphType::INVALID;
     }
@@ -58,10 +62,10 @@ std::shared_ptr<DeBruijnGraph> load_critical_dbg(const std::string &filename) {
 
         case Config::GraphType::HASH_FAST:
             return load_critical_graph_from_file<DBGHashFast>(filename);
-
         case Config::GraphType::BITMAP:
             return load_critical_graph_from_file<graph::DBGBitmap>(filename);
-
+        case Config::GraphType::SSHASH:
+            return load_critical_graph_from_file<graph::DBGSSHash>(filename);
         case Config::GraphType::INVALID:
             logger->error("Cannot load graph from file '{}', needs a valid file extension",
                           filename);

metagraph/src/graph/alignment/aligner_labeled.cpp

@@ -415,7 +415,7 @@ void LabeledExtender::call_alignments(score_t end_score,
                     }
 
                     if (label_it == label_end_it)
-                        throw std::bad_function_call();
+                        throw early_term();
 
                     if (c < *label_it)
                         return;
@@ -431,7 +431,7 @@ void LabeledExtender::call_alignments(score_t end_score,
                     }
                 }
             );
-        } catch (const std::bad_function_call&) {}
+        } catch (const early_term&) {}
     }
 
     if (alignment.label_coordinates.empty())

metagraph/src/graph/alignment/dbg_aligner.cpp

@@ -578,7 +578,7 @@ ::align_both_directions(std::string_view forward,
                 [&](Chain&& chain, score_t score) {
                     if (config_.num_alternative_paths <= 1
                             && finished_columns.size() == all_columns.size()) {
-                        throw std::bad_function_call();
+                        throw early_term();
                     }
 
                     double exact_match_fraction
@@ -595,7 +595,7 @@ ::align_both_directions(std::string_view forward,
 #endif
 
                     if (exact_match_fraction < config_.min_exact_match)
-                        throw std::bad_function_call();
+                        throw early_term();
 
                     extend_chain(chain[0].first.get_orientation() ? reverse : forward,
                                  chain[0].first.get_orientation() ? forward : reverse,
@@ -613,7 +613,7 @@ ::align_both_directions(std::string_view forward,
                 [&](Alignment::Column column) { return finished_columns.count(column); }
             );
             num_explored_nodes += this_num_explored;
-        } catch (const std::bad_function_call&) {}
+        } catch (const early_term&) {}
 
         for (Alignment &alignment : aggregator.get_alignments()) {
             if (alignment.get_score() < get_min_path_score(alignment))

metagraph/src/graph/representation/base/dbg_wrapper.hpp

@@ -73,7 +73,10 @@ class DBGWrapper : public DeBruijnGraph {
                             const std::function<bool()> &stop_early
                                 = [](){ return false; }) const override = 0;
 
-    virtual void call_kmers(const std::function<void(node_index, const std::string&)> &callback) const override = 0;
+    virtual void call_kmers(const std::function<void(node_index,
+                            const std::string&)> &callback,
+                            const std::function<bool()> &stop_early
+                                = [](){ return false; }) const override = 0;
 
   protected:
     std::shared_ptr<const Graph> graph_;

metagraph/src/graph/representation/base/sequence_graph.cpp

@@ -80,9 +80,33 @@ bool DeBruijnGraph::find(std::string_view sequence,
     return num_kmers_missing <= max_kmers_missing;
 }
 
-bool DeBruijnGraph::operator==(const DeBruijnGraph &) const {
-    throw std::runtime_error("Not implemented");
-    return false;
+void DeBruijnGraph
+::adjacent_outgoing_nodes(node_index node,
+                          const std::function<void(node_index)> &callback) const {
+    assert(node > 0 && node <= max_index());
+    call_outgoing_kmers(node, [&](auto child, char) { callback(child); });
+}
+
+void DeBruijnGraph
+::adjacent_incoming_nodes(node_index node,
+                          const std::function<void(node_index)> &callback) const {
+    assert(node > 0 && node <= max_index());
+    call_incoming_kmers(node, [&](auto parent, char) { callback(parent); });
+}
+
+bool DeBruijnGraph::operator==(const DeBruijnGraph &other) const {
+    bool not_equal = false;
+    other.call_kmers([&](node_index, const std::string &kmer) {
+        if (!find(kmer))
+            not_equal = true;
+    }, [&]() { return not_equal; });
+
+    call_kmers([&](node_index, const std::string &kmer) {
+        if (!other.find(kmer))
+            not_equal = true;
+    }, [&]() { return not_equal; });
+
+    return !not_equal;
 }
 
 void DeBruijnGraph::traverse(node_index start,
@@ -400,7 +424,12 @@ void DeBruijnGraph::call_unitigs(const CallPath &callback,
  * Traverse graph and iterate over all nodes
  */
 void DeBruijnGraph
-::call_kmers(const std::function<void(node_index, const std::string&)> &callback) const {
+::call_kmers(const std::function<void(node_index,
+             const std::string&)> &callback,
+             const std::function<bool()> &stop_early) const {
+    if (stop_early())
+        return;
+
     sdsl::bit_vector visited(max_index() + 1, false);
     std::stack<std::pair<node_index, std::string>> nodes;
 
@@ -444,7 +473,7 @@ ::call_kmers(const std::function<void(node_index, const std::string&)> &callback
                 sequence.push_back(next_c);
             }
         }
-    });
+    }, stop_early);
 }
 
 void DeBruijnGraph::print(std::ostream &out) const {

metagraph/src/graph/representation/base/sequence_graph.hpp

@@ -7,7 +7,6 @@
 #include <iostream>
 #include <memory>
 
-
 namespace utils {
     std::string make_suffix(const std::string &str, const std::string &suffix);
 } // namespace utils
@@ -16,6 +15,11 @@ namespace utils {
 namespace mtg {
 namespace graph {
 
+class early_term : public std::exception {
+  public:
+    const char* what() const noexcept override { return "Finished iteration"; }
+};
+
 class SequenceGraph {
   public:
     // Node indexes [1,...,max_index], but only num_nodes of them are real.
@@ -195,7 +199,8 @@ class DeBruijnGraph : public SequenceGraph {
                               size_t min_tip_size = 1,
                               bool kmers_in_single_form = false) const;
 
-    virtual void call_kmers(const std::function<void(node_index, const std::string&)> &callback) const;
+    virtual void call_kmers(const std::function<void(node_index, const std::string&)> &callback,
+                            const std::function<bool()> &stop_early = [](){ return false; }) const;
 
     virtual size_t outdegree(node_index) const = 0;
     virtual bool has_single_outgoing(node_index node) const { return outdegree(node) == 1; }
@@ -217,6 +222,13 @@ class DeBruijnGraph : public SequenceGraph {
     virtual void call_incoming_kmers(node_index kmer,
                                      const IncomingEdgeCallback &callback) const = 0;
 
+    // Given a node index, call the target nodes of all edges outgoing from it.
+    virtual void adjacent_outgoing_nodes(node_index node,
+                                         const std::function<void(node_index)> &callback) const override;
+    // Given a node index, call the source nodes of all edges incoming to it.
+    virtual void adjacent_incoming_nodes(node_index node,
+                                         const std::function<void(node_index)> &callback) const override;
+
     // Check whether graph contains fraction of nodes from the sequence
     virtual bool find(std::string_view sequence, double discovery_fraction = 1) const;
 

metagraph/src/graph/representation/bitmap/dbg_bitmap.cpp

@@ -196,20 +196,6 @@ size_t DBGBitmap::indegree(node_index node) const {
             - kmers_.rank1(first.data());
 }
 
-void DBGBitmap::adjacent_outgoing_nodes(node_index node,
-                                        const std::function<void(node_index)> &callback) const {
-    assert(node > 0 && node <= num_nodes());
-
-    call_outgoing_kmers(node, [&](node_index child, char) { callback(child); });
-}
-
-void DBGBitmap::adjacent_incoming_nodes(node_index node,
-                                        const std::function<void(node_index)> &callback) const {
-    assert(node > 0 && node <= num_nodes());
-
-    call_incoming_kmers(node, [&](node_index parent, char) { callback(parent); });
-}
-
 DBGBitmap::node_index DBGBitmap::to_node(const Kmer &kmer) const {
     uint64_t index = kmer.data() + 1;
 

metagraph/src/graph/representation/bitmap/dbg_bitmap.hpp

@@ -57,13 +57,6 @@ class DBGBitmap : public DeBruijnGraph {
     // Traverse the incoming edge
     node_index traverse_back(node_index node, char prev_char) const;
 
-    // Given a node index, call the target nodes of all edges outgoing from it.
-    void adjacent_outgoing_nodes(node_index node,
-                                 const std::function<void(node_index)> &callback) const;
-    // Given a node index, call the source nodes of all edges incoming to it.
-    void adjacent_incoming_nodes(node_index node,
-                                 const std::function<void(node_index)> &callback) const;
-
     size_t outdegree(node_index) const;
     bool has_single_outgoing(node_index) const;
     bool has_multiple_outgoing(node_index) const;