Theseus aligner

.gitlab-ci.yml

@@ -158,7 +158,7 @@ oldest-supported-compiler-job:
     GIT_SUBMODULE_STRATEGY: none
     # DO NOT change this version number without updating the README to reflect
     # the requirement bump.
-    COMPILER_VERSION: 9
+    COMPILER_VERSION: 10
 
 
 # We define one job to do the Docker container build

.gitmodules

@@ -133,3 +133,7 @@
 [submodule "deps/mimalloc"]
 	path = deps/mimalloc
 	url = https://github.com/microsoft/mimalloc.git
+[submodule "deps/theseus-lib"]
+	path = deps/theseus-lib
+	url = https://github.com/albertjimenezbl/theseus-lib.git
+	branch = vg-integration

BOTS.md

@@ -0,0 +1,75 @@
+# VG Project Notes
+
+## Building
+- New `.cpp` files auto-discovered
+- Build with `make -j8` or `make obj/whatever.o` to build just one .o.
+- You may be getting errors from `clangd`. If these errors seem spurious, stop and demand a `clangd` that works properly.
+
+## Testing
+
+### Running Bash-TAP Tests
+Use `prove -v` (not `bash`) to execute Bash-TAP tests. This provides proper test harness output and better error reporting.
+
+**Important**: Run `prove` from the `test/` directory:
+```bash
+cd test
+prove -v t/26_deconstruct.t
+```
+
+### Running Unit Tests
+To run all unit tests:
+```bash
+./bin/vg test
+```
+- `./bin/vg test "[tag]"` runs tests matching a tag
+
+#### Writing Unit Tests
+- Framework: Catch v2 (header-only)
+- Include: `#include "catch.hpp"` (in `src/unittest/catch.hpp`)
+- Macros: `TEST_CASE("name", "[tags]")`, `SECTION("name")`, `REQUIRE(cond)`
+- Namespace: `vg::unittest`
+- Directory: `src/unittest/`
+
+### Running All Tests
+```bash
+make test
+```
+
+## Writing Code
+
+### HandleGraph API
+The interfaces in libhandlegraph model a bidirected sequence graph (where nodes have DNA sequences and edges can connect to either the start or end of each involved node).
+
+#### Core types
+- `handle_t` - opaque 64-bit value
+- `nid_t` - node ID type
+- `edge_t` = `pair<handle_t, handle_t>`
+
+#### Key HandleGraph methods
+- `get_handle(nid_t, bool is_reverse=false)` → `handle_t`
+- `get_id(handle_t)` → `nid_t`
+- `get_is_reverse(handle_t)` → `bool`
+- `flip(handle_t)` → `handle_t` (toggle orientation)
+- `get_sequence(handle_t)` → `string` (in handle's orientation)
+- `follow_edges(handle_t, bool go_left, iteratee)` - iterate neighbors
+- `for_each_handle(iteratee, bool parallel=false)` - iterate all nodes
+- `for_each_edge(iteratee, bool parallel=false)` - iterate all edges
+- `has_edge(handle_t left, handle_t right)` → `bool`
+
+#### MutableHandleGraph additions
+- `create_handle(string seq)` / `create_handle(string seq, nid_t id)` → `handle_t`
+- `create_edge(handle_t left, handle_t right)`
+- `destroy_handle(handle_t)` / `destroy_edge(handle_t, handle_t)`
+
+#### HandleGraph algorithms
+- Things like `topological_sort.hpp` and copy_graph.hpp` are in `deps/libhandlegraph/src/include/handlegraph/algorithms`.
+
+#### bdsg::HashGraph
+- Header: `deps/libbdsg/bdsg/include/bdsg/hash_graph.hpp`
+- Implements MutablePathMutableHandleGraph
+- Go-to handlegraph implementation to use
+- In libbdsg
+
+### Utilities
+- `reverse_complement(string)` → `string` in src/utility.hpp
+

CLAUDE.md

@@ -0,0 +1 @@
+BOTS.md

Makefile

@@ -74,7 +74,7 @@ INCLUDE_FLAGS :=-I$(CWD)/$(INC_DIR) -I. -I$(CWD)/$(SRC_DIR) -I$(CWD)/$(UNITTEST_
 # These need to come before library search paths from LDFLAGS or we won't
 # prefer linking vg-installed dependencies over system ones.
 LD_LIB_DIR_FLAGS := -L$(CWD)/$(LIB_DIR)
-LD_LIB_FLAGS := -lvcflib -lwfa2 -ltabixpp -lgssw -lssw -lsublinearLS -lpthread -lncurses -lgcsa2 -lgbwtgraph -lgbwt -lkff -ldivsufsort -ldivsufsort64 -lvcfh -lraptor2 -lpinchesandcacti -l3edgeconnected -lsonlib -lfml -lstructures -lbdsg -lxg -lsdsl -lzstd -lhandlegraph -lcrypto
+LD_LIB_FLAGS := -lvcflib -lwfa2 -ltabixpp -lgssw -lssw -lsublinearLS -lpthread -lncurses -lgcsa2 -lgbwtgraph -lgbwt -lkff -ldivsufsort -ldivsufsort64 -lvcfh -lraptor2 -lpinchesandcacti -l3edgeconnected -lsonlib -lfml -lstructures -lbdsg -lxg -lsdsl -lzstd -ltheseus -lhandlegraph -lcrypto
 # We omit Boost Program Options for now; we find it in a platform-dependent way.
 # By default it has no suffix
 BOOST_SUFFIX=""
@@ -104,7 +104,8 @@ ifeq ($(shell uname -s),Darwin)
     LD_UTIL_RPATH_FLAGS=""
 
     # Homebrew installs a Protobuf that uses an Abseil that is built with C++17, so we need to build with at least C++17
-    CXX_STANDARD?=17
+	# C++20 for spaceship operator and ranges
+    CXX_STANDARD?=20
 
     # We may need libraries from Macports
     ifeq ($(shell if [ -d /opt/local/lib ];then echo 1;else echo 0;fi), 1)
@@ -229,8 +230,9 @@ else
     $(info Compiler $(CXX) is assumed to be GCC)
 
 	# gbwtgraph uses inline variables and our oldest supported compiler has
-	# C++17, so we should use C++17
-    CXX_STANDARD?=17
+	# C++17, so we should use at least C++17.
+    # C++20 for spaceship operator and ranges
+    CXX_STANDARD?=20
 
     # Set an rpath for vg and dependency utils to find installed libraries
     LD_UTIL_RPATH_FLAGS="-Wl,-rpath,$(CWD)/$(LIB_DIR)"
@@ -372,6 +374,8 @@ IPS4O_DIR=deps/ips4o
 BBHASH_DIR=deps/BBHash
 MIO_DIR=deps/mio
 ATOMIC_QUEUE_DIR=deps/atomic_queue
+THESEUS_DIR:=deps/theseus-lib
+INCLUDE_FLAGS += -I$(CWD)/$(THESEUS_DIR)/include
 
 # Dependencies that go into libvg's archive
 # These go in libvg but come from dependencies
@@ -408,6 +412,7 @@ LIB_DEPS += $(LIB_DIR)/libvgio.a
 LIB_DEPS += $(LIB_DIR)/libhandlegraph.a
 LIB_DEPS += $(LIB_DIR)/libbdsg.a
 LIB_DEPS += $(LIB_DIR)/libxg.a
+LIB_DEPS += $(LIB_DIR)/libtheseus.a
 ifneq ($(shell uname -s),Darwin)
     # On non-Mac (i.e. Linux), where ELF binaries are used, pull in libdw which
     # backward-cpp will use.
@@ -820,7 +825,7 @@ $(INC_DIR)/dynamic/dynamic.hpp: $(DYNAMIC_DIR)/include/dynamic/*.hpp $(DYNAMIC_D
 	+mkdir -p $(INC_DIR)/dynamic && cp -r $(CWD)/$(DYNAMIC_DIR)/include/dynamic/* $(INC_DIR)/dynamic/
 
 $(INC_DIR)/sparsehash/sparse_hash_map: $(wildcard $(SPARSEHASH_DIR)/**/*.cc) $(wildcard $(SPARSEHASH_DIR)/**/*.h)
-	+cd $(SPARSEHASH_DIR) && ./autogen.sh && LDFLAGS="$(LD_LIB_DIR_FLAGS) $(LDFLAGS)" ./configure --prefix=$(CWD) $(FILTER) && $(MAKE) $(FILTER) && $(MAKE) install
+	+cd $(SPARSEHASH_DIR) && ./autogen.sh && LDFLAGS="$(LD_LIB_DIR_FLAGS) $(LDFLAGS)" ./configure --prefix=$(CWD) $(FILTER) && $(MAKE) src/sparsehash/internal/sparseconfig.h $(FILTER) && $(MAKE) install-data $(FILTER)
 
 $(INC_DIR)/sparsepp/spp.h: $(wildcard $(SPARSEPP_DIR)/sparsepp/*.h)
 	+cp -r $(SPARSEPP_DIR)/sparsepp $(INC_DIR)/
@@ -934,6 +939,9 @@ $(LIB_DIR)/libxg.a: $(XG_DIR)/src/*.hpp $(XG_DIR)/src/*.cpp $(INC_DIR)/mmmultima
 	+$(CXX) $(INCLUDE_FLAGS) $(CXXFLAGS) $(CPPFLAGS) -fPIC -DNO_GFAKLUGE -c -o $(XG_DIR)/xg.o $(XG_DIR)/src/xg.cpp $(FILTER)
 	+ar rs $@ $(XG_DIR)/xg.o
 
+$(LIB_DIR)/libtheseus.a: $(LIB_DIR)/libhandlegraph.a $(wildcard $(THESEUS_DIR)/theseus/*.cpp) $(wildcard $(THESEUS_DIR)/include/theseus/*.h)
+	+cd $(THESEUS_DIR) && rm -Rf build && mkdir build && cd build && cmake -DCMAKE_C_COMPILER="$(CC)" -DCMAKE_CXX_COMPILER="$(CXX)" -DCMAKE_VERBOSE_MAKEFILE=ON -DCMAKE_C_FLAGS="-fPIC $(CFLAGS)" -DCMAKE_CXX_FLAGS="-fPIC $(CPPFLAGS)" -DCMAKE_BUILD_TYPE=Release -DENABLE_TESTS=OFF -DFETCHCONTENT_SOURCE_DIR_LIBHANDLEGRAPH=$(CWD)/$(LIBHANDLEGRAPH_DIR) .. $(FILTER) && $(MAKE) theseus $(FILTER) && cp libtheseus.a $(CWD)/$(LIB_DIR)/
+
 # Auto-git-versioning
 
 # Can be overridden from the environment to supply a version if none is on disk.
@@ -1155,6 +1163,7 @@ clean: clean-vcflib
 	cd $(DEP_DIR) && cd sublinear-Li-Stephens && $(MAKE) clean
 	cd $(DEP_DIR) && cd libhandlegraph && rm -Rf build CMakeCache.txt CMakeFiles
 	cd $(DEP_DIR) && cd libvgio && rm -Rf build CMakeCache.txt CMakeFiles
+	cd $(DEP_DIR) && cd theseus-lib && rm -Rf build
 	cd $(DEP_DIR) && cd raptor && cd build && find . -not \( -name '.gitignore' -or -name 'pkg.m4' \) -delete
 	# lru_cache is never built because it is header-only
 	# bash-tap is never built either

README.md

@@ -93,7 +93,7 @@ On other distros, or if you do not have root access, you will need to perform th
                          liblzma-dev liblz4-dev libffi-dev libcairo-dev libboost-all-dev \
                          libzstd-dev pybind11-dev python3-pybind11 libssl-dev kmc
 
-At present, you will need GCC version 9 or greater, with support for C++17, to compile vg. (Check your version with `gcc --version`.) GCC up to 11.4.0 is supported.
+At present, you will need GCC version 10 or greater, with support for C++20, to compile vg. (Check your version with `gcc --version`.) GCC up to 11.4.0 is supported.
 
 Other libraries may be required. Please report any build difficulties.
 

src/cactus.cpp

@@ -999,8 +999,8 @@ VG cactus_to_vg(stCactusGraph* cactus_graph) {
     return vg_graph;
 }
 
-VG cactusify(VG& graph) {
-    if (graph.size() == 0) {
+VG cactusify(const PathHandleGraph& graph) {
+    if (graph.get_node_count() == 0) {
         return VG();
     }
     auto parts = handle_graph_to_cactus(graph, unordered_set<string>());

src/cactus.hpp

@@ -46,7 +46,7 @@ VG cactus_to_vg(stCactusGraph* cactus_graph);
 
 // Convert vg into vg formatted cactus representation
 // Input graph must be sorted!
-VG cactusify(VG& graph);
+VG cactusify(const PathHandleGraph& graph);
 
 }
 

src/cluster.hpp

@@ -212,8 +212,8 @@ class MEMClusterer {
 
 protected:
 
-    class HitNode;
     class HitEdge;
+    class HitNode;
     class HitGraph;
     class DPScoreComparator;
 
@@ -232,7 +232,47 @@ class MEMClusterer {
     /// is closest to the optimal separation
     void deduplicate_cluster_pairs(vector<pair<pair<size_t, size_t>, int64_t>>& cluster_pairs, int64_t optimal_separation);
 };
+
+class MEMClusterer::HitEdge {
+public:
+    HitEdge(size_t to_idx, int32_t weight, int64_t distance) : to_idx(to_idx), weight(weight), distance(distance) {}
+    HitEdge() = default;
+    ~HitEdge() = default;
+
+    /// Index of the node that the edge points to
+    size_t to_idx;
 
+    /// Weight for dynamic programming
+    int32_t weight;
+
+    /// Estimated distance
+    int64_t distance;
+};
+
+class MEMClusterer::HitNode {
+public:
+    HitNode(const MaximalExactMatch& mem, pos_t start_pos, int32_t score) : mem(&mem), start_pos(start_pos), score(score) { }
+    HitNode() = default;
+    ~HitNode() = default;
+
+    const MaximalExactMatch* mem;
+
+    /// Position of GCSA hit in the graph
+    pos_t start_pos;
+
+    /// Score of the exact match this node represents
+    int32_t score;
+
+    /// Score used in dynamic programming
+    int32_t dp_score;
+
+    /// Edges from this node that are colinear with the read
+    vector<HitEdge> edges_from;
+
+    /// Edges to this node that are colinear with the read
+    vector<HitEdge> edges_to;
+};
+
 class MEMClusterer::HitGraph {
 public:
 
@@ -286,46 +326,6 @@ class MEMClusterer::HitGraph {
     UnionFind components;
 };
 
-class MEMClusterer::HitNode {
-public:
-    HitNode(const MaximalExactMatch& mem, pos_t start_pos, int32_t score) : mem(&mem), start_pos(start_pos), score(score) { }
-    HitNode() = default;
-    ~HitNode() = default;
-
-    const MaximalExactMatch* mem;
-
-    /// Position of GCSA hit in the graph
-    pos_t start_pos;
-
-    /// Score of the exact match this node represents
-    int32_t score;
-
-    /// Score used in dynamic programming
-    int32_t dp_score;
-
-    /// Edges from this node that are colinear with the read
-    vector<HitEdge> edges_from;
-
-    /// Edges to this node that are colinear with the read
-    vector<HitEdge> edges_to;
-};
-
-class MEMClusterer::HitEdge {
-public:
-    HitEdge(size_t to_idx, int32_t weight, int64_t distance) : to_idx(to_idx), weight(weight), distance(distance) {}
-    HitEdge() = default;
-    ~HitEdge() = default;
-
-    /// Index of the node that the edge points to
-    size_t to_idx;
-
-    /// Weight for dynamic programming
-    int32_t weight;
-
-    /// Estimated distance
-    int64_t distance;
-};
-
 struct MEMClusterer::DPScoreComparator {
 private:
     const vector<HitNode>& nodes;

src/gbwtgraph_helper.cpp

@@ -431,7 +431,7 @@
 void cache_payloads(
     const gbwtgraph::GBZ& gbz,
     const SnarlDistanceIndex& distance_index,
-    hash_map<nid_t, payload_t>& node_id_to_payload,
+    vg::hash_map<nid_t, payload_t>& node_id_to_payload,
     ZipCodeCollection* oversized_zipcodes,
     bool progress
 ) {
@@ -442,22 +442,37 @@
 
     const handlegraph::HandleGraph* graph_ptr = (const handlegraph::HandleGraph*) &gbz.graph;
 
+    double total_zipcode_time = 0.0, total_decoder_time = 0.0;
+    std::atomic<uint64_t> node_count = 0;
     gbz.graph.for_each_handle([&](const handle_t& handle) {
         nid_t node_id = gbz.graph.get_id(handle);
-        ZipCode zipcode;
         pos_t pos = make_pos_t(node_id, false, 0);
-        zipcode.fill_in_zipcode_from_pos(distance_index, pos, true, graph_ptr);
+        ZipCode zipcode;
+        zipcode.fill_in_zipcode_from_pos(distance_index, pos, false, graph_ptr);
+        zipcode.fill_in_full_decoder();
+        if (++node_count % 10000 == 0 && progress) {
+            double telapsed = gbwt::readTimer() - start;
+            #pragma omp critical (cerr)
+            std::cerr << "  Cached " << node_count << " nodes in " << telapsed << "s" << std::endl;
+        }
+
         payload_t payload = zipcode.get_payload_from_zip();
         if (payload == MIPayload::NO_CODE && oversized_zipcodes != nullptr) {
             // The zipcode is too large for the payload field.
             // Add it to the oversized zipcode list.
-            zipcode.fill_in_full_decoder();
-            size_t offset = oversized_zipcodes->size();
-            oversized_zipcodes->emplace_back(zipcode);
+            size_t offset;
+            #pragma omp critical (cache_payloads_zipcodes)
+            {
+                offset = oversized_zipcodes->size();
+                oversized_zipcodes->emplace_back(zipcode);
+            }
             payload = { 0, offset };
         }
-        node_id_to_payload.emplace(node_id, payload);
-    });
+        #pragma omp critical (cache_payloads_map)
+        {
+            node_id_to_payload.emplace(node_id, payload);
+        }
+    }, true);
 
     if (progress) {
         double seconds = gbwt::readTimer() - start;
@@ -519,8 +534,18 @@
     } else {
         // Cache payloads before building the index.
         // A zipcode only depends on the node id.
-        hash_map<nid_t, payload_t> node_id_to_payload;
+        vg::hash_map<nid_t, payload_t> node_id_to_payload;
         node_id_to_payload.reserve(gbz.graph.max_node_id() - gbz.graph.min_node_id());
+        // Re-preload the distance index right before use. find_frequent_kmers
+        // runs for a long time and may evict the mmap'd index pages from the OS
+        // page cache. We also preload eagerly right after loading the index (in
+        // minimizer_main.cpp) so the kernel treats those pages as recently-used;
+        // together the two preloads prevent cache_payloads from page-faulting on
+        // every node under the memory pressure of 32 parallel threads.
+        if (params.progress) {
+            std::cerr << "Preloading distance index";
+        }
+        distance_index->preload(true);
         cache_payloads(gbz, *distance_index, node_id_to_payload, oversized_zipcodes, params.progress);
 
         auto get_payload = [&](const pos_t& pos) -> const code_type* {