diff --git a/.gitignore b/.gitignore index d4ff01ddd3..6081fff992 100644 --- a/.gitignore +++ b/.gitignore @@ -21,6 +21,8 @@ database-cli *._* metagraph/build* metagraph/**/build* +metagraph/api/python/pymetagraph/pymetagraph/_pymetagraph_core* +*.egg-info .*.swp .ipynb_checkpoints *lsf.o* diff --git a/metagraph/CMakeLists.txt b/metagraph/CMakeLists.txt index ebb4c74487..87cbaf5962 100644 --- a/metagraph/CMakeLists.txt +++ b/metagraph/CMakeLists.txt @@ -20,6 +20,18 @@ include(CheckFunctionExists) include(CheckLibraryExists) include(CMakeListsHelpers.txt) +# Building the pymetagraph Python extension requires every static library in the +# link chain (including vendored deps and HTSLIB) to be position-independent. +# This is opt-in so the default static `metagraph` binary build is unaffected. +option(BUILD_PYTHON_BINDINGS "Build the pymetagraph pybind11 extension" OFF) +if(BUILD_PYTHON_BINDINGS) + set(CMAKE_POSITION_INDEPENDENT_CODE ON) + + # jemalloc can exhaust static TLS when the extension is imported. + # Disable it for this build directory; ordinary builds retain their allocator. + set(CMAKE_DISABLE_FIND_PACKAGE_Jemalloc ON) +endif() + enable_testing() include(GoogleTest) @@ -202,6 +214,7 @@ configure_file(CMakeListsKMC.txt.in KMC/CMakeLists.txt @ONLY) execute_process(COMMAND ${CMAKE_COMMAND} -G "${CMAKE_GENERATOR}" -DCMAKE_CXX_COMPILER_LAUNCHER=${CMAKE_CXX_COMPILER_LAUNCHER} -DCMAKE_C_COMPILER_LAUNCHER=${CMAKE_C_COMPILER_LAUNCHER} + -DCMAKE_POSITION_INDEPENDENT_CODE=${CMAKE_POSITION_INDEPENDENT_CODE} . RESULT_VARIABLE result WORKING_DIRECTORY ${PROJECT_BINARY_DIR}/KMC @@ -313,7 +326,12 @@ if(APPLE) get_filename_component(ICU_DATA_LIBRARY_PATH ${ICU_DATA_LIBRARY} DIRECTORY) endif() -set(Boost_USE_STATIC_LIBS ON) +# Distribution static Boost libraries may lack PIC. +if(BUILD_PYTHON_BINDINGS AND CMAKE_SYSTEM_NAME STREQUAL "Linux") + set(Boost_USE_STATIC_LIBS OFF) +else() + set(Boost_USE_STATIC_LIBS ON) +endif() find_package(Boost REQUIRED COMPONENTS iostreams) find_library(LIBDEFLATE_LIBRARY NAMES deflate HINTS ENV LIBRARY_PATH PATHS ${HOMEBREW_DIR}/lib REQUIRED) @@ -386,17 +404,29 @@ else() endif() include(ExternalProject) +# HTSLIB is built via autoconf/make, outside CMAKE_POSITION_INDEPENDENT_CODE's reach, +# so -fPIC has to be injected into its CFLAGS/CPPFLAGS by hand when building the +# Python extension (see BUILD_PYTHON_BINDINGS above). +set(HTSLIB_PIC_FLAG "") +if(BUILD_PYTHON_BINDINGS) + set(HTSLIB_PIC_FLAG "-fPIC ") +endif() +if(CMAKE_GENERATOR MATCHES "Makefiles") + set(HTSLIB_MAKE "$(MAKE)") +else() + find_program(HTSLIB_MAKE NAMES gmake make REQUIRED) +endif() ExternalProject_Add(HTSLIB SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/external-libraries/htslib PREFIX "${PROJECT_BINARY_DIR}/HTSLIB-build" CONFIGURE_COMMAND ${CMAKE_COMMAND} -E env "CC=${HTSLIB_CC}" - "CPPFLAGS=${OSX_ARCH_FLAGS}\ -I${CMAKE_CURRENT_SOURCE_DIR}/external-libraries/zlib" - "CFLAGS=${OSX_ARCH_FLAGS}\ -I${CMAKE_CURRENT_SOURCE_DIR}/external-libraries/zlib" + "CPPFLAGS=${HTSLIB_PIC_FLAG}${OSX_ARCH_FLAGS}\ -I${CMAKE_CURRENT_SOURCE_DIR}/external-libraries/zlib" + "CFLAGS=${HTSLIB_PIC_FLAG}${OSX_ARCH_FLAGS}\ -I${CMAKE_CURRENT_SOURCE_DIR}/external-libraries/zlib" "LDFLAGS=${OSX_ARCH_FLAGS}\ -L${PROJECT_BINARY_DIR}/external-libraries/zlib" /bin/sh -c "autoreconf -i ${CMAKE_CURRENT_SOURCE_DIR}/external-libraries/htslib && exec ${CMAKE_CURRENT_SOURCE_DIR}/external-libraries/htslib/configure --disable-libcurl --disable-bz2 --disable-lzma" - BUILD_COMMAND $(MAKE) -j${N} lib-static + BUILD_COMMAND ${HTSLIB_MAKE} -j${N} lib-static #INSTALL_COMMAND $(MAKE) install prefix=${HTSLIB_INSTALL} # this installs both static and dynamic # Install only static lib to avoid the dependency INSTALL_COMMAND @@ -434,6 +464,8 @@ list(REVERSE metagraph_core_files) # for compiling faster list(FILTER metagraph_core_files EXCLUDE REGEX ".*\\._.*") list(FILTER metagraph_core_files EXCLUDE REGEX ".*main.cpp") list(FILTER metagraph_core_files EXCLUDE REGEX ".*/cli/.*") +list(FILTER metagraph_core_files EXCLUDE REGEX ".*/pybindings/.*") + add_library(metagraph-core STATIC ${metagraph_core_files}) # Set up include directories for metagraph-core @@ -577,6 +609,46 @@ if(METAGRAPH_BUILD_EXECUTABLE) add_dependencies(metagraph link_target) endif() +#------------------- +# pymetagraph Python extension +#------------------- +if(BUILD_PYTHON_BINDINGS) + find_package(Python 3.9 REQUIRED COMPONENTS Interpreter Development.Module) + execute_process( + COMMAND ${Python_EXECUTABLE} -m pybind11 --cmakedir + OUTPUT_VARIABLE pybind11_DIR + OUTPUT_STRIP_TRAILING_WHITESPACE + COMMAND_ERROR_IS_FATAL ANY + ) + find_package(pybind11 CONFIG REQUIRED) + + pybind11_add_module(_pymetagraph_core src/pybindings/pymetagraph_bindings.cpp) + target_link_libraries(_pymetagraph_core PRIVATE metagraph-cli metagraph-core) + add_dependencies(_pymetagraph_core HTSLIB) + install(TARGETS _pymetagraph_core DESTINATION pymetagraph COMPONENT _pymetagraph_core) + if(CMAKE_SYSTEM_NAME STREQUAL "Linux") + # Use an inherited RPATH so Boost's own shared Boost dependencies also + # resolve from the wheel, without relying on a matching system Boost. + target_link_options(_pymetagraph_core PRIVATE "LINKER:--disable-new-dtags") + set_target_properties(_pymetagraph_core PROPERTIES INSTALL_RPATH "$ORIGIN/.libs") + install(CODE " + set(boost_library \"$\") + file(GET_RUNTIME_DEPENDENCIES + LIBRARIES \"\${boost_library}\" + RESOLVED_DEPENDENCIES_VAR boost_dependencies + UNRESOLVED_DEPENDENCIES_VAR missing_boost_dependencies + PRE_INCLUDE_REGEXES \"^libboost_\" + PRE_EXCLUDE_REGEXES \".*\") + if(missing_boost_dependencies) + message(FATAL_ERROR \"Unresolved Boost runtime libraries: \${missing_boost_dependencies}\") + endif() + file(INSTALL \"\${boost_library}\" \${boost_dependencies} + DESTINATION \"\$ENV{DESTDIR}\${CMAKE_INSTALL_PREFIX}/pymetagraph/.libs\" + TYPE SHARED_LIBRARY FOLLOW_SYMLINK_CHAIN) + " COMPONENT _pymetagraph_core) + endif() +endif() + #------------------- # install library diff --git a/metagraph/api/python/pymetagraph/README.md b/metagraph/api/python/pymetagraph/README.md new file mode 100644 index 0000000000..2efb524d86 --- /dev/null +++ b/metagraph/api/python/pymetagraph/README.md @@ -0,0 +1,121 @@ +# pymetagraph + +A [Readfish](https://github.com/LooseLab/readfish) `Aligner` plugin backed by a prebuilt +metagraph graph + annotation index. By default (`method="query"`) it reports label-presence hits +(which reference sequence(s) a read's k-mers match, above a discovery/presence-fraction threshold) +rather than real base-level alignments, so `r_st`/`r_en`/`strand` on the returned alignment +objects are placeholder values. Set `method="align"` for real seed-and-extend alignment with +actual coordinates and orientation instead — see [Usage](#usage) below. + +## Building + +Unlike a typical small pybind11 project, this is not an instant pip install: metagraph needs its +git submodules checked out and takes a while to compile from scratch. But it *is* a plain +`pip install`-driven build — pip's normal build isolation fetches `scikit-build-core` and +`pybind11` (this project's `[build-system].requires`) into a throwaway env on its own, so you +don't need to pre-install anything Python-side yourself. Recommended workflow: + +1. From the repo root, make sure submodules are initialized (see the top-level README). +2. `sdsl-lite` is pre-built separately (via its own `install.sh`, not by the main CMake project) + and the resulting `external-libraries/sdsl-lite/lib/libsdsl.a` is **not** position-independent + by default, which breaks the final link into a `.so` even with `BUILD_PYTHON_BINDINGS=ON`. It + only needs to be rebuilt once (PIC static objects link fine into the normal static `metagraph` + binary too, so this doesn't need to be undone for the default build): + ``` + cd metagraph/external-libraries/sdsl-lite/build + ./clean.sh + cmake -DCMAKE_POSITION_INDEPENDENT_CODE=ON \ + -DCMAKE_INSTALL_PREFIX=$(cd .. && pwd) .. + make -j$(nproc) && make install + ``` +3. Build and install — this is the only step most users need, and it never requires running + `cmake` yourself: + ``` + cd metagraph/api/python/pymetagraph + pip install -e . + ``` + Under the hood, pip's `scikit-build-core` build backend invokes CMake for you (with + `-DBUILD_PYTHON_BINDINGS=ON` already wired in via this project's `pyproject.toml`) inside a + throwaway build environment, so you don't need CMake installed in your active environment or + any familiarity with its flags. + + If you're in a conda environment that ships its own (shared-only) Boost, CMake may pick that up + ahead of the system Boost metagraph needs (static libs), failing with `boost_iostreams_FOUND to + FALSE`. Since pip's build backend is the one invoking CMake here, pass the override through the + standard `CMAKE_ARGS` environment variable, which `scikit-build-core` forwards to it: + ``` + CMAKE_ARGS="-DBoost_DIR=/usr/lib/x86_64-linux-gnu/cmake/Boost-1.83.0" pip install -e . + ``` + (adjust the path/version to whatever `find /usr/lib -name 'BoostConfig.cmake'` reports on your + system). + +### Faster iteration on the C++ side (advanced, driving CMake directly) + +The steps above are a *different, independent path* from what follows: `pip install -e .` never +invokes `cmake` on your command line, it happens inside pip's own build backend. The workflow +below is only for developers iterating on the C++ binding itself, who want incremental rebuilds +without re-running pip's build backend (and its temp-env setup) on every change. It talks to CMake +directly, so **it needs `-D` flags spelled out explicitly** — none of `pyproject.toml`'s settings +or pip's `CMAKE_ARGS` env var apply here, since pip isn't involved: +``` +cmake -S metagraph -B metagraph/api/python/pymetagraph/build -DBUILD_PYTHON_BINDINGS=ON +# add -DBoost_DIR=... too if you hit the same conda/Boost conflict as above +cmake --build metagraph/api/python/pymetagraph/build --target _pymetagraph_core -j$(nproc) +cd metagraph/api/python/pymetagraph +pip install -e . --no-build-isolation # requires scikit-build-core + pybind11 pre-installed +``` + +### Note: jemalloc is disabled for this extension + +`BUILD_PYTHON_BINDINGS=ON` also sets `CMAKE_DISABLE_FIND_PACKAGE_Jemalloc=ON`, so +`_pymetagraph_core` is built without jemalloc even if it's installed on the system. This is +deliberate: dlopen()-ing a library linked against jemalloc can fail at import time with `cannot +allocate memory in static TLS block` (a glibc static-TLS-space limitation), which would otherwise +require every user of this plugin to remember `LD_PRELOAD=libjemalloc.so.2` before launching +Readfish. Since this plugin is aimed at non-power-users, it trades away folly's jemalloc-backed +allocator optimizations for an extension that just imports cleanly with no extra steps. This only +affects `_pymetagraph_core` — the regular static `metagraph` CLI binary still picks up jemalloc +normally when built without `BUILD_PYTHON_BINDINGS`. + +## Usage + +Two unrelated query methods are available via the `method` setting: + +- `"query"` (default): label-presence hits via `AnnotatedDBG::get_top_labels` — the same + primitive backing the CLI's `query` command. Fast and tolerant of sequencing errors, but + `r_st`/`r_en`/`strand` are dummy values (`0`, `len(seq)`, `+1`); only `ctg` (the matched label) + is meaningful. + ```toml + [mapper_settings.pymetagraph] + input = "/path/to/graph.dbg" + annotator = "/path/to/graph.column.annodbg" + threads = 4 + method = "query" + num_top_labels = 1 + discovery_fraction = .4 + presence_fraction = .1 + ``` +- `"align"`: real seed-and-extend alignment via metagraph's alignment_redone module. Neither the + `metagraph align` CLI command nor the server's `/align` endpoint uses this engine. Slower, but `r_st`/`r_en`/`strand` reflect the + actual matched region and orientation of the read. + ```toml + [mapper_settings.pymetagraph] + input = "/path/to/graph.dbg" + annotator = "/path/to/graph.column.annodbg" + threads = 4 + method = "align" + seed_length = 19 + max_num_nodes_per_seq_char = 10.0 + min_exact_match = .7 + connect_anchors = true + extend_chains = true + ``` + +## Known limitations + +- `r_st`/`r_en`/`strand` are always dummy values (`0`, `len(seq)`, `+1`) under `method="query"`; + only `ctg` (the matched label) is meaningful there. `method="align"` reports real values. +- Exactly one annotator file is supported. +- Malformed configuration may, in rare cases, abort the process rather than raise a Python + exception, since metagraph's underlying CLI config parser calls `exit()` on some invalid inputs. + `Aligner.validate()` checks the common cases (missing files) up front. diff --git a/metagraph/api/python/pymetagraph/pymetagraph/__init__.py b/metagraph/api/python/pymetagraph/pymetagraph/__init__.py new file mode 100644 index 0000000000..39515cfb2b --- /dev/null +++ b/metagraph/api/python/pymetagraph/pymetagraph/__init__.py @@ -0,0 +1,7 @@ +# -*- coding: utf-8 -*- + +"""Readfish-compatible Aligner backed by a metagraph index.""" + +from .aligner import Aligner + +__all__ = ["Aligner"] diff --git a/metagraph/api/python/pymetagraph/pymetagraph/aligner.py b/metagraph/api/python/pymetagraph/pymetagraph/aligner.py new file mode 100644 index 0000000000..656e46752c --- /dev/null +++ b/metagraph/api/python/pymetagraph/pymetagraph/aligner.py @@ -0,0 +1,127 @@ +from __future__ import annotations + +import sys +from pathlib import Path +from typing import Dict, Iterable, List, Optional + +from readfish._config import Barcode, Region +from readfish.plugins.abc import AlignerABC +from readfish.plugins.utils import Result + +from . import _pymetagraph_core + + +class Aligner(AlignerABC): + """Readfish Aligner plugin backed by a metagraph graph + annotation index. + + Two unrelated query methods are available, picked via ``method``: + + - ``"query"`` (default): label-presence hits via + ``AnnotatedDBG::get_top_labels``. ``Alignment.r_st``/``r_en``/``strand`` + are placeholder values (``0``, ``len(seq)``, ``+1``) since only the + matched label is meaningful for accept/reject decisions in this mode. + - ``"align"``: real seed-and-extend alignment via metagraph's + alignment_redone module. Neither the ``metagraph align`` CLI command + nor the server's ``/align`` endpoint uses this engine. + ``r_st``/``r_en``/``strand`` reflect + the actual matched region and orientation of the read. + + :param debug_log: Filename (or ``"stdout"``/``"stderr"``) for debug logging. + :keyword input: Path to the metagraph graph file (``.dbg``). + :keyword annotator: Path to the metagraph annotation file. Required for + ``method="query"``; optional for ``method="align"`` (omitting it runs + plain graph alignment with no label lookups). + :keyword threads: Number of threads to use for querying (default: 1). + :keyword method: Either ``"query"`` (default) or ``"align"``. + :keyword num_top_labels: Max number of labels to consider per query. + Only used by ``method="query"``. + :keyword discovery_fraction: Min fraction of k-mers required to be present + in a label for it to count as a hit (default: 0.7). Only used by + ``method="query"``. + :keyword presence_fraction: Min fraction of k-mers required to be present + in the graph at all before querying labels (default: 0.0). Only used + by ``method="query"``. + :keyword seed_length: Minimum seed length. Only used by ``method="align"``. + :keyword max_num_nodes_per_seq_char: Max graph nodes to explore per query + character in each extension or gap fill (default: 5.0). Only used by + ``method="align"``. + :keyword max_num_seeds_per_locus: Drop the seeds shorter than k at a + query position that has more than this many (default: 1000). Only + used by ``method="align"``. + :keyword min_exact_match: Min fraction of nucleotides covered by seeds + required to align (default: 0.7). Only used by ``method="align"``. + :keyword connect_anchors: If ``True`` (default), traverse the graph to + align query regions falling between anchors. Only used by + ``method="align"``. + :keyword extend_chains: If ``True`` (default), perform ends-free + extension from the first/last anchors in a chain. Only used by + ``method="align"``. + """ + + def __init__(self, debug_log: Optional[str] = None, **kwargs): + if debug_log: + if debug_log == "stdout": + self.logfile = sys.stdout + elif debug_log == "stderr": + self.logfile = sys.stderr + else: + self.logfile = open(debug_log, "w") + else: + self.logfile = None + + if "input" not in kwargs: + raise AttributeError('Required argument "input" not found.') + + self.kwargs = kwargs + self.validate() + self.index = _pymetagraph_core.Index(**kwargs) + + def validate(self) -> None: + index_path = Path(self.kwargs["input"]) + if not index_path.is_file(): + raise FileNotFoundError(f"{index_path} does not exist") + + # "query" (the default) needs an annotation to have any labels to + # report; "align" can run against the bare graph (see _pymetagraph_core.Index). + method = self.kwargs.get("method", "query") + annotator = self.kwargs.get("annotator") + if not annotator: + if method != "align": + raise AttributeError('Required argument "annotator" not found.') + return + if not Path(annotator).is_file(): + raise FileNotFoundError(f"{annotator} does not exist") + + @property + def initialised(self) -> bool: + return True + + def describe(self, regions: List[Region], barcodes: Dict[str, Barcode]) -> str: + return ( + f"Using the pymetagraph plugin. Graph: {self.kwargs['input']}, " + f"annotator: {self.kwargs.get('annotator', 'none')}." + ) + + def map_reads(self, basecall_results: Iterable[Result]) -> Iterable[Result]: + skipped = [] + kept = [] + for result in basecall_results: + if result.seq: + kept.append(result) + else: + skipped.append(result) + + if kept: + alignments = self.index.query_batch([r.seq for r in kept]) + for result, alignment in zip(kept, alignments): + result.alignment_data = [] if alignment.ctg == "*" else [alignment] + yield result + + for result in skipped: + result.alignment_data = [] + yield result + + def disconnect(self) -> None: + self.index = None + if self.logfile and self.logfile not in (sys.stdout, sys.stderr): + self.logfile.close() diff --git a/metagraph/api/python/pymetagraph/pyproject.toml b/metagraph/api/python/pymetagraph/pyproject.toml new file mode 100644 index 0000000000..233d199857 --- /dev/null +++ b/metagraph/api/python/pymetagraph/pyproject.toml @@ -0,0 +1,28 @@ +[build-system] +requires = ["scikit-build-core", "pybind11"] +build-backend = "scikit_build_core.build" + +[project] +name = "pymetagraph" +version = "0.1.2" +description = "Readfish-compatible Aligner plugin backed by a metagraph index" +readme = "README.md" +authors = [ + { name = "ratschlab", email = "grlab@ratschlab.org" }, +] +requires-python = ">=3.9" +dependencies = ["readfish"] +classifiers = [ + "Development Status :: 3 - Alpha", + "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", + "Programming Language :: Python :: 3 :: Only", +] + +[tool.scikit-build] +build-dir = "build" +cmake.source-dir = "../../.." +cmake.args = ["-DBUILD_PYTHON_BINDINGS=ON"] +cmake.build-type = "Release" +wheel.packages = ["pymetagraph"] +build.targets = ["_pymetagraph_core"] +install.components = ["_pymetagraph_core"] diff --git a/metagraph/api/python/pymetagraph/tests/test_aligner.py b/metagraph/api/python/pymetagraph/tests/test_aligner.py new file mode 100644 index 0000000000..f20d5d0a63 --- /dev/null +++ b/metagraph/api/python/pymetagraph/tests/test_aligner.py @@ -0,0 +1,83 @@ +import os + +import pytest +from readfish.plugins.utils import Result + +from pymetagraph import Aligner + +FIXTURE_DIR = os.environ.get("PYMETAGRAPH_TEST_FIXTURE_DIR", "/tmp/pymetagraph_fixture") +GRAPH = os.path.join(FIXTURE_DIR, "idx.dbg") +ANNOTATOR = os.path.join(FIXTURE_DIR, "idx.column.annodbg") + +MATCHING_SEQ = ( + "GTTAACTTGCCGTCAGCCTTTTCTTTGACCTCTTCTTTCTGTTCATGTGTATTTGCTGTC" + "TCTTAGCCCAGACTTCCCGTGTCCTTTCCACCGGGCCTTTGAGAGGTCACAGGGTCTTGA" +) +EXPECTED_LABEL = ( + "ENST00000456328.2|ENSG00000223972.5|OTTHUMG00000000961.2|" + "OTTHUMT00000362751.1|DDX11L1-202|DDX11L1|1657|processed_transcript|" +) + + +@pytest.fixture +def aligner(): + a = Aligner(debug_log=None, input=GRAPH, annotator=ANNOTATOR) + yield a + a.disconnect() + + +@pytest.fixture +def align_aligner(): + a = Aligner(debug_log=None, input=GRAPH, annotator=ANNOTATOR, method="align") + yield a + a.disconnect() + + +def test_validate_missing_input(): + with pytest.raises(AttributeError): + Aligner(debug_log=None) + + +def test_initialised(aligner): + assert aligner.initialised is True + + +def test_describe(aligner): + description = aligner.describe(regions=[], barcodes={}) + assert GRAPH in description + + +def test_map_reads(aligner): + results = [ + Result(channel=1, read_id="matching", seq=MATCHING_SEQ), + Result(channel=2, read_id="no_seq", seq=""), + Result(channel=3, read_id="garbage", seq="N" * 100), + ] + + by_id = {r.read_id: r for r in aligner.map_reads(iter(results))} + + assert len(by_id) == 3 + assert by_id["matching"].alignment_data + assert by_id["matching"].alignment_data[0].ctg == EXPECTED_LABEL + assert by_id["no_seq"].alignment_data == [] + assert by_id["garbage"].alignment_data == [] + + +def test_map_reads_align_method(align_aligner): + results = [ + Result(channel=1, read_id="matching", seq=MATCHING_SEQ), + Result(channel=2, read_id="no_seq", seq=""), + Result(channel=3, read_id="garbage", seq="N" * 100), + ] + + by_id = {r.read_id: r for r in align_aligner.map_reads(iter(results))} + + assert len(by_id) == 3 + matching = by_id["matching"].alignment_data + assert matching + assert matching[0].ctg == EXPECTED_LABEL + # Real coordinates, unlike method="query"'s dummy 0/len(seq)/+1. + assert 0 <= matching[0].r_st < matching[0].r_en <= len(MATCHING_SEQ) + assert matching[0].strand in (1, -1) + assert by_id["no_seq"].alignment_data == [] + assert by_id["garbage"].alignment_data == [] diff --git a/metagraph/src/cli/align.cpp b/metagraph/src/cli/align.cpp index 3e25af0564..af23eb4c19 100644 --- a/metagraph/src/cli/align.cpp +++ b/metagraph/src/cli/align.cpp @@ -9,6 +9,7 @@ #include "graph/representation/canonical_dbg.hpp" #include "graph/alignment/dbg_aligner.hpp" #include "graph/alignment/aligner_labeled.hpp" +#include "graph/alignment_redone/aligner_config.hpp" #include "graph/annotated_dbg.hpp" #include "graph/graph_extensions/node_first_cache.hpp" #include "seq_io/sequence_io.hpp" @@ -68,6 +69,11 @@ DBGAlignerConfig initialize_aligner_config(const Config &config, return c; } +align_redone::DBGAlignerConfig initialize_redone_aligner_config(const Config &config, + const DeBruijnGraph &graph) { + return { initialize_aligner_config(config, graph) }; +} + void map_sequences_in_file(const std::string &file, const DeBruijnGraph &graph, const Config &config, diff --git a/metagraph/src/cli/align.hpp b/metagraph/src/cli/align.hpp index 92137022f3..af382cf2b0 100644 --- a/metagraph/src/cli/align.hpp +++ b/metagraph/src/cli/align.hpp @@ -12,6 +12,9 @@ class DeBruijnGraph; namespace align { struct DBGAlignerConfig; } // namespace align +namespace align_redone { +struct DBGAlignerConfig; +} // namespace align_redone } // namespace graph @@ -22,6 +25,11 @@ class Config; graph::align::DBGAlignerConfig initialize_aligner_config(const Config &config, const graph::DeBruijnGraph &graph); +// Options for the new alignment engine (graph/alignment_redone), which is not +// yet used by the align command. +graph::align_redone::DBGAlignerConfig +initialize_redone_aligner_config(const Config &config, const graph::DeBruijnGraph &graph); + int align_to_graph(Config *config); } // namespace cli diff --git a/metagraph/src/cli/build.cpp b/metagraph/src/cli/build.cpp index 8345d43122..8d59bacf47 100644 --- a/metagraph/src/cli/build.cpp +++ b/metagraph/src/cli/build.cpp @@ -248,7 +248,8 @@ int build_graph(Config *config) { } } 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)); + graph.reset(new DBGSSHash(files.at(0), config->k, config->graph_mode, + config->num_chars, config->tmp_dir)); if (files.size() > 1) { logger->error("DBGSSHash does not support multiple input files."); exit(1); diff --git a/metagraph/src/graph/alignment/aligner_config.hpp b/metagraph/src/graph/alignment/aligner_config.hpp index 459e2c1f05..393477e90f 100644 --- a/metagraph/src/graph/alignment/aligner_config.hpp +++ b/metagraph/src/graph/alignment/aligner_config.hpp @@ -53,10 +53,10 @@ struct DBGAlignerConfig { bool no_backtrack = false; bool seed_complexity_filter = true; - bool alignment_edit_distance; - int8_t alignment_match_score; - int8_t alignment_mm_transition_score; - int8_t alignment_mm_transversion_score; + bool alignment_edit_distance = false; + int8_t alignment_match_score = 2; + int8_t alignment_mm_transition_score = 3; + int8_t alignment_mm_transversion_score = 3; ScoreMatrix score_matrix; diff --git a/metagraph/src/graph/alignment/alignment.hpp b/metagraph/src/graph/alignment/alignment.hpp index 8e87cd43a7..dd8158e742 100644 --- a/metagraph/src/graph/alignment/alignment.hpp +++ b/metagraph/src/graph/alignment/alignment.hpp @@ -417,6 +417,11 @@ class AlignmentResults { std::vector alignments_; }; +// Return the string spelled by the path, skipping the first offset characters. +// The path may have disconnects (npos nodes, if it came from a chain alignment). +// Their characters are filled in from overlapping neighbouring nodes, and are +// '$' where no node covers them. If there is an invalid edge, this throws a +// runtime error. std::string spell_path(const DeBruijnGraph &graph, const std::vector &path, size_t offset = 0); diff --git a/metagraph/src/graph/alignment_redone/aligner_config.cpp b/metagraph/src/graph/alignment_redone/aligner_config.cpp new file mode 100644 index 0000000000..82e15894ab --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aligner_config.cpp @@ -0,0 +1,69 @@ +#include "aligner_config.hpp" + +#include "aln_cigar.hpp" + + +namespace mtg::graph::align_redone { + +DBGAlignerConfig::score_t DBGAlignerConfig +::score_cigar(std::string_view reference, + std::string_view query, + const Cigar &cigar) const { + assert(cigar.is_valid(reference, query)); + + if (cigar.empty()) + return 0; + + score_t score = (!cigar.get_clipping() ? left_end_bonus : 0) + + (!cigar.get_end_clipping() ? right_end_bonus : 0); + + auto ref_it = reference.begin(); + auto alt_it = query.begin(); + auto it = cigar.data().begin(); + if (it->first == Cigar::CLIPPED) { + alt_it += it->second; + ++it; + } + + for ( ; it != cigar.data().end(); ++it) { + const auto &op = *it; + switch (op.first) { + case Cigar::CLIPPED: { + alt_it += op.second; + } break; + case Cigar::MATCH: { + score += match_score(std::string_view(ref_it, op.second)); + ref_it += op.second; + alt_it += op.second; + } break; + case Cigar::MISMATCH: { + score += score_sequences(std::string_view(ref_it, op.second), + std::string_view(alt_it, op.second)); + ref_it += op.second; + alt_it += op.second; + } break; + case Cigar::INSERTION: { + score += gap_opening_penalty + (op.second - 1) * gap_extension_penalty; + alt_it += op.second; + } break; + case Cigar::DELETION: { + score += gap_opening_penalty + (op.second - 1) * gap_extension_penalty; + ref_it += op.second; + if (it >= cigar.data().begin() + 2 && (it - 2)->first == Cigar::DELETION + && (it - 1)->first == Cigar::NODE_INSERTION) { + score -= gap_opening_penalty - gap_extension_penalty; + } + } break; + case Cigar::NODE_INSERTION: { + score += gap_opening_penalty + (op.second - 1) * gap_extension_penalty; + } break; + } + } + + assert(ref_it == reference.end()); + assert(alt_it == query.end()); + + return score; +} + +} // namespace mtg::graph::align_redone diff --git a/metagraph/src/graph/alignment_redone/aligner_config.hpp b/metagraph/src/graph/alignment_redone/aligner_config.hpp new file mode 100644 index 0000000000..0e7c911db9 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aligner_config.hpp @@ -0,0 +1,25 @@ +#pragma once + +#include + +#include "graph/alignment/aligner_config.hpp" + + +namespace mtg::graph::align_redone { + +class Cigar; + +// The new engine shares the scoring and search options of the existing aligner +// and adds the parameters of its anchor chaining. +struct DBGAlignerConfig : public align::DBGAlignerConfig { + size_t max_dist_between_seeds = 400; + size_t gap_shrinking_factor = 4; + double chaining_algorithm_switch_cutoff = 1.0; + size_t chaining_bandwidth = 65; + + score_t score_cigar(std::string_view reference, + std::string_view query, + const Cigar &cigar) const; +}; + +} // namespace mtg::graph::align_redone diff --git a/metagraph/src/graph/alignment_redone/aln_chainer.hpp b/metagraph/src/graph/alignment_redone/aln_chainer.hpp new file mode 100644 index 0000000000..e4a3cad40d --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_chainer.hpp @@ -0,0 +1,373 @@ +#pragma once + +#include + +#include "aligner_config.hpp" +#include "aln_query.hpp" +#include "aln_match.hpp" +#include "annotation_buffer.hpp" + +#include "common/logger.hpp" + +namespace mtg::graph::align_redone { + +template +using ChainScores = std::vector>; + +template +using ScoreUpdater = std::function; + +template +using AnchorChain = std::vector>; + +template +using ExtensionStarter = std::function&, // chain + const std::vector& // score traceback + )>; + + +template +using AnchorConnector = std::function::iterator, + const ScoreUpdater& + )>; + +using AlignmentCallback = std::function; + +template +using AnchorExtender = std::function; + +template +struct ChainHash { + inline std::size_t operator()(const AnchorChain &chain) const { + uint64_t hash = 0; + + auto update_hash = [&](auto val) { + hash ^= val + 0x9e3779b9 + (hash << 6) + (hash >> 2); + }; + + for (const auto &[anchor_it, dist] : chain) { + std::for_each(anchor_it->get_path().begin(), anchor_it->get_path().end(), + update_hash); + update_hash(dist); + update_hash(anchor_it->get_coord()); + } + return hash; + } +}; + +template +struct ChainEqual { + inline std::size_t operator()(const AnchorChain &a, + const AnchorChain &b) const { + if (a.size() != b.size()) + return false; + + for (size_t i = 0; i < a.size(); ++i) { + if (a[i].second != b[i].second) + return false; + } + + for (size_t i = 0; i < a.size(); ++i) { + if (a[i].first->get_coord() != b[i].first->get_coord()) + return false; + + const Match &a_i = static_cast(*a[i].first); + const Match &b_i = static_cast(*b[i].first); + if (!(a_i == b_i)) + return false; + } + + return true; + } +}; + +template +bool extend_chain(const AnchorChain &chain, + const AnchorExtender &anchor_extender, + const AlignmentCallback &callback, + const std::function &terminate = []() { return false; }) { + assert(chain.size()); + + std::vector alns; + alns.emplace_back(*chain.back().first); + + for (auto it = chain.rbegin(); it + 1 != chain.rend(); ++it) { + std::vector next_alns; + for (auto&& aln : alns) { + assert(it->first->get_clipping() == aln.get_clipping()); + + size_t dist = it->second; + anchor_extender(it->first, (it + 1)->first, std::move(aln), dist, 0, + [&](Alignment&& next_aln) { next_alns.emplace_back(std::move(next_aln)); } + ); + } + + std::swap(next_alns, alns); + } + + size_t num_alns = 0; + + for (auto&& aln : alns) { + if (terminate()) + return num_alns; + + callback(std::move(aln)); + ++num_alns; + } + + return num_alns; +} + +template +struct AnchorLess { + bool operator()(const Anchor &a, const Anchor &b) const { + return std::make_tuple(a.get_orientation(), a.get_label_class(), b.get_end_clipping(), a.get_clipping()) + < std::make_tuple(b.get_orientation(), b.get_label_class(), a.get_end_clipping(), b.get_clipping()); + } +}; + +template +AnchorIt chain_anchors(const Query &query, + const DBGAlignerConfig &config, + AnchorIt begin, + AnchorIt end, + const AnchorConnector &anchor_connector, + const ExtensionStarter &extension_starter, + const AnchorExtender &anchor_extender, + const AlignmentCallback &callback = [](Alignment&&) {}, + const std::function &terminate = []() { return false; }, + AnnotationBuffer *anno_buffer = nullptr) { + if (terminate() || begin == end) + return end; + + ssize_t query_size = query.get_query().size(); + + assert(std::is_sorted(begin, end, AnchorLess())); + + ChainScores chain_scores; + chain_scores.reserve(end - begin); + std::transform(begin, end, std::back_inserter(chain_scores), [&](const Anchor &a) { + return std::make_tuple(score_match(a, config), + end, 0, a.get_spelling().size(), a.get_seed().size()); + }); + + // forward pass + ssize_t max_gap_between_anchors = std::min( + static_cast(config.max_dist_between_seeds), + query_size + ); + + auto forward_pass = [&](AnchorIt begin, + AnchorIt end, + typename ChainScores::iterator chain_scores) { + if (begin == end) + return; + + auto make_anchor_connector = [&](const AnchorIt i) { + return [&,i](DBGAlignerConfig::score_t score, const AnchorIt last, size_t dist) { + assert(last != i); + auto &[max_score, best_last, best_dist, dist_traversed, queried] = *(chain_scores + (i - begin)); + if (std::tie(score, best_dist) > std::tie(max_score, dist)) { + if (last != end) { + dist_traversed = std::get<3>(*(chain_scores + (last - begin))) + dist; + assert(i->get_seed().end() >= last->get_seed().end()); + queried = std::get<4>(*(chain_scores + (last - begin))) + (i->get_seed().end() - last->get_seed().end()); + } + + max_score = score; + best_last = last; + best_dist = dist; + return true; + } else { + return false; + } + }; + }; + + if (static_cast(end - begin) / query_size + <= config.chaining_algorithm_switch_cutoff) { + // if there are fewer seeds, this algorithm is faster + ssize_t b = max_gap_between_anchors; + ssize_t b_last; + DBGAlignerConfig::score_t best_score = std::get<0>(chain_scores[0]); + size_t traversed = std::get<3>(chain_scores[0]); + size_t queried = begin->get_seed().size(); + ssize_t best_cost = std::numeric_limits::max(); + DBGAlignerConfig::score_t match_score = config.match_score("A"); + do { + AnchorIt j = begin; + for (AnchorIt i = begin; i != end; ++i) { + j = std::find_if(j, end, [&](const Anchor &a) { + return i->get_seed().end() - a.get_seed().end() <= b; + }); + + // align anchors [j, i) to i + anchor_connector(*i, b, j, i, chain_scores + (j - begin), + make_anchor_connector(i)); + + if (std::get<0>(chain_scores[i - begin]) > best_score) { + AnchorIt last; + size_t added_dist; + std::tie(best_score, last, added_dist, traversed, queried) + = chain_scores[i - begin]; + } + } + b_last = b; + b *= config.gap_shrinking_factor; + + // best_score = 1/2(match_score(query_size + match_spelling_size - best_cost)) + // => 2*best_score/match_score - query_size - match_spelling_size = -best_cost + // => best_cost = match_spelling_size + query_size - 2*best_score/match_score + + best_cost = traversed + queried - 2 * best_score / match_score; + } while (best_cost > b_last); + } else { + common::logger->trace("Using fallback chaining with bandwidth {}", config.chaining_bandwidth); + // otherwise, use this algorithm + for (AnchorIt i = begin + 1; i != end; ++i) { + AnchorIt j = i - std::min(static_cast(i - begin), config.chaining_bandwidth); + anchor_connector(*i, std::numeric_limits::max(), j, i, + chain_scores + (j - begin), + make_anchor_connector(i)); + } + } + }; + + auto cur_begin = begin; + auto jt = chain_scores.begin(); + for (auto it = begin + 1; it != end; ++it) { + if (it->get_orientation() != (it - 1)->get_orientation() + || it->get_label_class() != (it - 1)->get_label_class()) { + forward_pass(cur_begin, it, jt); + jt += it - cur_begin; + cur_begin = it; + } + } + + forward_pass(cur_begin, end, jt); + + // backtracking + std::vector> best_chains; + best_chains.reserve(chain_scores.size()); + for (size_t i = 0; i < chain_scores.size(); ++i) { + const auto &[score, last, dist, traversed, queried] = chain_scores[i]; + + if (score > 0) { + assert(last <= end); + assert(last >= begin); + best_chains.emplace_back( + -score, + (begin + i)->get_end_clipping(), + (begin + i)->get_orientation(), + i + ); + } + } + + std::sort(best_chains.begin(), best_chains.end()); + + sdsl::bit_vector used(chain_scores.size()); + + tsl::hopscotch_set label_classes; + AnchorChain last_chain; + + auto flush_chains = [&]() { + if (last_chain.size()) { + if (label_classes.size() == 1) { + extend_chain(last_chain, anchor_extender, callback, terminate); + } else { + AnnotationBuffer::Columns merged_columns; + for (auto label_class : label_classes) { + if (merged_columns.empty()) { + merged_columns = anno_buffer->get_cached_column_set(label_class); + } else { + const auto &next_columns = anno_buffer->get_cached_column_set(label_class); + AnnotationBuffer::Columns next_merged; + std::set_union(merged_columns.begin(), merged_columns.end(), + next_columns.begin(), next_columns.end(), + std::back_inserter(next_merged)); + std::swap(next_merged, merged_columns); + } + } + auto merged_label_class = anno_buffer->cache_column_set(std::move(merged_columns)); + std::vector relabeled_anchors; + relabeled_anchors.reserve(last_chain.size()); + for (const auto &[it, dist] : last_chain) { + relabeled_anchors.emplace_back(*it).set_label_class(merged_label_class); + } + for (size_t i = 0; i < last_chain.size(); ++i) { + static_assert(std::is_same_v::iterator>); + last_chain[i].first = relabeled_anchors.begin() + i; + } + extend_chain(last_chain, anchor_extender, callback, terminate); + } + } + + last_chain.clear(); + label_classes.clear(); + }; + + for (auto [nscore, end_clipping, orientation, i] : best_chains) { + if (terminate()) + return end; + + if (used[i]) + continue; + + AnchorChain chain; + std::vector scores; + AnchorIt last_anchor = begin + i; + auto [score, last, dist, traversed, queried] = chain_scores[i]; + chain.emplace_back(last_anchor, dist); + assert(score == -nscore); + scores.emplace_back(score); + while (last != end) { + last_anchor = last; + std::tie(score, last, dist, traversed, queried) = chain_scores[last - begin]; + assert(last_anchor != end); + chain.emplace_back(last_anchor, dist); + scores.emplace_back(score); + if (used[last_anchor - begin]) { + chain.clear(); + scores.clear(); + break; + } + } + + std::reverse(chain.begin(), chain.end()); + std::reverse(scores.begin(), scores.end()); + + if (chain.size() && extension_starter(chain, scores)) { + if (!anno_buffer) { + extend_chain(chain, anchor_extender, callback, terminate); + } else { + if (label_classes.empty() || !ChainEqual()(last_chain, chain)) { + flush_chains(); + last_chain = chain; + } + label_classes.emplace(chain[0].first->get_label_class()); + } + for (const auto &[a_ptr, dist] : chain) { + used[a_ptr - begin] = true; + } + } + } + + flush_chains(); + + return end; +} + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_cigar.cpp b/metagraph/src/graph/alignment_redone/aln_cigar.cpp new file mode 100644 index 0000000000..ab7aee4f10 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_cigar.cpp @@ -0,0 +1,196 @@ +#include "aln_cigar.hpp" + +#include "graph/representation/succinct/boss.hpp" + +namespace mtg::graph::align_redone { + + +Cigar::Cigar(std::string_view cigar_str) { + std::string op_count; + for (char c : cigar_str) { + switch (c) { + case '=': + cigar_.emplace_back(Cigar::MATCH, std::stol(op_count)); + op_count.clear(); + break; + case 'X': + cigar_.emplace_back(Cigar::MISMATCH, std::stol(op_count)); + op_count.clear(); + break; + case 'I': + cigar_.emplace_back(Cigar::INSERTION, std::stol(op_count)); + op_count.clear(); + break; + case 'D': + cigar_.emplace_back(Cigar::DELETION, std::stol(op_count)); + op_count.clear(); + break; + case 'S': + cigar_.emplace_back(Cigar::CLIPPED, std::stol(op_count)); + op_count.clear(); + break; + case 'G': + cigar_.emplace_back(Cigar::NODE_INSERTION, std::stol(op_count)); + op_count.clear(); + break; + default: + op_count += c; + } + } +} + +std::string Cigar::to_string() const { + std::string cigar_string; + + for (const auto &pair : cigar_) { + cigar_string += std::to_string(pair.second) + opt_to_char(pair.first); + } + + return cigar_string; +} + +void Cigar::append(Operator op, LengthType num) { + if (!num) + return; + + if (cigar_.empty() || cigar_.back().first != op) { + cigar_.emplace_back(op, num); + } else { + cigar_.back().second += num; + } +} + +void Cigar::append(Cigar&& other) { + if (other.empty()) + return; + + append(other.cigar_.front().first, other.cigar_.front().second); + cigar_.insert(cigar_.end(), std::next(other.cigar_.begin()), other.cigar_.end()); +} + +bool Cigar::is_valid(std::string_view reference, std::string_view query) const { + auto ref_it = reference.begin(); + auto alt_it = query.begin(); + + for (size_t i = 0; i < cigar_.size(); ++i) { + const auto &op = cigar_[i]; + if (!op.second) { + std::cerr << "Empty operation found in CIGAR" << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + switch (op.first) { + case CLIPPED: { + if (op.second > static_cast(query.end() - alt_it)) { + std::cerr << "Query too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + alt_it += op.second; + } break; + case MATCH: + case MISMATCH: { + if (op.second > static_cast(reference.end() - ref_it)) { + std::cerr << "Reference too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + if (op.second > static_cast(query.end() - alt_it)) { + std::cerr << "Query too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + if (std::equal(ref_it, ref_it + op.second, alt_it) + == (op.first != Cigar::MATCH)) { + std::cerr << "Mismatch despite MATCH in CIGAR" << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + ref_it += op.second; + alt_it += op.second; + } break; + case INSERTION: { + if (i && cigar_[i - 1].first == DELETION) { + std::cerr << "INSERTION after DELETION" << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + if (op.second > static_cast(query.end() - alt_it)) { + std::cerr << "Query too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + alt_it += op.second; + } break; + case DELETION: { + if (i && cigar_[i - 1].first == INSERTION + && (ref_it == reference.end() || *ref_it != boss::BOSS::kSentinel)) { + std::cerr << "DELETION after INSERTION" << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + if (op.second > static_cast(reference.end() - ref_it)) { + std::cerr << "Reference too short after " + << Cigar::opt_to_char(op.first) << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + ref_it += op.second; + } break; + case NODE_INSERTION: { + // do nothing + } break; + } + } + + if (ref_it != reference.end()) { + std::cerr << "Reference end not reached" << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + if (alt_it != query.end()) { + std::cerr << "Query end not reached" << std::endl + << to_string() << std::endl + << std::string(get_clipping(), ' ') << reference << std::endl + << query << std::endl; + return false; + } + + return true; +} + +} // namespace mtg::graph::align diff --git a/metagraph/src/graph/alignment_redone/aln_cigar.hpp b/metagraph/src/graph/alignment_redone/aln_cigar.hpp new file mode 100644 index 0000000000..19e29e8d25 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_cigar.hpp @@ -0,0 +1,109 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace mtg::graph::align_redone { + +class Cigar { + public: + enum Operator : int8_t { + CLIPPED = 0, + MISMATCH, + MATCH, + DELETION, + INSERTION, + NODE_INSERTION + }; + + typedef uint32_t LengthType; + typedef std::pair value_type; + + Cigar(Operator op = CLIPPED, LengthType num = 0) + : cigar_(num ? 1 : 0, std::make_pair(op, num)) { } + + // See section 1.4 in https://samtools.github.io/hts-specs/SAMv1.pdf for + // a specification of the CIGAR string format. + // e.g., 3=1X2I3D for 3 matches, 1 mismatch, 2 insertions, 3 deletions + // The symbol 'G' is introduced to indicate the insertion of a graph node. + Cigar(std::string_view cigar_str); + + size_t size() const { return cigar_.size(); } + bool empty() const { return cigar_.empty(); } + + std::string to_string() const; + + void append(Operator op, LengthType num = 1); + void append(Cigar&& other); + + LengthType trim_clipping() { + if (cigar_.size() && cigar_.front().first == CLIPPED) { + LengthType ret_val = cigar_.front().second; + cigar_.erase(cigar_.begin(), cigar_.begin() + 1); + return ret_val; + } else { + return 0; + } + } + + LengthType trim_end_clipping() { + if (cigar_.size() && cigar_.back().first == CLIPPED) { + LengthType ret_val = cigar_.back().second; + cigar_.pop_back(); + return ret_val; + } else { + return 0; + } + } + + LengthType get_clipping() const { + return cigar_.size() && cigar_.front().first == CLIPPED ? cigar_.front().second : 0; + } + + LengthType get_end_clipping() const { + return cigar_.size() && cigar_.back().first == CLIPPED ? cigar_.back().second : 0; + } + + std::vector& data() { return cigar_; } + const std::vector& data() const { return cigar_; } + + bool operator==(const Cigar &other) const { return cigar_ == other.cigar_; } + bool operator!=(const Cigar &other) const { return !(*this == other); } + + size_t get_num_matches() const { + return std::accumulate(cigar_.begin(), cigar_.end(), size_t{0}, + [&](size_t old, const value_type &op) { + return old + (op.first == MATCH) * op.second; + }); + } + + size_t get_num_query() const { + return std::accumulate(cigar_.begin(), cigar_.end(), size_t{0}, + [&](size_t old, const value_type &op) { + return old + (op.first == MATCH || op.first == MISMATCH || op.first == INSERTION || op.first == CLIPPED) * op.second; + }); + } + + // Return true if the cigar is valid + bool is_valid(std::string_view reference, std::string_view query) const; + + bool is_exact_match(size_t query_size) const { + return cigar_.size() == 1 && cigar_.front() == value_type{ MATCH, query_size }; + } + + static constexpr char opt_to_char(Cigar::Operator op) { return op_str_[op]; } + + private: + static constexpr char op_str_[] = "SX=DIG"; + std::vector cigar_; +}; + + + + +} // namespace mtg::graph::align diff --git a/metagraph/src/graph/alignment_redone/aln_match.cpp b/metagraph/src/graph/alignment_redone/aln_match.cpp new file mode 100644 index 0000000000..8547dc202f --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_match.cpp @@ -0,0 +1,124 @@ +#include "aln_match.hpp" + +namespace mtg::graph::align_redone { + +using mtg::common::logger; + +bool operator==(const Match &a, const Match &b) { + return a.get_orientation() == b.get_orientation() + && a.get_query().begin() == b.get_query().begin() + && a.get_query().end() == b.get_query().end() + && a.get_clipping() == b.get_clipping() + && a.get_end_clipping() == b.get_end_clipping() + && a.get_end_trim() == b.get_end_trim() + && a.get_path() == b.get_path(); +} + +bool operator==(const Anchor &a, const Anchor &b) { + return static_cast(a) == static_cast(b) + && a.get_label_class() == b.get_label_class(); +} + +bool operator==(const Alignment &a, const Alignment &b) { + return static_cast(a) == static_cast(b) + && a.get_cigar() == b.get_cigar(); +} + +std::ostream& operator<<(std::ostream &out, const Match &a) { + Cigar cigar = a.generate_cigar(); + out << fmt::format("{}\t{}\t*\t{}\t{}\t{}", + !a.get_orientation() ? "+" : "-", + a.get_spelling(), + cigar.get_num_matches(), + cigar.to_string(), + a.get_end_trim()); + return out; +} + +size_t Anchor::trim_end() { + assert(get_spelling().size() == get_path_spelling().size() - get_end_trim()); + if (suffix_.empty()) + return 0; + + size_t num_trimmed = 0; + +#ifndef NDEBUG + std::string old_spelling(get_spelling()); +#endif + + while (suffix_.size() && path_.size() > 1) { + suffix_.pop_back(); + path_.pop_back(); + ++num_trimmed; + assert(get_spelling().size() == get_path_spelling().size() - get_end_trim()); + assert(get_spelling() == old_spelling); + } + + assert(get_spelling() == old_spelling); + + return num_trimmed; +} + +size_t Alignment::trim_end() { + assert(get_spelling().size() == get_path_spelling().size() - get_end_trim()); + if (end_trim_ == 0) + return 0; + +#ifndef NDEBUG + std::string old_spelling(get_spelling()); +#endif + + size_t num_trimmed = 0; + while (end_trim_ && path_.size() > 1) { + assert(path_spelling_.size()); + path_spelling_.pop_back(); + path_.pop_back(); + --end_trim_; + ++num_trimmed; + assert(get_spelling().size() == get_path_spelling().size() - get_end_trim()); + assert(get_spelling() == old_spelling); + } + + assert(get_spelling() == old_spelling); + + return num_trimmed; +} + +void Anchor::append(const Anchor &other, const DeBruijnGraph *graph) { + assert(query_.data() == other.query_.data()); + assert(query_.data() + query_.size() == other.query_.data() + other.query_.size()); + assert(label_class_ == other.label_class_); + + if (empty()) { + *this = other; + return; + } + + assert(other.seed_.data() >= seed_.data()); + + if (seed_.data() + seed_.size() >= other.seed_.data() + other.seed_.size()) + return; + + if (seed_.data() == other.seed_.data()) { + assert(coord_ == other.coord_); + *this = other; + return; + } + + auto next = path_.begin() + (other.seed_.data() - seed_.data()); + assert(next <= path_.end()); + + seed_ = std::string_view( + seed_.data(), + other.seed_.data() + other.seed_.size() - seed_.data() + ); + + path_.erase(next, path_.end()); + path_.insert(path_.end(), other.path_.begin(), other.path_.end()); + suffix_ = other.suffix_; + + assert(!graph || is_spelling_valid(*graph)); +} + +} // namespace mtg::graph::align + diff --git a/metagraph/src/graph/alignment_redone/aln_match.hpp b/metagraph/src/graph/alignment_redone/aln_match.hpp new file mode 100644 index 0000000000..828b11dc58 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_match.hpp @@ -0,0 +1,216 @@ +#pragma once + +#include +#include + +#include "aln_cigar.hpp" +#include "aligner_config.hpp" +#include "annotation_buffer.hpp" + +#include "graph/alignment/alignment.hpp" +#include "graph/annotated_dbg.hpp" +#include "annotation/int_matrix/base/int_matrix.hpp" +#include "common/logger.hpp" + +namespace mtg::graph::align_redone { + +using align::spell_path; + +class Match { + public: + using node_index = SequenceGraph::node_index; + using score_t = DBGAlignerConfig::score_t; + + static constexpr score_t ninf = DBGAlignerConfig::ninf; + + Match() : orientation_(false) {} + + virtual ~Match() {} + + std::string_view get_query() const { return query_; } + std::string_view get_seed() const { return seed_; } + + size_t size() const { return path_.size(); } + + const std::vector& get_path() const { return path_; } + + virtual std::string_view get_spelling() const = 0; + virtual std::string get_path_spelling() const = 0; + virtual std::string_view get_trim_spelling() const = 0; + + virtual Cigar generate_cigar() const = 0; + + bool get_orientation() const { return orientation_; } + + virtual size_t get_end_trim() const = 0; + virtual size_t trim_end() = 0; + + bool empty() const { return path_.empty(); } + + size_t get_clipping() const { return seed_.data() - query_.data(); } + size_t get_end_clipping() const { + return query_.data() + query_.size() - (seed_.data() + seed_.size()); + } + + bool is_spelling_valid(const DeBruijnGraph &graph) const { + return get_path_spelling() == spell_path(graph, path_); + } + + protected: + Match(std::string_view query, + size_t begin, + size_t end, + bool orientation, + std::vector&& path) + : query_(query), orientation_(orientation), seed_(query_.substr(begin, end - begin)), + path_(std::move(path)) {} + + std::string_view query_; + bool orientation_; + std::string_view seed_; + std::vector path_; +}; + +inline Match::score_t score_match(const Match &match, const DBGAlignerConfig &config) { + return config.score_cigar(match.get_spelling(), match.get_query(), match.generate_cigar()); +} + +bool operator==(const Match &a, const Match &b); + +std::ostream& operator<<(std::ostream &out, const Match &a); + +class Anchor : public Match { + public: + using label_class_t = AnnotationBuffer::label_class_t; + using coord_t = AnnotationBuffer::coord_t; + static constexpr label_class_t nannot = AnnotationBuffer::nannot; + static constexpr coord_t ncoord = AnnotationBuffer::ncoord; + + Anchor() : Match(), label_class_(nannot), coord_(ncoord) {} + + Anchor(std::string_view query, + size_t begin, + size_t end, + bool orientation, + std::vector&& path, + std::string&& suffix = "") + : Match(query, begin, end, orientation, std::move(path)), + label_class_(nannot), coord_(ncoord), suffix_(std::move(suffix)) { + assert(get_spelling() == get_seed()); + } + + std::string_view get_spelling() const override final { return seed_; } + std::string get_path_spelling() const override final { return std::string(seed_) + suffix_; } + std::string_view get_trim_spelling() const override final { return suffix_; } + + Cigar generate_cigar() const override final { + Cigar cigar; + cigar.append(Cigar::CLIPPED, get_clipping()); + cigar.append(Cigar::MATCH, seed_.size()); + cigar.append(Cigar::CLIPPED, get_end_clipping()); + return cigar; + } + + size_t get_end_trim() const override final { return suffix_.size(); } + size_t trim_end() override final; + + void append(const Anchor &other, + const DeBruijnGraph *graph = nullptr); + + label_class_t get_label_class() const { return label_class_; } + int64_t get_coord() const { return coord_; } + + void set_label_class(label_class_t label_class) { + assert(label_class); + label_class_ = label_class; + } + + void set_coord(int64_t coord) { coord_ = coord; } + + private: + label_class_t label_class_; + int64_t coord_; + std::string suffix_; +}; + +bool operator==(const Anchor &a, const Anchor &b); + +class Alignment : public Match { + public: + using label_class_t = Anchor::label_class_t; + using coord_t = Anchor::coord_t; + using CoordinateSet = AnnotationBuffer::CoordinateSet; + + Alignment() : Alignment(Anchor()) {} + + Alignment(const DeBruijnGraph &graph, + std::string_view query, + bool orientation, + std::vector&& path, + Cigar&& cigar, + std::string&& path_spelling, + size_t end_trim = 0, + label_class_t label_class = Anchor::nannot, + CoordinateSet&& coords = {}) + : Match(query, cigar.get_clipping(), query.size() - cigar.get_end_clipping(), + orientation, std::move(path)), + end_trim_(end_trim), + path_spelling_(std::move(path_spelling)), + cigar_(std::move(cigar)), + label_classes_(path_.size(), label_class), + coords_(std::move(coords)) { + assert(path_spelling_ == spell_path(graph, path_)); + assert(get_spelling().size() + get_end_trim() == path_spelling_.size()); + assert(cigar_.get_clipping() == get_clipping()); + assert(cigar_.get_end_clipping() == get_end_clipping()); + assert(label_class); + assert(label_class != Anchor::nannot || coords_.empty()); + } + + Alignment(const Anchor &anchor) + : Match(anchor.get_query(), + anchor.get_clipping(), + anchor.get_clipping() + anchor.get_seed().size(), + anchor.get_orientation(), + std::vector(anchor.get_path())), + end_trim_(anchor.get_end_trim()), + path_spelling_(anchor.get_path_spelling()), + cigar_(anchor.generate_cigar()), + label_classes_(path_.size(), anchor.get_label_class()) { + assert(get_spelling().size() + get_end_trim() == path_spelling_.size()); + if (anchor.get_coord() != Anchor::ncoord) + coords_.emplace_back(AnnotationBuffer::Tuple(1, anchor.get_coord())); + } + + const Cigar& get_cigar() const { return cigar_; } + Cigar generate_cigar() const override final { return cigar_; } + + size_t get_end_trim() const override final { return end_trim_; } + size_t trim_end() override final; + + std::string_view get_spelling() const override final { + return std::string_view(path_spelling_.c_str(), path_spelling_.size() - end_trim_); + } + + std::string get_path_spelling() const override final { return path_spelling_; } + std::string_view get_trim_spelling() const override final { + return std::string_view(path_spelling_.c_str() + path_spelling_.size() - end_trim_, + end_trim_); + } + + const std::vector& get_label_classes() const { return label_classes_; } + + void set_coords(CoordinateSet&& coords) { coords_ = std::move(coords); } + + private: + size_t end_trim_; + std::string path_spelling_; + Cigar cigar_; + std::vector label_classes_; + CoordinateSet coords_; +}; + +bool operator==(const Alignment &a, const Alignment &b); + +} // namespace mtg::graph::align + diff --git a/metagraph/src/graph/alignment_redone/aln_query.hpp b/metagraph/src/graph/alignment_redone/aln_query.hpp new file mode 100644 index 0000000000..14595d82db --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_query.hpp @@ -0,0 +1,34 @@ +#pragma once + +#include + +#include "graph/representation/base/sequence_graph.hpp" +#include "common/seq_tools/reverse_complement.hpp" + +namespace mtg::graph::align_redone { + +class Query { + public: + template + Query(const DeBruijnGraph &graph, Args&&... args) + : graph_(graph), fwd_(std::forward(args)...), rc_(generate_rc()) {} + + const std::string& get_query(bool orientation = false) const { + return !orientation ? fwd_ : rc_; + } + + const DeBruijnGraph& get_graph() const { return graph_; } + + private: + const DeBruijnGraph &graph_; + std::string fwd_; + std::string rc_; + + std::string generate_rc() { + std::string rc = fwd_; + reverse_complement(rc.begin(), rc.end()); + return rc; + } +}; + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/aln_seeder.cpp b/metagraph/src/graph/alignment_redone/aln_seeder.cpp new file mode 100644 index 0000000000..9f272fa093 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_seeder.cpp @@ -0,0 +1,2689 @@ +#include "aln_seeder.hpp" + +#include + +#include + +#include "aln_chainer.hpp" +#include "common/logger.hpp" +#include "common/vector_map.hpp" +#include "common/algorithms.hpp" +#include "graph/representation/succinct/boss.hpp" +#include "graph/representation/succinct/dbg_succinct.hpp" +#include "graph/representation/hash/dbg_sshash.hpp" +#include "graph/representation/canonical_dbg.hpp" + +namespace mtg::graph::align_redone { + +void align_query(const Query &query, + const Seeder &seeder, + const Extender &extender, + const std::function &callback, + bool connect_anchors_in_chain) { + align_query(query, seeder, + [&](Alignment&& base_path) { extender.extend(base_path, callback); }, + connect_anchors_in_chain); +} + +void align_query(const Query &query, + const Seeder &seeder, + const std::function &callback, + bool connect_anchors_in_chain) { + for (Alignment &path : seeder.get_inexact_anchors(connect_anchors_in_chain)) { + callback(std::move(path)); + } +} + +using TraverseCallback = std::function; +using EdgeCallback = std::function; +using Terminator = std::function; + +class AlignmentGraph { + public: + using node_index = DeBruijnGraph::node_index; + using label_class_t = Anchor::label_class_t; + + AlignmentGraph(const DeBruijnGraph &graph, + AnnotationBuffer *anno_buffer = nullptr, + label_class_t target = Anchor::nannot) + : graph_(graph), anno_buffer_(anno_buffer), target_(target) { + assert((target == Anchor::nannot) == !anno_buffer_); + } + + AnnotationBuffer* get_anno_buffer() const { return anno_buffer_; } + + std::pair traverse(node_index node, char c) const { + assert(c != boss::BOSS::kSentinel); + node_index next = graph_.traverse(node, c); + if (target_ == Anchor::nannot) + return std::make_pair(next, target_); + + if (!next) + return std::make_pair(next, 0); + + label_class_t next_labels = get_intersect_labels(next, target_); + assert(next_labels != Anchor::nannot); + if (!next_labels) + next = DeBruijnGraph::npos; + + return std::make_pair(next, next_labels); + } + + std::pair traverse_back(node_index node, char c) const { + assert(c != boss::BOSS::kSentinel); + node_index prev = graph_.traverse_back(node, c); + if (target_ == Anchor::nannot) + return std::make_pair(prev, target_); + + if (!prev) + return std::make_pair(prev, 0); + + label_class_t prev_labels = get_intersect_labels(prev, target_); + assert(prev_labels != Anchor::nannot); + if (!prev_labels) + prev = DeBruijnGraph::npos; + + return std::make_pair(prev, prev_labels); + } + + void traverse(node_index node, std::string_view seq, + const TraverseCallback &callback, + const Terminator &terminate = [](){ return false; }, + size_t min_dist_before_check_anno = 0) const { + if (target_ == Anchor::nannot) { + graph_.traverse(node, seq, + [&](node_index next) { + assert(has_labels(next, target_)); + callback(next, target_); + }, terminate); + return; + } + + std::vector forward_nodes; + graph_.traverse(node, seq, + [&](node_index next) { forward_nodes.emplace_back(next); }, + terminate); + if (forward_nodes.size() < min_dist_before_check_anno) + return; + + anno_buffer_->queue_path(forward_nodes); + anno_buffer_->fetch_queued_annotations(); + + label_class_t cur_target = target_; + for (node_index next : forward_nodes) { + cur_target = get_intersect_labels(next, cur_target); + assert(cur_target != Anchor::nannot); + if (cur_target) { + assert(has_labels(next, cur_target)); + callback(next, cur_target); + } else { + return; + } + } + } + + void traverse_back(node_index node, std::string_view prefix_seq, + const TraverseCallback &callback, + const Terminator &terminate = [](){ return false; }, + size_t min_dist_before_check_anno = 0) const { + if (target_ == Anchor::nannot) { + graph_.traverse_back(node, prefix_seq, + [&](node_index prev) { + assert(has_labels(prev, target_)); + callback(prev, target_); + }, + terminate); + return; + } + + std::vector backward_nodes; + graph_.traverse_back(node, prefix_seq, + [&](node_index prev) { backward_nodes.emplace_back(prev); }, + terminate); + if (backward_nodes.size() < min_dist_before_check_anno) + return; + + anno_buffer_->queue_path(std::vector(backward_nodes.rbegin(), + backward_nodes.rend())); + anno_buffer_->fetch_queued_annotations(); + + label_class_t cur_target = target_; + for (node_index prev : backward_nodes) { + cur_target = get_intersect_labels(prev, cur_target); + assert(cur_target != Anchor::nannot); + if (cur_target) { + assert(has_labels(prev, cur_target)); + callback(prev, cur_target); + } else { + return; + } + } + } + + void call_outgoing_kmers(node_index node, const EdgeCallback &callback) const { + if (target_ == Anchor::nannot) { + graph_.call_outgoing_kmers(node, [&](node_index next, char c) { callback(next, c, target_); }); + } else { + std::vector forward_n; + std::vector forward_c; + graph_.call_outgoing_kmers(node, [&](node_index next, char c) { + if (c != boss::BOSS::kSentinel) { + forward_n.emplace_back(next); + forward_c.emplace_back(c); + anno_buffer_->queue_path(std::vector{ next }); + } + }); + anno_buffer_->fetch_queued_annotations(); + + for (size_t i = 0; i < forward_n.size(); ++i) { + node_index next = forward_n[i]; + char c = forward_c[i]; + label_class_t next_target = get_intersect_labels(next, target_); + assert(next_target != Anchor::nannot); + if (next_target) { + assert(has_labels(next, next_target)); + callback(next, c, next_target); + } + } + } + } + + void call_incoming_kmers(node_index node, const EdgeCallback &callback) const { + if (target_ == Anchor::nannot) { + graph_.call_incoming_kmers(node, [&](node_index prev, char c) { callback(prev, c, target_); }); + } else { + std::vector backward_n; + std::vector backward_c; + graph_.call_incoming_kmers(node, [&](node_index prev, char c) { + if (c != boss::BOSS::kSentinel) { + backward_n.emplace_back(prev); + backward_c.emplace_back(c); + anno_buffer_->queue_path(std::vector{ prev }); + } + }); + anno_buffer_->fetch_queued_annotations(); + + for (size_t i = 0; i < backward_n.size(); ++i) { + node_index prev = backward_n[i]; + char c = backward_c[i]; + label_class_t prev_target = get_intersect_labels(prev, target_); + assert(prev_target != Anchor::nannot); + if (prev_target) { + assert(has_labels(prev, prev_target)); + callback(prev, c, prev_target); + } + } + } + } + + bool has_single_outgoing(node_index node) const { + if (target_ == Anchor::nannot) + return graph_.has_single_outgoing(node); + + size_t outdegree = 0; + adjacent_outgoing_nodes(node, [&](node_index, label_class_t) { ++outdegree; }); + return outdegree == 1; + } + + bool has_single_incoming(node_index node) const { + if (target_ == Anchor::nannot) { + return graph_.has_single_incoming(node); + } else { + size_t indegree = 0; + adjacent_incoming_nodes(node, [&](node_index, label_class_t) { ++indegree; }); + return indegree == 1; + } + } + + void adjacent_outgoing_nodes(node_index node, const TraverseCallback &callback) const { + if (target_ == Anchor::nannot) { + graph_.adjacent_outgoing_nodes(node, [&](node_index next) { callback(next, target_); }); + } else { + // TODO: replace this once we can guarantee that we won't get + // dummy nodes + call_outgoing_kmers(node, [&](node_index next, char, label_class_t next_target) { + callback(next, next_target); + }); + } + } + + void adjacent_incoming_nodes(node_index node, const TraverseCallback &callback) const { + if (target_ == Anchor::nannot) { + graph_.adjacent_incoming_nodes(node, [&](node_index prev) { callback(prev, target_); }); + } else { + // TODO: replace this once we can guarantee that we won't get + // dummy nodes + call_incoming_kmers(node, [&](node_index prev, char, label_class_t prev_target) { + callback(prev, prev_target); + }); + } + } + + bool has_labels(node_index node, label_class_t target_labels_id) const { + if (!anno_buffer_) { + assert(target_labels_id == Anchor::nannot); + return true; + } + + assert(target_labels_id != Anchor::nannot); + label_class_t node_labels_id = get_label_class(node); + if (node_labels_id == Anchor::nannot) + return false; + + if (node_labels_id == target_labels_id) + return true; + + const auto &node_labels = anno_buffer_->get_cached_column_set(node_labels_id); + const auto &target_labels = anno_buffer_->get_cached_column_set(target_labels_id); + + return utils::count_intersection(node_labels.begin(), node_labels.end(), + target_labels.begin(), target_labels.end()) == target_labels.size(); + } + + private: + const DeBruijnGraph &graph_; + AnnotationBuffer* const anno_buffer_; + label_class_t target_; + + label_class_t get_label_class(node_index node) const { + assert(anno_buffer_); + + label_class_t node_labels_id = anno_buffer_->get_labels_id(node); + if (node_labels_id == Anchor::nannot) { + anno_buffer_->queue_path(std::vector{ node }); + anno_buffer_->fetch_queued_annotations(); + node_labels_id = anno_buffer_->get_labels_id(node); + } + assert(node_labels_id != Anchor::nannot); + return node_labels_id; + } + + label_class_t get_intersect_labels(node_index node, label_class_t target_labels_id) const { + assert(anno_buffer_); + assert(target_labels_id != Anchor::nannot); + + label_class_t node_labels_id = get_label_class(node); + assert(node_labels_id != Anchor::nannot); + + if (node_labels_id == target_labels_id) + return node_labels_id; + + const auto &target_labels = anno_buffer_->get_cached_column_set(target_labels_id); + assert(std::find(target_labels.begin(), target_labels.end(), AnnotationBuffer::ncolumn) == target_labels.end()); + assert(std::all_of(target_labels.begin(), target_labels.end(), + [&](auto c) { return c < anno_buffer_->get_annotator().get_label_encoder().size(); })); + + const auto &node_labels = anno_buffer_->get_cached_column_set(node_labels_id); + assert(std::find(node_labels.begin(), node_labels.end(), AnnotationBuffer::ncolumn) == node_labels.end()); + assert(std::all_of(node_labels.begin(), node_labels.end(), + [&](auto c) { return c < anno_buffer_->get_annotator().get_label_encoder().size(); })); + + AnnotationBuffer::Columns columns; + std::set_intersection(node_labels.begin(), node_labels.end(), + target_labels.begin(), target_labels.end(), + std::back_inserter(columns)); + + return columns.size() ? anno_buffer_->cache_column_set(std::move(columns)) : 0; + } +}; + +std::vector ExactSeeder::get_anchors() const { + const DeBruijnGraph &graph = query_.get_graph(); + std::vector anchors; + + if (query_.get_query().size() < config_.min_seed_length) + return anchors; + + std::vector> mapped_nodes(2); + std::vector>> max_seeds(2); + + for (bool orientation : { false, true }) { + std::string_view this_query = query_.get_query(orientation); + auto &nodes = mapped_nodes[orientation]; + nodes = map_to_nodes_sequentially(graph, this_query); + + auto &cur_max_seeds = max_seeds[orientation]; + cur_max_seeds.resize(this_query.size() - config_.min_seed_length + 1); + + for (size_t begin = 0; begin < nodes.size(); ++begin) { + Match::node_index node = nodes[begin]; + if (node != DeBruijnGraph::npos) { + cur_max_seeds[begin][node] = Anchor( + this_query, + begin, begin + graph.get_k(), + orientation, + std::vector{ node } + ); + } + } + } + + for (bool orientation : { false, true }) { + std::string_view this_query = query_.get_query(orientation); + std::string_view rc_query = query_.get_query(!orientation); + auto &nodes = mapped_nodes[orientation]; + auto &cur_max_seeds = max_seeds[orientation]; + + if (config_.min_seed_length < graph.get_k()) { + const DBGSuccinct *dbg_succ = nullptr; + const auto *canonical = dynamic_cast(&graph); + if (canonical) { + dbg_succ = dynamic_cast(&canonical->get_graph()); + assert(!dynamic_cast(&canonical->get_graph())); + } else { + dbg_succ = dynamic_cast(&graph); + } + + if (dbg_succ && (!orientation || canonical)) { + const auto &boss = dbg_succ->get_boss(); + auto encoded = boss.encode(this_query); + + size_t k = graph.get_k(); + for (size_t i = 0; i + config_.min_seed_length <= this_query.size(); ++i) { + if (std::find(encoded.begin() + i, encoded.begin() + i + config_.min_seed_length, + boss.alph_size) != encoded.begin() + i + config_.min_seed_length) { + continue; + } + + using edge_index = boss::BOSS::edge_index; + edge_index first; + edge_index last; + auto it = encoded.begin(); + std::tie(first, last, it) = boss.index_range( + encoded.begin() + i, + encoded.begin() + i + config_.min_seed_length + ); + + size_t match_size = it - (encoded.begin() + i); + if (match_size < config_.min_seed_length) + continue; + + first = first == boss.select_last(1) ? 1 : boss.pred_last(first - 1) + 1; + + if (canonical && orientation) { + if (dbg_succ->in_graph(last)) { + size_t i_rc = this_query.size() - (i + match_size); + dbg_succ->adjacent_incoming_nodes(last, [&](DeBruijnGraph::node_index node) { + node = canonical->reverse_complement(node); + bool inserted = true; + auto it = max_seeds[!orientation][i_rc].find(node); + if (it != max_seeds[!orientation][i_rc].end() && match_size <= it->second.get_seed().size()) + return; + + std::string spelling = graph.get_node_sequence(node); + if (spelling.find(boss::BOSS::kSentinel) != std::string::npos) + return; + + assert(std::equal(rc_query.begin() + i_rc, + rc_query.begin() + i_rc + match_size, + spelling.begin(), spelling.begin() + match_size)); + + std::tie(it, inserted) = max_seeds[!orientation][i_rc].try_emplace(node, Anchor()); + it.value() = Anchor(rc_query, + i_rc, i_rc + match_size, + !orientation, + std::vector{ node }, + spelling.substr(match_size)); + assert(it->second.get_path_spelling().find(boss::BOSS::kSentinel) == std::string::npos); + assert(it->second.get_path_spelling().size() == graph.get_k()); + assert(it->second.get_spelling().size() + it->second.get_end_trim() == graph.get_k()); + assert(it->second.is_spelling_valid(graph)); + }); + } + } + + std::vector> traverse; + std::string suffix; + bool is_exact_match = true; + traverse.emplace_back(first, last, suffix, is_exact_match); + while (traverse.size()) { + std::tie(first, last, suffix, is_exact_match) = traverse.back(); + assert(boss.get_last(last)); + traverse.pop_back(); + + size_t cur_match_size = match_size + suffix.size(); + + assert(cur_match_size < k); + if (!orientation && cur_match_size == k - 1) { + assert(first == last || !boss.get_last(first)); + assert(boss.succ_last(first) == last); + for (edge_index node = first; node <= last; ++node) { + if (i < nodes.size() && node == nodes[i]) + continue; + + if (dbg_succ->in_graph(node)) { + char c = boss.decode(boss.get_W(node) % boss.alph_size); + if (c != boss::BOSS::kSentinel) { + auto first_mm = std::mismatch(this_query.begin() + i + match_size, this_query.end(), + suffix.begin(), suffix.end()).second; + bool last_char_matches = i + k <= this_query.size() && first_mm == suffix.end() && c == this_query[i + k - 1]; + size_t full_match_size = match_size + (first_mm - suffix.begin()) + last_char_matches; + + is_exact_match &= last_char_matches; + + if (!is_exact_match) { + auto [it, inserted] = cur_max_seeds[i].try_emplace(node, Anchor()); + if (full_match_size > it->second.get_seed().size()) { + it.value() = Anchor(this_query, + i, i + full_match_size, + orientation, + std::vector{ node }, + (suffix + c).substr(full_match_size - match_size)); + assert(it->second.get_path_spelling().find(boss::BOSS::kSentinel) == std::string::npos); + assert(it->second.get_path_spelling().size() == graph.get_k()); + assert(it->second.get_spelling().size() + it->second.get_end_trim() == graph.get_k()); + assert(it->second.is_spelling_valid(graph)); + } + } + } + } + } + + continue; + } + + boss.call_tighter_ranges(first, last, + [&](edge_index next_first, edge_index next_last, boss::BOSS::TAlphabet s) { + char c = boss.decode(s); + bool next_is_exact_match = is_exact_match && i + cur_match_size < this_query.size() && this_query[i + cur_match_size] == c; + // The reverse-complement pass only creates anchors from + // exact matches, which are complete at length k + bool keep = orientation + ? next_is_exact_match && cur_match_size + 1 < k + : !canonical || next_is_exact_match + || i >= nodes.size() || nodes[i] == DeBruijnGraph::npos; + if (keep) { + traverse.emplace_back(next_first, next_last, suffix + c, + next_is_exact_match); + } + if (canonical && orientation && next_is_exact_match) { + size_t next_match_size = cur_match_size + 1; + size_t i_rc = this_query.size() - (i + next_match_size); + + std::vector succs; + for (edge_index e = boss.succ_W(first, s); e <= last; e = e + 1 < boss.num_edges() ? boss.succ_W(e + 1, s) : boss.num_edges() + 1) { + succs.emplace_back(e); + } + boss::BOSS::TAlphabet splus = s + boss.alph_size; + for (edge_index e = boss.succ_W(first, splus); e <= last; e = e + 1 < boss.num_edges() ? boss.succ_W(e + 1, splus) : boss.num_edges() + 1) { + succs.emplace_back(e); + } + + for (edge_index e : succs) { + assert(e); + auto node = canonical->reverse_complement(e); + if (!graph.in_graph(node)) + continue; + + bool inserted = true; + auto it = max_seeds[!orientation][i_rc].find(node); + if (it != max_seeds[!orientation][i_rc].end() && next_match_size <= it->second.get_seed().size()) + continue; + + std::string spelling = graph.get_node_sequence(node); + if (spelling.find(boss::BOSS::kSentinel) != std::string::npos) + continue; + + assert(std::equal(rc_query.begin() + i_rc, + rc_query.begin() + i_rc + next_match_size, + spelling.begin(), spelling.begin() + next_match_size)); + + std::tie(it, inserted) = max_seeds[!orientation][i_rc].try_emplace(node, Anchor()); + it.value() = Anchor(rc_query, + i_rc, i_rc + next_match_size, + !orientation, + std::vector{ node }, + spelling.substr(next_match_size)); + assert(it->second.get_path_spelling().find(boss::BOSS::kSentinel) == std::string::npos); + assert(it->second.get_path_spelling().size() == graph.get_k()); + assert(it->second.get_spelling().size() + it->second.get_end_trim() == graph.get_k()); + assert(it->second.is_spelling_valid(graph)); + } + } + } + ); + } + } + } + } + } + + std::vector coverage(2, sdsl::bit_vector(query_.get_query(false).size())); + for (auto &cur_max_seeds : max_seeds) { + for (auto &bucket : cur_max_seeds) { + // As in the existing aligner, a query position with more than + // max_num_seeds_per_locus seeds shorter than k keeps none of them + auto is_short = [&](const Anchor &a) { return a.get_seed().size() < graph.get_k(); }; + bool drop_short = static_cast(std::count_if(bucket.begin(), bucket.end(), + [&](const auto &entry) { return is_short(entry.second); })) + > config_.max_num_seeds_per_locus; + + for (auto it = bucket.begin(); it != bucket.end(); ++it) { + if (drop_short && is_short(it->second)) + continue; + + auto &cov = coverage[it->second.get_orientation()]; + std::fill(cov.begin() + it->second.get_clipping(), + cov.begin() + it->second.get_clipping() + it->second.get_seed().size(), + true); + anchors.emplace_back(std::move(it.value())); + assert(anchors.back().get_path_spelling().find(boss::BOSS::kSentinel) == std::string::npos); + assert(anchors.back().get_path_spelling().size() == graph.get_k()); + assert(anchors.back().get_spelling().size() + anchors.back().get_end_trim() == graph.get_k()); + } + } + + if (graph.get_mode() != DeBruijnGraph::BASIC) + break; + } + + size_t num_covered = std::max(sdsl::util::cnt_one_bits(coverage[false]), + sdsl::util::cnt_one_bits(coverage[true])); + + if (static_cast(num_covered) / query_.get_query(false).size() < config_.min_exact_match) + anchors.clear(); + + return anchors; +} + +std::vector LabeledSeeder::get_anchors() const { + auto anchors = ExactSeeder::get_anchors(); + + common::logger->trace("Annotating anchors"); + + for (const auto &anchor : anchors) { + assert(anchor.get_path_spelling().find(boss::BOSS::kSentinel) == std::string::npos); + anno_buffer_.queue_path(anchor.get_path()); + } + + anno_buffer_.fetch_queued_annotations(); + + tsl::hopscotch_map coverages; + for (auto &anchor : anchors) { + assert(anchor.get_path().size() == 1); + if (const auto *labels_it = anno_buffer_.get_labels(anchor.get_path()[0])) { + for (AnnotationBuffer::Column c : *labels_it) { + auto [it, inserted] = coverages.try_emplace(c, sdsl::bit_vector()); + if (inserted) + it.value() = sdsl::bit_vector(query_.get_query(false).size()); + + std::fill(it.value().begin() + anchor.get_clipping(), + it.value().begin() + anchor.get_clipping() + anchor.get_seed().size(), + true); + } + } + } + + for (auto it = coverages.begin(); it != coverages.end(); ++it) { + size_t num_covered = sdsl::util::cnt_one_bits(it->second); + + if (static_cast(num_covered) / query_.get_query(false).size() < config_.min_exact_match) { + // if the coverage is too low, clear out the vector + it.value() = sdsl::bit_vector(); + } else { + // otherwise, we don't need the coverage anymore, so just make it tiny + it.value() = sdsl::bit_vector(1); + } + } + + std::vector labeled_anchors; + labeled_anchors.reserve(anchors.size()); + for (auto &anchor : anchors) { + auto [labels_it, coord_it] = anno_buffer_.get_labels_and_coords(anchor.get_path()[0]); + if (labels_it) { + // Copy out of the buffer before calling cache_column_set()/cache_column() + // below: those mutate anno_buffer_'s underlying column_sets_ storage and + // can reallocate it, which would invalidate labels_it mid-loop. + auto labels = *labels_it; + for (size_t i = 0; i < labels.size(); ++i) { + if (coverages[labels[i]].size()) { + // only pick labels that have sufficient coverage + anchor.set_label_class(anno_buffer_.cache_column(labels[i])); + if (coord_it) { + const auto &coords = (*coord_it)[i]; + for (auto coord : coords) { + labeled_anchors.emplace_back(anchor).set_coord(coord); + assert(labeled_anchors.back().get_path_spelling().find(boss::BOSS::kSentinel) == std::string::npos); + } + } else { + labeled_anchors.emplace_back(anchor); + assert(labeled_anchors.back().get_path_spelling().find(boss::BOSS::kSentinel) == std::string::npos); + } + } + } + } + anchor = Anchor(); + } + + return labeled_anchors; +} + +template +class OffsetVector { + public: + OffsetVector() : offset_(0) {} + + template + OffsetVector(size_t offset, Args&&... args) : offset_(offset), null_(T()), data_(std::forward(args)...) {} + + size_t size() const { return offset_ + data_.size(); } + bool empty() const { return offset_ == 0 && data_.empty(); } + + void resize(size_t new_size, const T &val = T()) { + if (new_size >= offset_) { + data_.resize(new_size - offset_, val); + return; + } + + new_size -= data_.size(); + data_.clear(); + assert(new_size <= offset_); + offset_ -= new_size; + } + + size_t offset() const { return offset_; } + + template + T& set(size_t i, Args&&... args) { + get(i) = T(std::forward(args)...); + return data_[i - offset_]; + } + + T& get(size_t i) { + if (i < offset_) { + data_.insert(data_.begin(), offset_ - i, null_); + offset_ = i; + } + + if (i - offset_ >= data_.size()) + data_.resize(i - offset_ + 1); + + return data_[i - offset_]; + } + + T& operator[](size_t i) { + assert(i >= offset_); + return data_[i - offset_]; + } + + const T& operator[](size_t i) const { + return i >= offset_ ? data_[i - offset_] : null_; + } + + using const_iterator = typename std::vector::const_iterator; + const_iterator begin() const { return data_.begin(); } + const_iterator end() const { return data_.end(); } + const_iterator cbegin() const { return data_.cbegin(); } + const_iterator cend() const { return data_.cend(); } + + private: + size_t offset_; + const T null_; + std::vector data_; +}; + +template +using WaveFront = OffsetVector>; // S(cost) + +template +using ScoreTable = std::vector>; + +DBGAlignerConfig::score_t cost_to_score(size_t cost, + size_t query_size, + size_t match_size, + DBGAlignerConfig::score_t match_score) { + DBGAlignerConfig::score_t double_score = match_score * (query_size + match_size) - cost; + assert(double_score % 2 == 0); + return double_score / 2; +} + +template +const Tuple* get_bucket(const ScoreTable &table, + size_t cost, + size_t query_dist, + DeBruijnGraph::node_index node) { + if (cost >= table.size()) + return nullptr; + + if (query_dist >= table[cost].size() || query_dist < table[cost].offset()) + return nullptr; + + const auto &table_slot = table[cost][query_dist]; + auto it = table_slot.find(node); + if (it == table_slot.end()) + return nullptr; + + return &it->second; +}; + +// dist, num_ops, last_node, last_op, last_char_of_this_node, num_matches +using SMap = std::tuple; + +// dist, num_ops +using EMap = std::tuple; + +// dist, num_ops, last_node +using FMap = std::tuple; + +template +void align_impl(const std::function &has_single_incoming, + const std::function &traverse_back, + const std::function &call_incoming_kmers, + const DBGAlignerConfig &config, + const DeBruijnGraph::node_index start_node, + size_t start_num_matches, + const StrItr query_window_begin, + const StrItr query_window_end, + size_t max_dist, + const std::function&&, Cigar&&, size_t, Anchor::label_class_t)> &callback, + const std::function &start_backtrack, + const std::function &terminate_branch, + const std::function &terminate, + Anchor::label_class_t start_target) { + using node_index = DeBruijnGraph::node_index; + + if (query_window_begin == query_window_end) + return; + + size_t query_size = query_window_end - query_window_begin; + + const auto query_window_rbegin = std::make_reverse_iterator(query_window_end); + const auto query_window_rend = std::make_reverse_iterator(query_window_begin); + + // derived from 2.4.1 in https://doi.org/10.1101/2022.01.12.476087 (v1) + DBGAlignerConfig::score_t match_score = config.match_score("A"); + DBGAlignerConfig::score_t mismatch_score = config.score_sequences("A", "T"); + + assert(config.gap_opening_penalty <= config.gap_extension_penalty); + assert(config.gap_extension_penalty < 0); + assert(mismatch_score < match_score); + + ssize_t mismatch_cost = (match_score - mismatch_score) * 2; + ssize_t gap_ext = match_score - 2 * config.gap_extension_penalty; + ssize_t gap_opn = match_score - 2 * config.gap_opening_penalty; + + assert(mismatch_cost > 0); + assert(gap_ext > 0); + assert(gap_opn >= gap_ext); + + assert(mismatch_cost % 2 == 0); + assert(gap_ext % 2 == 0); + assert(gap_opn % 2 == 0); + + common::logger->trace("x: {}\to: {}\te: {}", mismatch_cost, gap_opn, gap_ext); + + // S(cost, query_dist, node) = best_dist + ScoreTable S; // best cost + ScoreTable E; // best cost (last operation is insertion) + ScoreTable F; // best cost (last operation is deletion) + + auto set_value = [&]( + auto &table, + size_t cost, + size_t query_dist, + node_index node, + size_t dist, + node_index last_node, + size_t num_ops, + Cigar::Operator last_op, + char c, + size_t num_matches, + Anchor::label_class_t target) -> size_t { + assert(target); + assert(last_op == Cigar::MATCH + || last_op == Cigar::MISMATCH + || last_op == Cigar::INSERTION + || last_op == Cigar::DELETION); + bool inserted = false; + if (cost >= table.size()) { + inserted = true; + table.resize(cost + 1); + } + + if (query_dist >= table[cost].size()) { + inserted = true; + table[cost].resize(query_dist + 1); + } + + auto &bucket = table[cost].get(query_dist)[node]; + using table_t = std::decay_t; + static_assert(std::is_same_v + || std::is_same_v + || std::is_same_v); + + inserted |= dist > std::get<0>(bucket); + + if constexpr(std::is_same_v) { + inserted |= std::get<3>(bucket) == Cigar::CLIPPED || std::make_pair(dist, num_matches) > std::make_pair(std::get<0>(bucket), std::get<5>(bucket)); + } + + if constexpr(std::is_same_v) + inserted |= std::get<2>(bucket) == DeBruijnGraph::npos; + + if constexpr(std::is_same_v) + inserted |= std::get<1>(bucket) == 0; + + if (inserted) { + std::get<0>(bucket) = dist; + std::get<1>(bucket) = num_ops; + + if constexpr(std::is_same_v + || std::is_same_v) { + assert(target); + assert(dist >= num_ops); + std::get<2>(bucket) = last_node; + } + + if constexpr(std::is_same_v) { + assert(c != '\0'); + std::get<3>(bucket) = c; + std::get<4>(bucket) = target; + } + + if constexpr(std::is_same_v) { + std::get<3>(bucket) = last_op; + std::get<4>(bucket) = c; + std::get<5>(bucket) = num_matches; + std::get<6>(bucket) = target; + } + + inserted = true; + } + + if constexpr(std::is_same_v) { + assert(std::get<3>(bucket) == Cigar::MATCH + || std::get<3>(bucket) == Cigar::MISMATCH + || std::get<3>(bucket) == Cigar::INSERTION + || std::get<3>(bucket) == Cigar::DELETION); + } + + return inserted; + }; + + size_t num_explored_nodes = 0; + size_t max_explored_dist = 0; + { + size_t query_dist = 0; + DeBruijnGraph::node_index node = start_node; + SMap data(0, 0, start_node, Cigar::MATCH, *(query_window_rbegin - 1), start_num_matches, start_target); + auto &[best_dist, last_ext, last_node, last_op, c, num_matches, cur_target] = data; + + + if (!has_single_incoming(node, best_dist, start_target)) { + set_value( + S, + 0, + query_dist, + node, + best_dist, + start_node, + best_dist, + Cigar::MATCH, + best_dist > 0 ? *query_window_rbegin : '\0', + num_matches, + cur_target + ); + } + + traverse_back(node, query_window_begin, query_window_end, best_dist, cur_target, + [&](DeBruijnGraph::node_index prev, AlignmentGraph::label_class_t prev_target) { + auto &[best_dist, last_ext, last_node, last_op, c, num_matches, cur_target] = data; + ++num_explored_nodes; + ++max_explored_dist; + ++best_dist; + ++query_dist; + ++last_ext; + ++num_matches; + node = prev; + cur_target = prev_target; + + if (!has_single_incoming(node, best_dist, cur_target)) { + set_value( + S, + 0, + query_dist, + node, + best_dist, + start_node, + best_dist, + Cigar::MATCH, + best_dist > 0 ? *query_window_rbegin : '\0', + num_matches, + cur_target + ); + } + }, + [&]() { return terminate_branch(0, data, query_dist, node) || std::get<0>(data) == max_dist; } + ); + + + set_value( + S, + 0, + query_dist, + node, + best_dist, + start_node, + best_dist, + Cigar::MATCH, + best_dist > 0 ? *query_window_rbegin : '\0', + num_matches, + cur_target + ); + } + + + // As in the existing aligner, once this many nodes were explored, stop at + // the first entry that is not on the best path: entries are visited by + // increasing cost, so that is one which neither reaches further into the + // query than any before it, nor ties the cost of the first that got there + double max_explored_nodes = config.max_nodes_per_seq_char * query_size; + size_t max_query_dist = 0; + size_t max_query_dist_cost = 0; + + for (size_t cost = 0; cost < S.size(); ++cost) { + bool done = false; + using HeapEntry = std::tuple; + std::priority_queue queue; + for (size_t query_dist = S[cost].offset(); query_dist < S[cost].size(); ++query_dist) { + for (const auto &[node, data] : S[cost][query_dist]) { + size_t best_dist = std::get<0>(data); + queue.emplace(best_dist, query_dist, node); + } + } + + size_t best_dist; + size_t query_dist; + DeBruijnGraph::node_index node; + while (queue.size()) { + std::tie(best_dist, query_dist, node) = queue.top(); + queue.pop(); + assert(query_dist <= query_size); + + auto it = S[cost][query_dist].find(node); + if (it != S[cost][query_dist].end()) { + // outdated entry + if (std::get<0>(it->second) != best_dist) + continue; + + bool on_best_path = query_dist > max_query_dist + || (query_dist == max_query_dist && cost == max_query_dist_cost); + if (num_explored_nodes >= max_explored_nodes && !on_best_path) { + common::logger->trace("Node limit reached after {} nodes", num_explored_nodes); + done = true; + break; + } + if (query_dist > max_query_dist) { + max_query_dist = query_dist; + max_query_dist_cost = cost; + } + + Cigar::Operator last_op = std::get<3>(it->second); + size_t num_matches = std::get<5>(it->second); + Anchor::label_class_t cur_target = std::get<6>(it->second); + assert(last_op != Cigar::CLIPPED); + + + if (terminate(cost, it->second, query_dist, node)) { + done = true; + continue; + } + + if (terminate_branch(cost, it->second, query_dist, node)) + continue; + + // forward creation of insertions + if (!done && query_dist < query_size) { + // extend a previous insertion + if (const auto *ins_ext_bucket = get_bucket(E, cost, query_dist, node)) { + auto [last_dist, last_num_ops] = *ins_ext_bucket; + assert(last_num_ops > 0); + assert(query_dist >= last_num_ops); + ssize_t next_ext_cost = cost + gap_ext; + if (!terminate_branch(next_ext_cost, SMap(last_dist, 0, node, Cigar::INSERTION, '\0', 0, Anchor::nannot), query_dist + 1, node) + && set_value(E, next_ext_cost, query_dist + 1, node, last_dist, node, last_num_ops + 1, Cigar::INSERTION, '\0', 0, Anchor::nannot) + && set_value(S, next_ext_cost, query_dist + 1, node, last_dist, node, 0, Cigar::INSERTION, '\0', 0, cur_target)) { + } + } + + // open an insertion + if (last_op != Cigar::DELETION) { + ssize_t next_opn_cost = cost + gap_opn; + if (!terminate_branch(next_opn_cost, SMap(best_dist, 0, node, Cigar::INSERTION, '\0', 0, Anchor::nannot), query_dist + 1, node) + && set_value(E, next_opn_cost, query_dist + 1, node, best_dist, node, 1, Cigar::INSERTION, '\0', 0, Anchor::nannot) + && set_value(S, next_opn_cost, query_dist + 1, node, best_dist, node, 0, Cigar::INSERTION, '\0', 0, cur_target)) { + } + } + } + + if (best_dist < max_dist) { + std::vector> prevs; + call_incoming_kmers(node, [&](auto prev, char c, AlignmentGraph::label_class_t prev_target) { + ++num_explored_nodes; + if (c != boss::BOSS::kSentinel) + prevs.emplace_back(prev, c, prev_target); + }, best_dist, cur_target); + + if (prevs.empty()) + continue; + + max_explored_dist = std::max(max_explored_dist, best_dist + 1); + + // deletion extension + if (!done) { + if (const auto *del_ext_bucket = get_bucket(F, cost, query_dist, node)) { + auto [last_dist, last_num_ops, last_node, del_char, last_target] = *del_ext_bucket; + assert(last_node != DeBruijnGraph::npos); + assert(last_num_ops > 0); + assert(last_dist >= last_num_ops); + if (last_dist < max_dist) { + ssize_t next_ext_cost = cost + gap_ext; + std::vector> last_prevs; + if (best_dist == last_dist) { + last_prevs = prevs; + } else { + call_incoming_kmers(node, [&](auto prev, char c, auto prev_target) { + ++num_explored_nodes; + if (c != boss::BOSS::kSentinel) + last_prevs.emplace_back(prev, c, prev_target); + }, last_dist, last_target); + } + + for (const auto &[prev, c, prev_target] : last_prevs) { + if (!terminate_branch(next_ext_cost, SMap(last_dist + 1, 0, node, Cigar::DELETION, c, 0, prev_target), query_dist, prev) + && set_value(F, next_ext_cost, query_dist, prev, last_dist + 1, node, last_num_ops + 1, Cigar::DELETION, c, 0, prev_target) + && set_value(S, next_ext_cost, query_dist, prev, last_dist + 1, node, 0, Cigar::DELETION, c, 0, prev_target)) { + } + } + } + } + + // deletion open + if (last_op != Cigar::INSERTION) { + ssize_t next_opn_cost = cost + gap_opn; + for (const auto &[prev, c, prev_target] : prevs) { + if (!terminate_branch(next_opn_cost, SMap(best_dist + 1, 0, node, Cigar::DELETION, c, 0, prev_target), query_dist, prev) + && set_value(F, next_opn_cost, query_dist, prev, best_dist + 1, node, 1, Cigar::DELETION, c, 0, prev_target) + && set_value(S, next_opn_cost, query_dist, prev, best_dist + 1, node, 0, Cigar::DELETION, c, 0, prev_target)) { + } + } + } + } + + if (query_dist < query_size) { + // match + auto local_query_window_begin = query_window_begin; + auto local_query_window_end = query_window_end; + local_query_window_end -= query_dist; + + std::reverse_iterator local_query_window_rbegin = std::make_reverse_iterator(local_query_window_end); + std::reverse_iterator local_query_window_rend = std::make_reverse_iterator(local_query_window_begin); + + for (size_t i = 0; i < prevs.size(); ++i) { + DeBruijnGraph::node_index prev = std::get<0>(prevs[i]); + char c = std::get<1>(prevs[i]); + AlignmentGraph::label_class_t prev_target = std::get<2>(prevs[i]); + + SMap data(best_dist + 1, 1, node, Cigar::MATCH, c, num_matches + 1, prev_target); + auto &[cur_best, num_ops, last_node, op, cur_c, cur_num_matches, cur_target] = data; + size_t cur_query_dist = query_dist + 1; + size_t cur_cost = cost; + if (c != *local_query_window_rbegin) { + cur_cost += mismatch_cost; + op = Cigar::MISMATCH; + cur_num_matches = 0; + } else if (last_op == Cigar::MATCH) { + continue; + } + + if (!has_single_incoming(prev, cur_best, cur_target)) { + if (set_value( + S, + cur_cost, + cur_query_dist, + prev, + cur_best, + node, + cur_query_dist - query_dist, + op, + c, + cur_num_matches, + cur_target + )) { + if (cur_cost == cost) { + queue.emplace(cur_best, cur_query_dist, prev); + } + } + } + + traverse_back( + prev, + local_query_window_rend.base(), + (local_query_window_rbegin + 1).base(), + cur_best, cur_target, + [&](DeBruijnGraph::node_index pprev, AlignmentGraph::label_class_t pprev_target) { + auto &[cur_best, num_ops, last_node, op, cur_c, cur_num_matches, cur_target] = data; + prev = pprev; + ++num_explored_nodes; + ++cur_best; + max_explored_dist = std::max(max_explored_dist, cur_best); + ++cur_query_dist; + ++cur_num_matches; + cur_target = pprev_target; + + if (!has_single_incoming(prev, cur_best, cur_target)) { + if (set_value( + S, + cur_cost, + cur_query_dist, + prev, + cur_best, + node, + cur_query_dist - query_dist, + op, + c, + cur_num_matches, + cur_target + )) { + if (cur_cost == cost) { + queue.emplace(cur_best, cur_query_dist, prev); + } + } + } + }, + [&]() { return terminate_branch(cur_cost, data, cur_query_dist, prev) || std::get<0>(data) == max_dist; } + ); + + if (set_value( + S, + cur_cost, + cur_query_dist, + prev, + cur_best, + node, + cur_query_dist - query_dist, + op, + c, + cur_num_matches, + cur_target + )) { + if (cur_cost == cost) { + queue.emplace(cur_best, cur_query_dist, prev); + } + } + } + } + } + } + } + + if (done) + break; + } + + // backtrack + for (size_t start_cost = 0; start_cost < S.size(); ++start_cost) { + std::vector> starts; + for (size_t start_query_dist = S[start_cost].offset(); start_query_dist < S[start_cost].size(); ++start_query_dist) { + assert(start_query_dist <= query_size); + const auto &bucket = S[start_cost][start_query_dist]; + for (const auto &[node, data] : bucket) { + const auto &[dist, last_ext, last_node, last_op, mismatch_char, num_matches, cur_target] = data; + assert(last_op != Cigar::CLIPPED); + if (start_backtrack(start_cost, data, start_query_dist, node)) + starts.emplace_back(dist, static_cast(dist) - start_query_dist, start_query_dist, node); + } + } + + auto traverse_back_char = [&traverse_back](DeBruijnGraph::node_index node, char c, size_t dist, Anchor::label_class_t cur_target) { + std::string c_str(1, c); + StrItr begin; + StrItr end; + static_assert(std::is_same_v + || std::is_same_v>); + if constexpr(std::is_same_v) { + begin = c_str.c_str(); + end = begin + 1; + } else if constexpr(std::is_same_v>) { + begin = std::make_reverse_iterator(c_str.c_str() + 1); + end = std::make_reverse_iterator(c_str.c_str()); + } + + DeBruijnGraph::node_index prev = DeBruijnGraph::npos; + traverse_back(node, begin, end, dist, cur_target, + [&](DeBruijnGraph::node_index pprev, AlignmentGraph::label_class_t prev_target) { + prev = pprev; + cur_target = prev_target; + }, + []() { return false; }); + return std::make_pair(prev, cur_target); + }; + + for (auto [dist, diag, query_dist, node] : starts) { + size_t cost = start_cost; + + const auto *bt = get_bucket(S, cost, query_dist, node); + assert(bt); + + std::vector path; + Cigar cigar; + std::string spelling; + Anchor::label_class_t final_target = std::get<6>(*bt); + Anchor::label_class_t cur_target = final_target; + + do { + assert(std::get<0>(*bt) == dist); + + if (std::get<3>(*bt) == Cigar::INSERTION) { + assert(std::get<1>(*bt) == 0); + assert(cost > 0); + const auto *et = get_bucket(E, cost, query_dist, node); + assert(et); + + auto [ins_dist, num_ins] = *et; + assert(query_dist >= num_ins); + assert(num_ins > 0); + assert(ins_dist == dist); + + size_t prev_cost = cost - gap_opn - (num_ins - 1) * gap_ext; + assert(prev_cost < cost); + size_t prev_query_dist = query_dist - num_ins; + + const auto *check_bt = get_bucket(S, prev_cost, prev_query_dist, node); + assert(check_bt); + assert(std::get<0>(*check_bt) == dist); + assert(std::get<3>(*check_bt) != Cigar::CLIPPED); + assert(std::get<6>(*check_bt) == cur_target); + + query_dist = prev_query_dist; + cost = prev_cost; + bt = check_bt; + cigar.append(Cigar::INSERTION, num_ins); + } else if (std::get<3>(*bt) == Cigar::DELETION) { + assert(std::get<1>(*bt) == 0); + assert(cost > 0); + const auto *ft = get_bucket(F, cost, query_dist, node); + assert(ft); + + auto [del_dist, num_del, prev_node, del_char, del_target] = *ft; + assert(prev_node != DeBruijnGraph::npos); + assert(del_char != '\0'); + assert(del_dist == dist); + assert(dist >= num_del); + assert(del_target == cur_target); + + cigar.append(Cigar::DELETION, num_del); + + size_t prev_cost = cost; + size_t prev_dist = dist; + while (num_del > 1) { + path.emplace_back(node); + spelling += del_char; + + prev_cost -= gap_ext; + --prev_dist; + node = prev_node; + --num_del; + + const auto *prev_ft = get_bucket(F, prev_cost, query_dist, node); + assert(prev_ft); + assert(std::get<0>(*prev_ft) == prev_dist); + assert(std::get<1>(*prev_ft) == num_del); + ft = prev_ft; + prev_node = std::get<2>(*ft); + del_char = std::get<3>(*ft); + cur_target = std::get<4>(*ft); + } + + path.emplace_back(node); + spelling += del_char; + + prev_cost -= gap_opn; + --prev_dist; + + assert(prev_cost < cost); + + const auto *check_bt = get_bucket(S, prev_cost, query_dist, prev_node); + assert(check_bt); + assert(std::get<0>(*check_bt) == prev_dist); + assert(std::get<3>(*check_bt) != Cigar::CLIPPED); + assert(std::get<6>(*check_bt) == cur_target); + + dist = prev_dist; + cost = prev_cost; + node = prev_node; + bt = check_bt; + } else if (std::get<3>(*bt) == Cigar::MATCH || std::get<3>(*bt) == Cigar::MISMATCH) { + auto [cur_dist, num_match, last_node, last_op, mismatch_char, num_matches, mm_target] = *bt; + + if (!num_match) { + assert(last_op == Cigar::MATCH); + bt = nullptr; + break; + } + + assert(query_dist >= num_match); + assert(dist >= num_match); + assert(cur_dist == dist); + assert(cur_target == mm_target); + + node_index traverse_node = last_node; + + cigar.append(Cigar::MATCH, num_match - 1); + cigar.append(last_op); + size_t prev_dist = dist - num_match; + size_t prev_query = query_dist - num_match; + + // reconstruct the backwards traversal from last_node to node + auto it = query_window_rbegin + prev_query; + std::ignore = query_window_rend; + assert(it != query_window_rend); + size_t cur_path_size = path.size(); + size_t num_traversed = 0; + if (std::get<3>(*bt) == Cigar::MISMATCH) { + assert(mismatch_char != *it); + std::tie(traverse_node, cur_target) = traverse_back_char(traverse_node, mismatch_char, prev_dist, cur_target); + assert(traverse_node != DeBruijnGraph::npos); + path.emplace_back(traverse_node); + spelling += mismatch_char; + ++it; + --num_match; + ++num_traversed; + } else { + assert(mismatch_char == *it); + } + + std::reverse_iterator it_end = it + num_match; + assert(it_end <= query_window_rend); + size_t old_path_size = path.size(); + traverse_back(traverse_node, it_end.base(), it.base(), + prev_dist + num_traversed, + cur_target, + [&](DeBruijnGraph::node_index prev, AlignmentGraph::label_class_t prev_target) { + traverse_node = prev; + cur_target = prev_target; + path.emplace_back(traverse_node); + }, + []() { return false; }); + std::ignore = old_path_size; + assert(path.size() == old_path_size + num_match); + assert(traverse_node == node); + spelling += std::string(it, it_end); + num_traversed += num_match; + num_match = 0; + + std::reverse(path.begin() + cur_path_size, path.end()); + std::reverse(spelling.begin() + cur_path_size, spelling.end()); + + size_t prev_cost = cost - (last_op == Cigar::MATCH ? 0 : mismatch_cost); + assert(prev_cost <= cost); + + const auto *check_bt = bt; + check_bt = nullptr; + if (prev_cost || prev_query || prev_dist) { + check_bt = get_bucket(S, prev_cost, prev_query, last_node); + assert(check_bt); + assert(std::get<0>(*check_bt) == prev_dist); + assert(std::get<3>(*check_bt) != Cigar::CLIPPED); + assert(std::get<6>(*check_bt) == cur_target); + } else { + assert(last_node == start_node); + } + + dist = prev_dist; + cost = prev_cost; + query_dist = prev_query; + node = last_node; + + if (!check_bt) + break; + + bt = check_bt; + } else { + assert(false && "Different operation found"); + throw std::runtime_error("Inf. loop"); + } + + } while (bt); + + assert(dist == 0); + assert(cost == 0); + assert(query_dist == 0); + callback(std::move(spelling), std::move(path), std::move(cigar), start_cost, final_target); + } + } + common::logger->trace("Explored {} nodes. q: {}\td: {} / {}", + num_explored_nodes, query_size, max_explored_dist, max_dist); +} + +template +bool overlap_with_diff(const T1 &tuple1, const T2 &tuple2, ssize_t diff) { + auto a_begin = tuple1.begin(); + const auto a_end = tuple1.end(); + auto b_begin = tuple2.begin(); + const auto b_end = tuple2.end(); + + assert(std::is_sorted(a_begin, a_end)); + assert(std::is_sorted(b_begin, b_end)); + + while (a_begin != a_end && b_begin != b_end) { + if (*a_begin + diff == *b_begin) + return true; + + if (*a_begin + diff < *b_begin) { + ++a_begin; + } else { + ++b_begin; + } + } + + return false; +} + +Anchor::label_class_t subset_coords(AnnotationBuffer &anno_buffer, + Anchor::label_class_t source, + const AnnotationBuffer::CoordinateSet &source_coords, + Anchor::label_class_t target, + DeBruijnGraph::node_index target_node, + ssize_t coord_offset) { + assert(source); + assert(target); + assert(source != AnnotationBuffer::nannot); + assert(target != AnnotationBuffer::nannot); + + const auto &source_columns = anno_buffer.get_cached_column_set(source); + assert(source_columns.size() == source_coords.size()); + + const auto &target_columns = anno_buffer.get_cached_column_set(target); + const auto [target_columns_full, target_coords_full] = anno_buffer.get_labels_and_coords(target_node); + assert(target_columns_full); + assert(target_columns_full->size() >= target_columns.size()); + assert(target_coords_full); + assert(target_coords_full->size() == target_columns_full->size()); + +#ifndef NDEBUG + auto it = target_columns.begin(); +#endif + + AnnotationBuffer::Columns intersection; + utils::match_indexed_values( + source_columns.begin(), source_columns.end(), source_coords.begin(), + target_columns_full->begin(), target_columns_full->end(), target_coords_full->begin(), + [&](AnnotationBuffer::Column c, const AnnotationBuffer::Tuple &source_coords, const AnnotationBuffer::Tuple &target_coords) { + assert(source_coords.size()); + assert(target_coords.size()); +#ifndef NDEBUG + it = std::find(it, target_columns.end(), c); + assert(it != target_columns.end()); +#endif + if (overlap_with_diff(source_coords, target_coords, coord_offset)) + intersection.emplace_back(c); + } + ); + + assert(it + 1 == target_columns.end()); + + if (intersection.size() == target_columns.size()) + return target; + + return anno_buffer.cache_column_set(std::move(intersection)); +} + +void align_bwd(const DeBruijnGraph &base_graph, + const DBGAlignerConfig &config, + const DeBruijnGraph::node_index start_node, + size_t num_matches, + std::string_view query_window, + size_t max_dist, + const std::function&&, Cigar&&, size_t, Anchor::label_class_t)> &callback, + const std::function &start_backtrack, + const std::function &terminate_branch, + const std::function &terminate, + AnnotationBuffer *anno_buffer = nullptr, + Anchor::label_class_t target = Anchor::nannot, + const AnnotationBuffer::CoordinateSet &coords = {}) { + assert(!anno_buffer == (target == Anchor::nannot)); + assert(!anno_buffer || anno_buffer->node_has_labels(start_node, target)); + assert(!anno_buffer + || (!anno_buffer->has_coordinates() && coords.empty()) + || anno_buffer->get_cached_column_set(target).size() == coords.size()); + + using StrItr = std::string_view::iterator; + + align_impl( + [&](DeBruijnGraph::node_index node, size_t dist, Anchor::label_class_t cur_target) { + assert(cur_target); + AlignmentGraph aln_graph(base_graph, anno_buffer, cur_target); + if (coords.empty()) { + return aln_graph.has_single_incoming(node); + } else { + assert(anno_buffer); + size_t indegree = 0; + size_t next_dist = dist + 1; + aln_graph.call_incoming_kmers( + node, + [&](DeBruijnGraph::node_index prev, char, Anchor::label_class_t prev_target) { + if (subset_coords(*anno_buffer, target, coords, prev_target, prev, + -static_cast(next_dist))) + ++indegree; + } + ); + return indegree == 1; + } + }, + [&](DeBruijnGraph::node_index node, + StrItr begin, StrItr end, + size_t dist, Anchor::label_class_t cur_target, + const auto &callback, const auto &terminate) { + + assert(cur_target); + if (begin != end) { + AlignmentGraph aln_graph(base_graph, anno_buffer, cur_target); + if (coords.empty()) { + aln_graph.traverse_back( + node, + std::string_view(&*begin, end - begin), + callback, terminate + ); + } else { + assert(anno_buffer); + bool stop_early = false; + size_t next_dist = dist + 1; + aln_graph.traverse_back( + node, + std::string_view(&*begin, end - begin), + [&](DeBruijnGraph::node_index prev, Anchor::label_class_t prev_target) { + if (subset_coords(*anno_buffer, target, coords, prev_target, prev, + -static_cast(next_dist++))) { + callback(prev, prev_target); + } else { + stop_early = true; + } + }, + [&]() { return stop_early || terminate(); } + ); + } + } + }, + [&](DeBruijnGraph::node_index node, + const auto &callback, + size_t dist, Anchor::label_class_t cur_target) { + + assert(cur_target); + + AlignmentGraph aln_graph(base_graph, anno_buffer, cur_target); + if (coords.empty()) { + aln_graph.call_incoming_kmers(node, callback); + } else { + size_t next_dist = dist + 1; + aln_graph.call_incoming_kmers( + node, + [&](DeBruijnGraph::node_index prev, char c, Anchor::label_class_t prev_target) { + if (Anchor::label_class_t prev_target_s = subset_coords(*anno_buffer, target, coords, prev_target, prev, + -static_cast(next_dist))) + callback(prev, c, prev_target_s); + } + ); + } + }, + config, + start_node, + num_matches, + query_window.begin(), query_window.end(), + max_dist, + [&](auto&& spelling, auto&& path, auto&& cigar, size_t cost, Anchor::label_class_t final_target) { + assert(final_target); + assert(std::all_of(path.begin(), path.end(), + [&](auto node) { + return AlignmentGraph(base_graph, anno_buffer, target).has_labels(node, final_target); + })); + callback(std::move(spelling), std::move(path), std::move(cigar), cost, final_target); + }, + start_backtrack, + terminate_branch, + terminate, + target + ); +} + +void align_fwd(const DeBruijnGraph &base_graph, + const DBGAlignerConfig &config, + const DeBruijnGraph::node_index start_node, + size_t num_matches, + std::string_view query_window, + size_t max_dist, + const std::function&&, Cigar&&, size_t, Anchor::label_class_t)> &callback, + const std::function &start_backtrack, + const std::function &terminate_branch, + const std::function &terminate, + std::string_view suffix = "", + AnnotationBuffer *anno_buffer = nullptr, + Anchor::label_class_t target = Anchor::nannot, + const AnnotationBuffer::CoordinateSet &coords = {}) { + assert(!anno_buffer == (target == Anchor::nannot)); + assert(!anno_buffer || anno_buffer->node_has_labels(start_node, target)); + assert(!anno_buffer + || (!anno_buffer->has_coordinates() && coords.empty()) + || anno_buffer->get_cached_column_set(target).size() == coords.size()); + + using StrItr = std::reverse_iterator; + + align_impl( + [&](DeBruijnGraph::node_index node, size_t dist, Anchor::label_class_t cur_target) { + assert(cur_target); + assert(dist > suffix.size() || node == start_node); + if (dist < suffix.size()) + return true; + + AlignmentGraph aln_graph(base_graph, anno_buffer, cur_target); + if (coords.empty()) { + return aln_graph.has_single_outgoing(node); + } else { + assert(anno_buffer); + size_t outdegree = 0; + assert(dist >= suffix.size()); + size_t next_dist = dist - suffix.size() + 1; + aln_graph.call_outgoing_kmers( + node, + [&](DeBruijnGraph::node_index next, char, Anchor::label_class_t next_target) { + if (subset_coords(*anno_buffer, target, coords, next_target, next, next_dist)) + ++outdegree; + } + ); + return outdegree == 1; + } + }, + [&](DeBruijnGraph::node_index node, + StrItr rbegin, StrItr rend, + size_t dist, Anchor::label_class_t cur_target, + const auto &callback, const auto &terminate) { + + assert(cur_target); + assert(dist > suffix.size() || node == start_node); + + auto begin = rend.base(); + auto end = rbegin.base(); + assert(begin <= end); + + for (size_t i = dist; !terminate() && i < suffix.size() && begin != end; ++i, ++begin) { + if (*begin == suffix[i]) { + callback(node, cur_target); + } else { + return; + } + } + + if (begin != end) { + AlignmentGraph aln_graph(base_graph, anno_buffer, cur_target); + if (coords.empty()) { + aln_graph.traverse( + node, std::string_view(&*begin, end - begin), + callback, terminate + ); + } else { + assert(anno_buffer); + bool stop_early = false; + assert(dist >= suffix.size()); + size_t next_dist = dist - suffix.size() + 1; + aln_graph.traverse( + node, + std::string_view(&*begin, end - begin), + [&](DeBruijnGraph::node_index next, Anchor::label_class_t next_target) { + if (subset_coords(*anno_buffer, target, coords, next_target, next, next_dist++)) { + callback(next, next_target); + } else { + stop_early = true; + } + }, + [&]() { return stop_early || terminate(); } + ); + } + } + }, + [&](DeBruijnGraph::node_index node, const auto &callback, + size_t dist, Anchor::label_class_t cur_target) { + + assert(cur_target); + if (dist < suffix.size()) { + assert(node == start_node); + callback(node, suffix[dist], cur_target); + } else { + AlignmentGraph aln_graph(base_graph, anno_buffer, cur_target); + if (coords.empty()) { + aln_graph.call_outgoing_kmers(node, callback); + } else { + assert(dist >= suffix.size()); + size_t next_dist = dist - suffix.size() + 1; + aln_graph.call_outgoing_kmers( + node, + [&](DeBruijnGraph::node_index next, char c, Anchor::label_class_t next_target) { + if (Anchor::label_class_t next_target_s = subset_coords(*anno_buffer, target, coords, next_target, next, next_dist)) + callback(next, c, next_target_s); + } + ); + } + } + }, + config, + start_node, + num_matches, + query_window.rbegin(), query_window.rend(), + max_dist, + [&](auto&& spelling, auto&& path, auto&& cigar, size_t cost, Anchor::label_class_t final_target) { + assert(final_target); + path.resize(path.size() - std::min(path.size(), suffix.size())); + assert(std::all_of(path.begin(), path.end(), + [&](auto node) { + return AlignmentGraph(base_graph, anno_buffer, final_target).has_labels(node, target); + })); + std::reverse(spelling.begin(), spelling.end()); + std::reverse(path.begin(), path.end()); + std::reverse(cigar.data().begin(), cigar.data().end()); + callback(std::move(spelling), std::move(path), std::move(cigar), cost, final_target); + }, + start_backtrack, + terminate_branch, + terminate, + target + ); +} + +void Extender::extend(const Alignment &aln, + const std::function &callback) const { + // Chains whose anchors were not connected keep placeholder nodes between + // them, which extension cannot start from. Report them unextended. + const auto &aln_path = aln.get_path(); + if (std::find(aln_path.begin(), aln_path.end(), DeBruijnGraph::npos) != aln_path.end()) { + callback(Alignment(aln)); + return; + } + + DBGAlignerConfig::score_t match_score = config_.match_score("A"); + + auto set_coords = [&](Alignment &aln) { + assert(aln.get_label_classes().size()); + if (auto anno_buffer = make_aln_graph(aln.get_label_classes()[0]).get_anno_buffer()) { + if (anno_buffer->has_coordinates()) { + auto coords = anno_buffer->get_label_path_coords(aln.get_path(), aln.get_label_classes()); + assert(coords.size()); + aln.set_coords(std::move(coords)); + } + } + }; + + std::vector fwd_exts; + if (!aln.get_end_clipping()) { + fwd_exts.emplace_back(aln); + } else { + std::string_view query_window = aln.get_query(); + query_window.remove_prefix(aln.get_clipping() + aln.get_seed().size()); + + DBGAlignerConfig::score_t aln_score = score_match(aln, config_); + auto get_score = [&](size_t cost, size_t dist, size_t query_dist) { + return aln_score + cost_to_score(cost, query_dist, dist, match_score) + + (query_dist == query_window.size() ? config_.right_end_bonus : 0); + }; + + auto it = std::make_reverse_iterator(aln.get_cigar().data().end()); + auto end = std::make_reverse_iterator(aln.get_cigar().data().begin()); + std::ignore = end; + assert(it != end); + if (it->first == Cigar::CLIPPED) { + ++it; + assert(it != end); + } + size_t num_matches = it->first == Cigar::MATCH ? it->second : 0; + size_t max_dist = 0; + DBGAlignerConfig::score_t best_score = aln_score; + size_t max_query_dist = 0; + Anchor::label_class_t target = aln.get_label_classes().back(); + auto *anno_buffer = make_aln_graph(target).get_anno_buffer(); + align_fwd( + query_.get_graph(), + config_, aln.get_path().back(), num_matches, + query_window, std::numeric_limits::max(), + [&](auto&& spelling, auto&& path, auto&& cigar, size_t, Anchor::label_class_t final_target) { + size_t new_end_trim = aln.get_end_trim() - std::min(aln.get_end_trim(), spelling.size()); + + assert(path.empty() == (spelling.size() <= aln.get_end_trim())); + + std::string new_spelling = spelling.size() > aln.get_end_trim() + ? std::string(aln.get_spelling()) + spelling + : aln.get_path_spelling(); + + auto ext_path = aln.get_path(); + ext_path.insert(ext_path.end(), path.begin(), path.end()); + + size_t query_dist = cigar.get_num_query(); + assert(query_dist <= query_window.size()); + cigar.append(Cigar::CLIPPED, query_window.size() - query_dist); + + Cigar ext_cigar = aln.get_cigar(); + ext_cigar.trim_end_clipping(); + ext_cigar.append(std::move(cigar)); + + fwd_exts.emplace_back( + query_.get_graph(), + aln.get_query(), + aln.get_orientation(), + std::move(ext_path), + std::move(ext_cigar), + std::move(new_spelling), + new_end_trim, + final_target + ); + + set_coords(fwd_exts.back()); + + + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // start backtrack + size_t dist = std::get<0>(data); + if (dist == 0 && query_dist == 0) + return false; + + if (query_dist < max_query_dist) + return false; + + // Entries the node budget left unexplored can score higher + // than the best explored one; do not start from them. + auto score = get_score(cost, dist, query_dist); + return score == best_score; + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // terminate branch + if (query_dist > query_window.size()) + return true; + + size_t dist = std::get<0>(data); + auto score = get_score(cost, dist, query_dist); + return score <= 0 || config_.xdrop < best_score - score; + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // terminate + size_t dist = std::get<0>(data); + max_dist = std::max(dist, max_dist); + auto score = get_score(cost, dist, query_dist); + if (score >= best_score) { + best_score = score; + max_query_dist = std::max(max_query_dist, query_dist); + } + return query_dist == query_window.size(); + }, + aln.get_trim_spelling(), + anno_buffer, + target, + anno_buffer ? anno_buffer->get_label_path_coords({ aln.get_path().back() }, + { target }) + : AnnotationBuffer::CoordinateSet() + ); + + if (fwd_exts.empty()) + fwd_exts.emplace_back(aln); + } + + for (auto &fwd_ext : fwd_exts) { + fwd_ext.trim_end(); + if (!fwd_ext.get_clipping()) { + set_coords(fwd_ext); + callback(std::move(fwd_ext)); + continue; + } + + std::string_view query_window = fwd_ext.get_query(); + query_window.remove_suffix(fwd_ext.get_end_clipping() + fwd_ext.get_seed().size()); + + DBGAlignerConfig::score_t fwd_ext_score = score_match(fwd_ext, config_); + auto get_score = [&](size_t cost, size_t dist, size_t query_dist) { + return fwd_ext_score + cost_to_score(cost, query_dist, dist, match_score) + + (query_dist == query_window.size() ? config_.left_end_bonus : 0); + }; + + DBGAlignerConfig::score_t best_score = fwd_ext_score; + + auto it = aln.get_cigar().data().begin(); + auto end = aln.get_cigar().data().end(); + std::ignore = end; + assert(it != end); + if (it->first == Cigar::CLIPPED) { + ++it; + assert(it != end); + } + size_t num_matches = it->first == Cigar::MATCH ? it->second : 0; + + bool found = false; + size_t max_dist = 0; + size_t max_query_dist = 0; + auto *anno_buffer = make_aln_graph(fwd_ext.get_label_classes()[0]).get_anno_buffer(); + align_bwd( + query_.get_graph(), + config_, fwd_ext.get_path()[0], num_matches, + query_window, std::numeric_limits::max(), + [&](auto&& spelling, auto&& path, auto&& cigar, size_t, Anchor::label_class_t final_target) { + path.insert(path.end(), fwd_ext.get_path().begin(), fwd_ext.get_path().end()); + size_t query_dist = cigar.get_num_query(); + assert(query_dist <= query_window.size()); + Cigar ext_cigar(Cigar::CLIPPED, query_window.size() - query_dist); + ext_cigar.append(std::move(cigar)); + + Cigar aln_cigar = fwd_ext.get_cigar(); + aln_cigar.trim_clipping(); + ext_cigar.append(std::move(aln_cigar)); + + Alignment next_aln(query_.get_graph(), + fwd_ext.get_query(), + fwd_ext.get_orientation(), + std::move(path), + std::move(ext_cigar), + spelling + fwd_ext.get_path_spelling(), + fwd_ext.get_end_trim(), + final_target + ); + next_aln.trim_end(); + + set_coords(next_aln); + callback(std::move(next_aln)); + found = true; + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // start backtrack + size_t dist = std::get<0>(data); + + if (dist == 0 && query_dist == 0) + return false; + + if (query_dist < max_query_dist) + return false; + + // Entries the node budget left unexplored can score higher + // than the best explored one; do not start from them. + auto score = get_score(cost, dist, query_dist); + return score == best_score; + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // terminate branch + if (query_dist > query_window.size()) + return true; + + size_t dist = std::get<0>(data); + auto score = get_score(cost, dist, query_dist); + return score <= 0 || config_.xdrop < best_score - score; + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // terminate + size_t dist = std::get<0>(data); + max_dist = std::max(dist, max_dist); + auto score = get_score(cost, dist, query_dist); + if (score >= best_score) { + best_score = score; + max_query_dist = std::max(max_query_dist, query_dist); + } + return query_dist == query_window.size(); + }, + anno_buffer, + fwd_ext.get_label_classes()[0], + anno_buffer ? anno_buffer->get_label_path_coords({ fwd_ext.get_path()[0] }, + { fwd_ext.get_label_classes()[0] }) + : AnnotationBuffer::CoordinateSet() + ); + + if (!found) { + set_coords(fwd_ext); + callback(std::move(fwd_ext)); + } + } +} + +std::vector ExactSeeder::get_inexact_anchors(bool align) const { + common::logger->trace("Seeding"); + std::vector anchors = get_anchors(); + + const auto &graph = query_.get_graph(); + const auto *sshash = dynamic_cast(&graph); + + auto full_it = std::partition(anchors.begin(), anchors.end(), [](const auto &a) { return a.get_end_trim(); }); + std::sort(full_it, anchors.end(), AnchorLess()); + if (anchors.end() - full_it > 1) { + common::logger->trace("Merging {} anchors", anchors.size()); + assert(full_it->get_spelling().size() + full_it->get_end_trim() == graph.get_k()); + for (auto it = full_it; it + 1 != anchors.end(); ++it) { + auto &a_i = *it; + auto &a_j = *(it + 1); + assert(!a_i.get_end_trim()); + assert(!a_j.get_end_trim()); + assert(a_j.get_spelling().size() + a_j.get_end_trim() == graph.get_k()); + if (a_i.get_orientation() != a_j.get_orientation() || a_i.get_label_class() != a_j.get_label_class()) + continue; + + if (a_i.get_seed().size() < config_.max_seed_length + && a_i.get_end_clipping() == a_j.get_end_clipping() + 1 + && graph.has_single_outgoing(a_i.get_path().back()) + && graph.has_single_incoming(a_j.get_path()[0])) { + // merge + a_i.append(a_j, &graph); + std::swap(a_i, a_j); + a_i = Anchor(); + } + } + + full_it = std::remove_if(full_it, anchors.end(), [](const auto &a) { return a.get_path().empty(); }); + common::logger->trace("Merged {} anchors down to {}", anchors.size(), full_it - anchors.begin()); + anchors.erase(full_it, anchors.end()); + } + + using AnchorIt = std::vector::iterator; + std::sort(anchors.begin(), anchors.end(), AnchorLess()); + + + std::vector alignments; + + if (anchors.empty()) + return alignments; + + ssize_t k = graph.get_k(); + + assert(!dynamic_cast(this) || std::all_of(anchors.begin(), anchors.end(), + [](const auto &a) { return a.get_label_class() != Anchor::nannot; })); + + + bool first_chain = true; + DBGAlignerConfig::score_t first_chain_score = DBGAlignerConfig::ninf; + size_t num_chains = 0; + + std::string dummy(query_.get_query().size(), '$'); + + auto *anno_buffer = make_aln_graph(anchors[0].get_label_class()).get_anno_buffer(); + + bool ext_success = false; + Anchor::score_t last_chain_score = 0; + common::logger->trace("Chaining"); + chain_anchors(query_, config_, anchors.begin(), anchors.end(), + [this,k,sshash,match_score=config_.match_score("A"),&dummy]( + const Anchor &a_j, + ssize_t max_dist, + AnchorIt rbegin, + AnchorIt rend, + typename ChainScores::iterator rchain_scores, + const ScoreUpdater &update_score) { + auto rchain_scores_end = rchain_scores + (rend - rbegin); + + auto begin = std::make_reverse_iterator(rend); + auto end = std::make_reverse_iterator(rbegin); + auto chain_scores = std::make_reverse_iterator(rchain_scores_end); + + std::string_view query_j = a_j.get_seed(); + const DBGAlignerConfig::score_t &score_j = std::get<0>(*(rchain_scores_end)); + + auto try_chain = [&](auto it, auto chain_scores) { + const Anchor &a_i = *it; + assert(a_i.get_label_class() == a_j.get_label_class()); + assert(a_i.get_orientation() == a_j.get_orientation()); + + if (a_j.get_clipping() + k > query_.get_query(a_i.get_orientation()).size()) + return; + + auto get_sshash_distance = [&](DeBruijnGraph::node_index node, DeBruijnGraph::node_index target_node) -> size_t { + uint64_t contig = sshash->get_sequence_id(node); + uint64_t target_contig = sshash->get_sequence_id(target_node); + if (contig == target_contig) { + if (node < sshash->reverse_complement(node)) { + // canonical strand + if (target_node >= node) { +#ifndef NDEBUG + if (target_node != node) { + std::vector path(target_node - node + 1); + std::iota(path.begin(), path.end(), node); + assert(path.back() == target_node); + assert(spell_path(*sshash, path) != ""); + } +#endif + return target_node - node; + } + } else if (node >= target_node && target_node > sshash->reverse_complement(target_node) && node > sshash->reverse_complement(node)) { + // rc strand +#ifndef NDEBUG + if (target_node != node) { + std::vector path(node - target_node + 1); + std::iota(path.begin(), path.end(), target_node); + assert(path.back() == node); + std::reverse(path.begin(), path.end()); + assert(spell_path(*sshash, path) != ""); + } +#endif + return node - target_node; + } + } + + return std::numeric_limits::max(); + }; + + auto graph = make_aln_graph(a_i.get_label_class()); + + std::string_view query_i = a_i.get_seed(); + + DBGAlignerConfig::score_t dist = query_j.end() - query_i.end(); + + if (dist < 0 || query_i.begin() >= query_j.begin()) + return; + + DBGAlignerConfig::score_t base_score = std::get<0>(*chain_scores); + + if (dist == 0 && a_j.get_end_trim()) { + assert(a_j.get_path().size() == 1); + if (base_score <= score_j) + return; + + if (a_i.get_coord() != Anchor::ncoord && a_j.get_coord() != Anchor::ncoord) { + size_t a_i_pos = a_i.get_coord() + a_i.get_seed().size(); + size_t a_j_pos = a_j.get_coord() + a_j.get_seed().size(); + if (a_i_pos == a_j_pos) + update_score(base_score, it, 0); + + return; + } + + std::string_view a_j_trim_spelling = a_j.get_trim_spelling(); + if (a_i.get_end_trim()) { + std::string_view a_i_trim_spelling = a_i.get_trim_spelling(); + if (!std::equal(a_i_trim_spelling.begin(), a_i_trim_spelling.end(), + a_j_trim_spelling.begin(), a_j_trim_spelling.begin() + a_i_trim_spelling.size())) + return; + a_j_trim_spelling = a_j_trim_spelling.substr(a_i_trim_spelling.size()); + } + DeBruijnGraph::node_index node = a_i.get_path().back(); + size_t traverse = sshash + ? get_sshash_distance(node, a_j.get_path()[0]) + : std::numeric_limits::max(); + if (traverse == std::numeric_limits::max() || traverse != a_j_trim_spelling.size()) { + traverse = 0; + graph.traverse(a_i.get_path().back(), a_j_trim_spelling, + [&](DeBruijnGraph::node_index next, Anchor::label_class_t next_target) { + assert(next_target == a_i.get_label_class()); + ++traverse; + node = next; + }, + []() { return false; }, + a_j_trim_spelling.size()); + assert(traverse != a_j_trim_spelling.size() || node == a_j.get_path()[0]); + } else { + node = a_j.get_path()[0]; + } + if (traverse == a_j_trim_spelling.size()) { + update_score(base_score, it, 0); + } + } + + auto added_seq_start = std::max(query_j.begin(), query_i.end()); + std::string_view added_seq(added_seq_start, query_j.end() - added_seq_start); + + DBGAlignerConfig::score_t score = base_score + config_.match_score(added_seq) + + (config_.right_end_bonus * !a_j.get_end_clipping()); + + if (score <= score_j) + return; + + if (a_i.get_coord() != Anchor::ncoord && a_j.get_coord() != Anchor::ncoord) { + size_t a_i_pos = a_i.get_coord() + a_i.get_seed().size(); + size_t a_j_pos = a_j.get_coord() + a_j.get_seed().size(); + if (a_j_pos > a_i_pos) { + size_t traversed = a_j_pos - a_i_pos; + size_t gap = std::abs(static_cast(traversed) - static_cast(dist)); + if (gap) { + score -= (0.01 * config_.min_seed_length * gap + log2l(static_cast(gap)) / 2) * match_score; + } + + if (score > score_j) + update_score(score, it, traversed); + } + + return; + } + + size_t i_seed_size = k - a_i.get_end_trim(); + auto last_node_i_query = query_i.end() - i_seed_size; + + if (last_node_i_query >= query_j.begin()) { + // there should be overlapping nodes + size_t j_start = last_node_i_query - query_j.begin(); + if (a_i.get_path().size() >= j_start + && a_j.get_path().size() >= j_start + && std::equal(a_j.get_path().begin(), a_j.get_path().begin() + j_start, + a_i.get_path().end() - j_start, a_i.get_path().end())) { + update_score(score, it, dist); + } + return; + } + + if (query_i.end() > query_j.begin()) { + // overlapping nucleotides + assert(query_j.data() > query_i.data()); + if (!a_i.get_end_trim() && !a_j.get_end_trim() && query_i.end() - query_j.begin() == k - 1) { + // no need to check + update_score(score, it, dist); + } + return; + } + + if (sshash) { + size_t traversed = get_sshash_distance(a_i.get_path().back(), a_j.get_path()[0]); + if (traversed > 0 + && traversed != std::numeric_limits::max() + && traversed + a_i.get_path().size() + query_j.size() > query_i.size() + 1) { + traversed += a_i.get_path().size() - 1 - query_i.size() + query_j.size(); + size_t gap = std::abs(static_cast(traversed) - static_cast(dist)); + if (gap) { + score += config_.gap_opening_penalty + (gap - 1) * config_.gap_extension_penalty; + } + size_t nmismatch = query_j.data() - (query_i.data() + query_i.size()); + + // don't double-penalize insertions + if (dist > static_cast(traversed) && gap <= nmismatch) + nmismatch -= gap; + + std::string_view query_trim_window(query_i.data() + query_i.size(), nmismatch); + auto cur_mismatch = config_.score_sequences(query_trim_window, + std::string_view(dummy.c_str(), nmismatch)); + + + score += cur_mismatch; + + if (score > score_j) + update_score(score, it, traversed); + } + } + + if (a_i.get_end_trim() && a_i.get_clipping() + a_i.get_seed().size() + a_i.get_end_trim() >= a_j.get_clipping()) { + // overlap in end parts + std::string_view a_i_trim = a_i.get_trim_spelling(); + + size_t olap = std::min(a_i.get_query().size(), + a_i.get_clipping() + a_i.get_seed().size() + a_i.get_end_trim()) - a_j.get_clipping(); + auto a_j_spelling = a_j.get_path_spelling().substr(0, k); + assert(a_i_trim.size() >= olap); + for (size_t del = 0; del < olap; ++del) { + assert(a_i_trim.size() >= olap - del); + + // if the character before a_j also matches, then skip this + // since we should work with an earlier seed instead + assert(a_j.get_clipping()); + if (a_i_trim.size() >= olap - del + 1 && *(a_i_trim.end() - olap + del - 1) == *(a_j.get_seed().data() - 1)) + continue; + + if (std::equal(a_i_trim.end() - olap + del, a_i_trim.end(), a_j_spelling.begin(), a_j_spelling.begin() + olap - del)) { + DeBruijnGraph::node_index node = a_i.get_path().back(); + size_t traversed = sshash + ? get_sshash_distance(node, a_j.get_path()[0]) + : std::numeric_limits::max(); + + if (traversed == std::numeric_limits::max() || traversed != a_j_spelling.size() - olap + del) { + traversed = 0; + graph.traverse(node, a_j_spelling.substr(olap - del), + [&](DeBruijnGraph::node_index next, Anchor::label_class_t next_target) { + assert(next_target == a_i.get_label_class()); + node = next; + ++traversed; + }, + []() { return false; }, + a_j_spelling.size() - olap + del + ); + assert(traversed != a_j_spelling.size() - olap + del || node == a_j.get_path()[0]); + } else { + node = a_j.get_path()[0]; + } + + if (traversed == a_j_spelling.size() - olap + del) { + traversed += a_i.get_path().size() - 1; + assert(traversed >= del); + assert(static_cast(traversed - del) == query_j.begin() - query_i.begin()); + DBGAlignerConfig::score_t cur_score = score; + size_t nmismatch = query_j.begin() - query_i.end(); + std::string_view query_trim_window(query_i.data() + query_i.size(), nmismatch); + auto cur_mismatch = config_.score_sequences(query_trim_window, + del == 0 ? std::string_view(a_i_trim.begin(), nmismatch) + : std::string_view(dummy.c_str(), nmismatch)); + cur_score += cur_mismatch; + if (del) + cur_score += config_.gap_opening_penalty + (del - 1) * config_.gap_extension_penalty; + if (cur_score > score_j) { + update_score(cur_score, it, traversed - query_i.size() + query_j.size()); + break; + } + } + } + } + } + + if (!a_i.get_end_trim() && query_i.end() != query_j.begin()) { + // detect insertions + + // don't try inserting if the previous character matches + if (!graph.traverse_back(a_j.get_path()[0], *(query_j.data() - 1)).first) { + DeBruijnGraph::node_index node = a_i.get_path().back(); + size_t traversed = sshash + ? get_sshash_distance(node, a_j.get_path()[0]) + : std::numeric_limits::max(); + if (traversed == std::numeric_limits::max() || traversed != query_i.size() - a_i.get_path().size() + 1) { + traversed = 0; + auto a_j_spelling = a_j.get_path_spelling().substr(0, k); + graph.traverse(node, a_j_spelling, + [&](DeBruijnGraph::node_index next, Anchor::label_class_t next_target) { + assert(next_target == a_i.get_label_class()); + node = next; + ++traversed; + }, + []() { return false; }, + query_i.size() - a_i.get_path().size() + 1 + ); + assert(traversed != query_i.size() - a_i.get_path().size() + 1 || node == a_j.get_path()[0]); + } else { + node = a_j.get_path()[0]; + } + if (traversed == query_i.size() - a_i.get_path().size() + 1) { + size_t query_dist = a_j.get_clipping() - (a_i.get_clipping() + a_i.get_path().size() - 1); + assert(query_dist >= traversed); + size_t gap = query_dist - traversed; + if (gap > 0) + score += config_.gap_opening_penalty + (gap - 1) * config_.gap_extension_penalty; + + if (score > score_j) { + update_score(score, it, traversed - (query_i.size() - a_i.get_path().size() + 1) + query_j.size()); + } + } + } + } + + // TODO: if a_i.get_end_trim() != 0, then we need to deal with + // a combination of D from the trim, X to the trim, then I of the query + }; + + + for (auto rit = begin; rit != end; ++rit, ++chain_scores) { + try_chain(rit.base() - 1, chain_scores); + } + + + }, + [&](const AnchorChain &chain, const std::vector &score_traceback) { + assert(chain.size()); + + ext_success = false; + + bool ret_val = (!chain[0].first->get_clipping() && !chain[0].first->get_end_clipping()) + || chain.size() > 1 + || std::any_of(chain.begin(), chain.end(), + [this](const auto &a) { return a.first->get_seed().size() > config_.min_seed_length + || a.first->get_path().size() > 1; }); + + + + if (!ret_val) + return false; + + last_chain_score = score_traceback.back(); + ret_val = score_traceback.back() >= first_chain_score * config_.rel_score_cutoff; + + if (ret_val) { + common::logger->trace("Chain\t{}\t{}", score_traceback.back(), ret_val); + } + + + return ret_val; + }, + [this,align,anno_buffer,match_score=config_.match_score("A")]( + AnchorIt last, + AnchorIt next, + Alignment&& aln, + size_t next_to_last_dist, + DBGAlignerConfig::score_t score_up_to_now, + const AlignmentCallback &callback) { + size_t traversal_dist = next_to_last_dist - last->get_seed().size() + next->get_seed().size(); + const auto &graph = query_.get_graph(); + if (next->get_seed().begin() + next->get_path().size() == last->get_seed().begin() && traversal_dist == next->get_path().size()) { + // simple extension +#ifndef NDEBUG + size_t k = graph.get_k(); + if (graph.traverse(next->get_path().back(), aln.get_path_spelling()[k-1]) != aln.get_path()[0]) { + std::cerr << "\t" << next->get_path().back() << " > " << aln.get_path_spelling()[k-1] << "\t" << graph.traverse(next->get_path()[0], aln.get_path_spelling()[k-1]) << " vs. " << aln.get_path()[0] << std::endl; + assert(false); + } +#endif + std::vector path; + path.resize(aln.get_path().size() + next->get_path().size()); + std::copy(next->get_path().begin(), next->get_path().end(), path.begin()); + std::copy(aln.get_path().begin(), aln.get_path().end(), path.begin() + next->get_path().size()); + + Cigar cigar(Cigar::CLIPPED, next->get_clipping()); + + Cigar aln_cigar = aln.get_cigar(); + aln_cigar.trim_clipping(); + + cigar.append(Cigar::MATCH, next->get_path().size()); + cigar.append(std::move(aln_cigar)); + + Alignment next_aln(graph, + aln.get_query(), + aln.get_orientation(), + std::move(path), + std::move(cigar), + std::string(next->get_seed().substr(0, next->get_path().size())) + aln.get_path_spelling(), + aln.get_end_trim(), + next->get_label_class()); + callback(std::move(next_aln)); + return; + } + + if (!align) { + // connect anchors without actually aligning or traversing + size_t query_dist = last->get_clipping() - next->get_clipping(); + std::vector path; + path.resize(aln.get_path().size() + traversal_dist); + std::copy(aln.get_path().rbegin(), aln.get_path().rend(), path.rbegin()); + + Cigar cigar(Cigar::CLIPPED, next->get_clipping()); + std::string spelling; + + if (next->get_seed().end() > last->get_seed().begin()) { + // overlap + assert(traversal_dist == query_dist); + // TODO: don't bother traversing, trust that the connection worked. + // keep intermediate nodes as DeBruijnGraph::npos +#ifndef NDEBUG + DeBruijnGraph::node_index node = *(path.rbegin() + aln.get_path().size() - 1); + for (auto it = next->get_seed().rend() - query_dist; it != next->get_seed().rend(); ++it) { + node = graph.traverse_back(node, *it); + assert(graph.in_graph(node)); + } + assert(node == next->get_path()[0]); +#endif + cigar.append(Cigar::MATCH, query_dist); + assert(next->get_query().begin() + next->get_clipping() + query_dist <= next->get_query().end()); + spelling += std::string(next->get_seed().data(), next->get_seed().data() + query_dist); + } else { + cigar.append(Cigar::MATCH, next->get_seed().size()); + spelling += next->get_seed(); + if (traversal_dist == query_dist) { + // no gap + assert(query_dist > next->get_seed().size()); + cigar.append(Cigar::MISMATCH, query_dist - next->get_seed().size()); + spelling += std::string(query_dist - next->get_seed().size(), '$'); + } else if (traversal_dist > query_dist) { + // deletion + // TODO: which first, mismatch or deletion? + cigar.append(Cigar::DELETION, traversal_dist - query_dist); + spelling += std::string(traversal_dist - query_dist, '$'); + if (query_dist > next->get_seed().size()) { + cigar.append(Cigar::MISMATCH, query_dist - next->get_seed().size()); + spelling += std::string(query_dist - next->get_seed().size(), '$'); + } + } else { + // insertion + // TODO: which first, mismatch or insert? + cigar.append(Cigar::INSERTION, query_dist - traversal_dist); + if (traversal_dist > next->get_seed().size()) { + cigar.append(Cigar::MISMATCH, traversal_dist - next->get_seed().size()); + spelling += std::string(traversal_dist - next->get_seed().size(), '$'); + } + } + + if (next->get_end_trim()) { + assert(next->get_path().size() == 1); + std::string_view next_trim_spelling = next->get_trim_spelling(); + auto end = next_trim_spelling.begin() + std::min(next_trim_spelling.size(), spelling.size() - next->get_seed().size()); + assert(next->get_seed().size() + (end - next_trim_spelling.begin()) <= spelling.size()); + std::copy(next_trim_spelling.begin(), end, + spelling.begin() + next->get_seed().size()); + } + } + + std::copy(next->get_path().begin(), next->get_path().end(), path.begin()); + + Cigar aln_cigar = aln.get_cigar(); + aln_cigar.trim_clipping(); + cigar.append(std::move(aln_cigar)); + + Alignment next_aln(graph, + aln.get_query(), + aln.get_orientation(), + std::move(path), + std::move(cigar), + spelling + aln.get_path_spelling(), + aln.get_end_trim(), + next->get_label_class()); + callback(std::move(next_aln)); + return; + } + + DeBruijnGraph::node_index target_node = next->get_path().back(); + assert(traversal_dist > next->get_path().size() - 1); + traversal_dist -= next->get_path().size() - 1; + + std::string_view query_window( + next->get_seed().data() + next->get_path().size() - 1, + aln.get_seed().begin() - next->get_seed().begin() - (next->get_path().size() - 1) + ); + + bool aln_found = false; + align_bwd( + query_.get_graph(), + config_, + aln.get_path()[0], + last->get_seed().size(), + query_window, + traversal_dist, + [&](auto&& spelling, auto&& bt_path, auto&& bt_cigar, size_t cost, Anchor::label_class_t final_target) { // callback + assert(bt_path[0] == next->get_path().back()); + bt_path.insert(bt_path.begin(), next->get_path().begin(), next->get_path().end() - 1); + + Cigar cigar(Cigar::CLIPPED, next->get_clipping()); + cigar.append(Cigar::MATCH, next->get_path().size() - 1); + cigar.append(std::move(bt_cigar)); + + Cigar aln_cigar = aln.get_cigar(); + aln_cigar.trim_clipping(); + cigar.append(std::move(aln_cigar)); + bt_path.insert(bt_path.end(), aln.get_path().begin(), aln.get_path().end()); + + Alignment next_aln(query_.get_graph(), + aln.get_query(), + aln.get_orientation(), + std::move(bt_path), + std::move(cigar), + std::string(next->get_spelling().substr(0, next->get_path().size() - 1)) + spelling + aln.get_path_spelling(), + aln.get_end_trim(), + final_target); + aln_found = true; + callback(std::move(next_aln)); + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // start_backtracking + const auto &[dist, last_ext, last_node, last_op, mismatch_char, num_matches, cur_target] = data; + + if (aln_found || dist == 0 || query_dist == 0) + return false; + + bool ret_val = query_dist == query_window.size() && node == target_node; + if (ret_val) + common::logger->trace("Final cost: {}", cost); + + return ret_val; + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // terminate branch + const auto &[dist, last_ext, last_node, last_op, mismatch_char, num_matches, cur_target] = data; + + return dist > traversal_dist || query_dist > query_window.size(); + }, + [&](size_t cost, const SMap &data, size_t query_dist, DeBruijnGraph::node_index node) { + // terminate + std::ignore = cost; + std::ignore = data; + std::ignore = query_dist; + return aln_found || (query_dist == query_window.size() && node == target_node); + }, + anno_buffer, + last->get_label_class(), + anno_buffer ? anno_buffer->get_label_path_coords({ aln.get_path()[0] }, + { last->get_label_class() }) + : AnnotationBuffer::CoordinateSet() + ); + // The node budget can stop the gap fill before it reaches the + // previous anchor. The chain then yields no alignment. + if (!aln_found) + common::logger->trace("Failed to connect anchors at query positions {} and {}", + next->get_clipping(), aln.get_clipping()); + }, + [this,&num_chains,&alignments,&ext_success,&last_chain_score,&first_chain,&first_chain_score](Alignment&& aln) { + aln.trim_end(); + + assert(aln.get_label_classes().size()); + if (auto anno_buffer = make_aln_graph(aln.get_label_classes()[0]).get_anno_buffer()) { + if (anno_buffer->has_coordinates()) { + auto coords = anno_buffer->get_label_path_coords(aln.get_path(), aln.get_label_classes()); + assert(coords.size()); + aln.set_coords(std::move(coords)); + } + } + + if (!ext_success) { + ++num_chains; + ext_success = true; + if (first_chain) { + first_chain_score = last_chain_score; + first_chain = false; + } + } + + if (common::get_verbose()) { + common::logger->trace("Aln: {}\t{}\t{}\t{}", + aln.get_orientation() ? "-" : "+", + score_match(aln, config_), aln.get_cigar().to_string(), aln.get_label_classes()[0]); + } + alignments.emplace_back(std::move(aln)); + }, + []() { return false; }, + anno_buffer + ); + + + return alignments; +} + +AlignmentGraph ExactSeeder::make_aln_graph(Anchor::label_class_t) const { + return AlignmentGraph(query_.get_graph()); +} + +AlignmentGraph LabeledSeeder::make_aln_graph(Anchor::label_class_t target) const { + return AlignmentGraph(query_.get_graph(), &anno_buffer_, target); +} + +AlignmentGraph Extender::make_aln_graph(Anchor::label_class_t) const { + return AlignmentGraph(query_.get_graph()); +} + +AlignmentGraph LabeledExtender::make_aln_graph(Anchor::label_class_t target) const { + return AlignmentGraph(query_.get_graph(), &anno_buffer_, target); +} + + +} // namespace mtg::graph::align diff --git a/metagraph/src/graph/alignment_redone/aln_seeder.hpp b/metagraph/src/graph/alignment_redone/aln_seeder.hpp new file mode 100644 index 0000000000..d9013bc097 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/aln_seeder.hpp @@ -0,0 +1,102 @@ +#pragma once + +#include + +#include "aln_match.hpp" +#include "aln_query.hpp" +#include "aligner_config.hpp" +#include "annotation_buffer.hpp" + +namespace mtg::graph::align_redone { + +class AlignmentGraph; + +class Seeder { + public: + Seeder(const Query &query, const DBGAlignerConfig &config) + : query_(query), config_(config) {} + + virtual ~Seeder() {} + + virtual std::vector get_anchors() const = 0; + virtual std::vector get_inexact_anchors(bool align = true) const = 0; + + const Query& get_query() const { return query_; } + + protected: + const Query &query_; + const DBGAlignerConfig &config_; +}; + +class ExactSeeder : public Seeder { + public: + template + ExactSeeder(Args&&... args) : Seeder(std::forward(args)...) {} + + virtual ~ExactSeeder() {} + + std::vector get_anchors() const override; + std::vector get_inexact_anchors(bool align = true) const override; + + protected: + virtual AlignmentGraph make_aln_graph(Anchor::label_class_t) const; +}; + +class LabeledSeeder : public ExactSeeder { + public: + template + LabeledSeeder(AnnotationBuffer &anno_buffer, Args&&... args) + : ExactSeeder(std::forward(args)...), + anno_buffer_(anno_buffer) {} + + std::vector get_anchors() const override; + + protected: + virtual AlignmentGraph make_aln_graph(Anchor::label_class_t target) const override; + + private: + AnnotationBuffer &anno_buffer_; +}; + +class Extender { + public: + Extender(const Query &query, const DBGAlignerConfig &config) + : query_(query), config_(config) {} + + virtual ~Extender() {} + + void extend(const Alignment &aln, + const std::function &callback) const; + + protected: + const Query &query_; + const DBGAlignerConfig &config_; + + virtual AlignmentGraph make_aln_graph(Anchor::label_class_t) const; +}; + +class LabeledExtender : public Extender { + public: + template + LabeledExtender(AnnotationBuffer &anno_buffer, Args&&... args) + : Extender(std::forward(args)...), + anno_buffer_(anno_buffer) {} + + private: + AnnotationBuffer &anno_buffer_; + + virtual AlignmentGraph make_aln_graph(Anchor::label_class_t target) const override; +}; + +void align_query(const Query &query, + const Seeder &seeder, + const Extender &extender, + const std::function &callback, + bool connect_anchors_in_chain = true); + +void align_query(const Query &query, + const Seeder &seeder, + const std::function &callback, + bool connect_anchors_in_chain = true); + +} // namespace mtg::graph::align \ No newline at end of file diff --git a/metagraph/src/graph/alignment_redone/annotation_buffer.cpp b/metagraph/src/graph/alignment_redone/annotation_buffer.cpp new file mode 100644 index 0000000000..d260fc8e19 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/annotation_buffer.cpp @@ -0,0 +1,281 @@ +#include "annotation_buffer.hpp" +#include "graph/annotation_key_mapper.hpp" + +#include + +#include "aln_match.hpp" + +#include "graph/representation/canonical_dbg.hpp" +#include "annotation/binary_matrix/base/binary_matrix.hpp" +#include "common/utils/template_utils.hpp" +#include "common/algorithms.hpp" + +namespace mtg { +namespace graph { +namespace align_redone { + +using mtg::common::logger; + +typedef annot::matrix::BinaryMatrix::Row Row; +typedef annot::matrix::BinaryMatrix::Column Column; + +AnnotationBuffer::AnnotationBuffer(const DeBruijnGraph &graph, const Annotator &annotator) + : graph_(graph), + annotator_(annotator), + multi_int_(dynamic_cast(&annotator_.get_matrix())), + canonical_(dynamic_cast(&graph_)), + column_sets_({ {} }) { + if (multi_int_ && graph_.get_mode() != DeBruijnGraph::BASIC) { + multi_int_ = nullptr; + logger->warn("Coordinates not supported when aligning to CANONICAL " + "or PRIMARY mode graphs"); + } +} + +void AnnotationBuffer::fetch_queued_annotations() { + assert(graph_.get_mode() != DeBruijnGraph::PRIMARY + && "PRIMARY graphs must be wrapped into CANONICAL"); + + std::vector queued_nodes; + std::vector queued_rows; + + const AnnotationKeyMapper key_mapper(graph_); + + for (const auto &path : queued_paths_) { + assert(std::find(path.begin(), path.end(), DeBruijnGraph::npos) == path.end()); + assert(path == map_to_nodes_sequentially(graph_, spell_path(graph_, path))); + const std::string sequence = key_mapper.needs_sequence() ? spell_path(graph_, path) : ""; + auto base_path = key_mapper.keys(path, sequence); + assert(base_path.size() == path.size()); + assert(std::find(base_path.begin(), base_path.end(), DeBruijnGraph::npos) + == base_path.end()); + + for (size_t i = 0; i < path.size(); ++i) { + Row row = AnnotatedDBG::graph_to_anno_index(base_path[i]); + auto [it, inserted] = node_to_cols_.try_emplace(base_path[i], nannot); + if (inserted) { + assert(path[i] == base_path[i] || !node_to_cols_.count(path[i])); + if (path[i] != base_path[i]) + node_to_cols_.emplace(path[i], nannot); + + queued_rows.push_back(row); + queued_nodes.emplace_back(path[i]); + if (has_coordinates()) { + label_coords_.emplace_back(); + if (path[i] != base_path[i]) + label_coords_.emplace_back(); + } + } else if (base_path[i] != path[i] && it->second != nannot) { + auto [jt, jnserted] = node_to_cols_.try_emplace(path[i], it->second); + if (jnserted && has_coordinates()) + label_coords_.emplace_back(label_coords_[it - node_to_cols_.begin()]); + } + } + } + + queued_paths_.clear(); + + if (queued_nodes.empty()) + return; + + auto push_node_labels = [&](auto node_it, auto row_it, auto&& labels, CoordinateSet *coords) { + assert(node_it != queued_nodes.end()); + assert(row_it != queued_rows.end()); + node_index node = *node_it; + node_index base_node = AnnotatedDBG::anno_to_graph_index(*row_it); + + assert(node_to_cols_.count(node)); + assert(node_to_cols_.count(base_node)); + + assert(node_to_cols_[node] == nannot); + assert(node_to_cols_[base_node] == nannot); + + label_class_t label_i = cache_column_set(std::move(labels)); + auto it = node_to_cols_.find(base_node); + it.value() = label_i; + if (has_coordinates()) { + assert(coords); + label_coords_[it - node_to_cols_.begin()] = *coords; + } + + if (node != base_node) { + auto jt = node_to_cols_.find(node); + jt.value() = label_i; + if (has_coordinates()) { + assert(coords); + label_coords_[jt - node_to_cols_.begin()] = std::move(*coords); + } + } + }; + + auto node_it = queued_nodes.begin(); + auto row_it = queued_rows.begin(); + if (has_coordinates()) { + assert(multi_int_); + // extract both labels and coordinates, then store them separately + for (auto&& row_tuples : multi_int_->get_row_tuples(queued_rows)) { + std::sort(row_tuples.begin(), row_tuples.end(), utils::LessFirst()); + Columns labels; + labels.reserve(row_tuples.size()); + + CoordinateSet all_coords; + all_coords.reserve(row_tuples.size()); + + for (auto&& [label, coords] : row_tuples) { + labels.push_back(label); + all_coords.emplace_back(coords.begin(), coords.end()); + } + push_node_labels(node_it++, row_it++, std::move(labels), &all_coords); + } + } else { + for (auto&& labels : annotator_.get_matrix().get_rows(queued_rows)) { + std::sort(labels.begin(), labels.end()); + push_node_labels(node_it++, row_it++, std::move(labels), nullptr); + } + } + +#ifndef NDEBUG + for (const auto &[node, val] : node_to_cols_) { + assert(val != nannot); + } +#endif +} + +auto AnnotationBuffer::get_labels_id_and_coords(node_index node) const + -> std::pair { + std::pair ret_val { nannot, nullptr }; + + if (canonical_) + node = canonical_->get_base_node(node); + + auto it = node_to_cols_.find(node); + + // if the node hasn't been seen before, or if its annotations haven't + // been fetched, return nothing + if (it == node_to_cols_.end() || it->second == nannot) + return ret_val; + + ret_val.first = it->second; + + if (has_coordinates()) { + assert(static_cast(it - node_to_cols_.begin()) < label_coords_.size()); + ret_val.second = &label_coords_[it - node_to_cols_.begin()]; + } + + return ret_val; +} + +auto AnnotationBuffer::get_labels_and_coords(node_index node) const + -> std::pair { + auto [labels_i, coords] = get_labels_id_and_coords(node); + return std::make_pair( + labels_i < column_sets_.size() ? &column_sets_.data()[labels_i] : nullptr, + coords + ); +} + +bool AnnotationBuffer::node_has_labels(node_index node, label_class_t target) const { + assert(target != nannot); + + auto node_full_columns = get_labels(node); + assert(node_full_columns); + + const Columns &target_columns = get_cached_column_set(target); + return utils::count_intersection(node_full_columns->begin(), node_full_columns->end(), + target_columns.begin(), target_columns.end()) == target_columns.size(); +} + +auto AnnotationBuffer +::get_label_path_coords(const std::vector &path, + const std::vector &label_path) const + -> CoordinateSet { + assert(path.size() == label_path.size()); + + if (!has_coordinates() || path.empty()) + return {}; + + label_class_t source = label_path[0]; + assert(source); + assert(source != nannot); + assert(std::all_of(label_path.begin() + 1, + label_path.end(), + [&](label_class_t t) { return t == source; })); + assert(node_has_labels(path[0], source)); + + CoordinateSet ret_val; + + auto [source_full_columns, source_full_coords] = get_labels_and_coords(path[0]); + assert(source_full_columns); + assert(source_full_coords); + assert(source_full_coords->size() == source_full_columns->size()); + + auto [target_full_columns, target_full_coords] = get_labels_and_coords(path.back()); + assert(target_full_columns); + assert(target_full_coords); + assert(target_full_coords->size() == target_full_columns->size()); + + const Columns &source_columns = get_cached_column_set(source); + assert(source_columns.size()); + auto it_s = source_columns.begin(); + + size_t coord_dist = path.size() - 1; + utils::match_indexed_values( + source_full_columns->begin(), source_full_columns->end(), source_full_coords->begin(), + target_full_columns->begin(), target_full_columns->end(), target_full_coords->begin(), + [&](Column c, const Tuple &source_c, const Tuple &target_c) { + assert(source_c.size()); + assert(target_c.size()); + if (it_s != source_columns.end() && c == *it_s) { + Tuple &out = ret_val.emplace_back(); + utils::set_intersection(source_c.begin(), source_c.end(), + target_c.begin(), target_c.end(), + std::back_inserter(out), + coord_dist); + assert(out.size()); + ++it_s; + } + } + ); + + assert(it_s == source_columns.end()); + + return ret_val; +} + +std::string AnnotationBuffer::generate_column_set_str(label_class_t i, size_t spelling_size) const { + if (i == nannot) + return "*"; + + std::ostringstream ostd; + const Columns &columns = get_cached_column_set(i); + if (columns.empty()) + return "*"; + + assert(columns[0] != ncolumn); + assert(columns[0] < annotator_.get_label_encoder().size()); + ostd << annotator_.get_label_encoder().decode(columns[0]); + if (label_coords_.size()) { + assert(label_coords_.size() == node_to_cols_.size()); + for (auto coord : label_coords_[i][0]) { + ostd << ":" << coord + 1 << "-" << coord + spelling_size; + } + } + + for (Column c = 1; c < columns.size(); ++c) { + assert(columns[c] != ncolumn); + assert(columns[c] < annotator_.get_label_encoder().size()); + ostd << ";" << annotator_.get_label_encoder().decode(columns[c]); + if (label_coords_.size()) { + assert(label_coords_.size() == node_to_cols_.size()); + for (auto coord : label_coords_[i][c]) { + ostd << ":" << coord + 1 << "-" << coord + spelling_size; + } + } + } + + return ostd.str(); +} + +} // namespace align +} // namespace graph +} // namespace mtg diff --git a/metagraph/src/graph/alignment_redone/annotation_buffer.hpp b/metagraph/src/graph/alignment_redone/annotation_buffer.hpp new file mode 100644 index 0000000000..4ec5e37646 --- /dev/null +++ b/metagraph/src/graph/alignment_redone/annotation_buffer.hpp @@ -0,0 +1,116 @@ +#pragma once + +#include "graph/annotated_dbg.hpp" +#include "annotation/int_matrix/base/int_matrix.hpp" +#include "common/vector_set.hpp" +#include "common/vector_map.hpp" +#include "common/hashers/hash.hpp" + +namespace mtg::graph { + +class CanonicalDBG; + +namespace align_redone { + +// caches queried annotations to speed up next queries (labels with or w/o coordinates) +class AnnotationBuffer { + public: + using label_class_t = size_t; + using coord_t = int64_t; + + using Annotator = AnnotatedDBG::Annotator; + using node_index = DeBruijnGraph::node_index; + + using Column = annot::matrix::BinaryMatrix::Column; + using Tuple = SmallVector; + using Columns = annot::matrix::BinaryMatrix::SetBitPositions; + using CoordinateSet = Vector; + + // dummy index for an unfetched annotations + static constexpr label_class_t nannot = std::numeric_limits::max(); + static constexpr Column ncolumn = std::numeric_limits::max(); + static constexpr coord_t ncoord = std::numeric_limits::max(); + + AnnotationBuffer(const DeBruijnGraph &graph, const Annotator &annotator); + + void queue_path(const std::vector &path) { + queued_paths_.push_back(path); + } + + void queue_path(std::vector&& path) { + queued_paths_.push_back(std::move(path)); + } + + // fetch annotations for the queued nodes from the buffer and reset the buffer + void fetch_queued_annotations(); + + bool has_coordinates() const { return multi_int_; } + + // Get the annotations and coordinates of a node if they have been fetched. + // The returned pointers are valid until next fetch_queued_annotations(). + std::pair + get_labels_and_coords(node_index node) const; + + std::pair + get_labels_id_and_coords(node_index node) const; + + // get the labels of a node if they have been fetched + inline const Columns* get_labels(node_index node) const { + return get_labels_and_coords(node).first; + } + + // get the labels of a node if they have been fetched + inline label_class_t get_labels_id(node_index node) const { + return get_labels_id_and_coords(node).first; + } + + const Annotator& get_annotator() const { return annotator_; } + + size_t num_nodes_buffered() const { return node_to_cols_.size(); } + size_t num_column_sets() const { return column_sets_.size(); } + + // This method lets the caller push additional column sets to dictionary + // `column_sets_`. These column sets can later be fetched with `get_cached_column_set()` + template + inline label_class_t cache_column_set(Args&&... args) { + auto it = column_sets_.emplace(std::forward(args)...).first; + assert(std::is_sorted(it->begin(), it->end())); + return it - column_sets_.begin(); + } + + inline label_class_t cache_column(Column column) { + return cache_column_set(Columns(1, column)); + } + + // Fetch a label set given its index returned by `cache_column_set()` + inline const Columns& get_cached_column_set(label_class_t i) const { + assert(i < column_sets_.size()); + return column_sets_.data()[i]; + } + + std::string generate_column_set_str(label_class_t i, size_t spelling_size) const; + + CoordinateSet get_label_path_coords(const std::vector &path, + const std::vector &label_path) const; + + bool node_has_labels(node_index node, label_class_t target) const; + + private: + const DeBruijnGraph &graph_; + const Annotator &annotator_; + const annot::matrix::MultiIntMatrix *multi_int_; + const CanonicalDBG *canonical_; + + // keep a unique set of annotation rows + // the first element is the empty label set + VectorSet column_sets_; + // map node to index in |column_sets_| + VectorMap node_to_cols_; + // coordinate sets for all nodes in |node_to_cols_| in the same order + std::vector label_coords_; + // buffer of paths to later querying with fetch_queued_annotations() + std::vector> queued_paths_; +}; + +} // namespace align +} // namespace mtg::graph diff --git a/metagraph/src/graph/annotation_key_mapper.cpp b/metagraph/src/graph/annotation_key_mapper.cpp index 5b12caac20..6693fbd1ee 100644 --- a/metagraph/src/graph/annotation_key_mapper.cpp +++ b/metagraph/src/graph/annotation_key_mapper.cpp @@ -5,12 +5,14 @@ #include #include "representation/canonical_dbg.hpp" +#include "representation/hash/dbg_sshash.hpp" namespace mtg::graph { AnnotationKeyMapper::AnnotationKeyMapper(const DeBruijnGraph &graph) : graph_(graph), mode_(graph.get_mode()), - canonical_(dynamic_cast(&graph)) { + canonical_(dynamic_cast(&graph)), + sshash_(dynamic_cast(&graph)) { if (mode_ == DeBruijnGraph::PRIMARY) throw std::invalid_argument("Wrap PRIMARY graphs in CanonicalDBG before mapping annotations"); } @@ -18,6 +20,8 @@ AnnotationKeyMapper::AnnotationKeyMapper(const DeBruijnGraph &graph) auto AnnotationKeyMapper::canonical_key(node_index node, std::string_view kmer) const -> node_index { if (canonical_) return canonical_->get_base_node(node); + if (sshash_) + return sshash_->get_base_node(node); assert(kmer.size() == graph_.get_k()); node_index result = DeBruijnGraph::npos; graph_.map_to_nodes(kmer, [&](node_index key) { result = key; }); diff --git a/metagraph/src/graph/annotation_key_mapper.hpp b/metagraph/src/graph/annotation_key_mapper.hpp index 84a98a327f..c43f9b2fc6 100644 --- a/metagraph/src/graph/annotation_key_mapper.hpp +++ b/metagraph/src/graph/annotation_key_mapper.hpp @@ -9,6 +9,7 @@ namespace mtg::graph { class CanonicalDBG; +class DBGSSHash; // Map oriented graph nodes to the node IDs whose rows carry their annotation. // IDs remain one-based here; AnnotatedDBG::graph_to_anno_index converts to rows. @@ -18,7 +19,7 @@ class AnnotationKeyMapper { explicit AnnotationKeyMapper(const DeBruijnGraph &graph); bool needs_sequence() const { - return mode_ == DeBruijnGraph::CANONICAL && !canonical_; + return mode_ == DeBruijnGraph::CANONICAL && !canonical_ && !sshash_; } // For a nonzero node, supply its k-mer spelling when needs_sequence() is true. node_index key(node_index node, std::string_view kmer = {}) const { @@ -36,6 +37,7 @@ class AnnotationKeyMapper { const DeBruijnGraph &graph_; DeBruijnGraph::Mode mode_; const CanonicalDBG *canonical_; + const DBGSSHash *sshash_; }; } // namespace mtg::graph diff --git a/metagraph/src/graph/representation/base/sequence_graph.cpp b/metagraph/src/graph/representation/base/sequence_graph.cpp index 3ff04fc7a3..0446dd14cc 100644 --- a/metagraph/src/graph/representation/base/sequence_graph.cpp +++ b/metagraph/src/graph/representation/base/sequence_graph.cpp @@ -118,17 +118,39 @@ bool DeBruijnGraph::operator==(const DeBruijnGraph &other) const { } void DeBruijnGraph::traverse(node_index start, - const char *begin, - const char *end, + std::string_view seq, const std::function &callback, const std::function &terminate) const { assert(start >= 1 && start <= max_index()); - assert(end >= begin); if (terminate()) return; - for (; begin != end && !terminate(); ++begin) { - start = traverse(start, *begin); + for (char c : seq) { + if (terminate()) + return; + + start = traverse(start, c); + + if (!in_graph(start)) + return; + + callback(start); + } +} + +void DeBruijnGraph::traverse_back(node_index start, + std::string_view prefix_seq, + const std::function &callback, + const std::function &terminate) const { + assert(start >= 1 && start <= max_index()); + if (terminate()) + return; + + for (auto begin = prefix_seq.rbegin(); begin != prefix_seq.rend(); ++begin) { + if (terminate()) + return; + + start = traverse_back(start, *begin); if (!in_graph(start)) return; diff --git a/metagraph/src/graph/representation/base/sequence_graph.hpp b/metagraph/src/graph/representation/base/sequence_graph.hpp index c0bf4bfb57..7a87b6f9d2 100644 --- a/metagraph/src/graph/representation/base/sequence_graph.hpp +++ b/metagraph/src/graph/representation/base/sequence_graph.hpp @@ -175,11 +175,19 @@ class DeBruijnGraph : public SequenceGraph { // when the sequence is exhausted. // In canonical mode, non-canonical k-mers are NOT mapped to canonical ones. virtual void traverse(node_index start, - const char *begin, - const char *end, + std::string_view seq, const std::function &callback, const std::function &terminate = [](){ return false; }) const; + // Given a starting node and a sequence of edge labels, traverse the graph + // backward starting from the end of prefix_seq. The traversal is terminated + // once terminate() returns true or when the sequence is exhausted. + // In canonical mode, non-canonical k-mers are NOT mapped to canonical ones. + virtual void traverse_back(node_index start, + std::string_view prefix_seq, + const std::function &callback, + const std::function &terminate = [](){ return false; }) const; + typedef std::function&)> CallPath; /** diff --git a/metagraph/src/graph/representation/canonical_dbg.hpp b/metagraph/src/graph/representation/canonical_dbg.hpp index c0841d6eec..e0e49e9b7a 100644 --- a/metagraph/src/graph/representation/canonical_dbg.hpp +++ b/metagraph/src/graph/representation/canonical_dbg.hpp @@ -40,6 +40,13 @@ class CanonicalDBG : public DBGWrapper { return node <= offset_ ? node : node - offset_; } + virtual bool in_graph(node_index node) const override final { + if (node == npos || node > max_index()) + return false; + + return graph_->in_graph(get_base_node(node)); + } + /** * Methods from DeBruijnGraph */ diff --git a/metagraph/src/graph/representation/hash/dbg_sshash.cpp b/metagraph/src/graph/representation/hash/dbg_sshash.cpp index ddbfe8a145..69dddb5df6 100644 --- a/metagraph/src/graph/representation/hash/dbg_sshash.cpp +++ b/metagraph/src/graph/representation/hash/dbg_sshash.cpp @@ -37,13 +37,6 @@ static constexpr uint16_t bits_per_char = sshash::dna_uint_kmer_t::bit const std::string DBGSSHash::alphabet_ = kmer::KmerExtractor2Bit().alphabet; -constexpr DeBruijnGraph::node_index sshash_to_graph_index(uint64_t idx) { - return idx + 1; -} -constexpr uint64_t graph_index_to_sshash(DeBruijnGraph::node_index idx) { - return idx - 1; -} - size_t DBGSSHash::dict_size() const { return std::visit([](const auto& d) { return d.num_kmers(); }, dict_); } @@ -80,7 +73,7 @@ generate_succ_pred(const DBGSSHash& graph, uint64_t num_nodes) { [&](const auto& dict) { using kmer_t = typename std::decay_t::kmer_type; - ProgressBar progress_bar(graph.num_nodes(), "Checking succ/pred", + ProgressBar progress_bar(num_nodes, "Checking succ/pred", std::cerr, !common::get_verbose()); static const uint64_t BS = 1'048'576; @@ -107,7 +100,7 @@ generate_succ_pred(const DBGSSHash& graph, uint64_t num_nodes) { // a k-mer reached from its predecessor was already found at // this node by the neighbour lookup, so check only the others if (dict.lookup(kmer, with_rc).kmer_id - != graph_index_to_sshash(node)) + != DBGSSHash::graph_index_to_sshash(node)) has_duplicates.store(true, std::memory_order_relaxed); } @@ -123,7 +116,7 @@ generate_succ_pred(const DBGSSHash& graph, uint64_t num_nodes) { break; } else { found_fw = true; - next = sshash_to_graph_index(nb.forward[i].kmer_id); + next = DBGSSHash::sshash_to_graph_index(nb.forward[i].kmer_id); found_i = i; if (nb.forward[i].kmer_orientation < 0) next = graph.reverse_complement(next); @@ -150,7 +143,7 @@ generate_succ_pred(const DBGSSHash& graph, uint64_t num_nodes) { break; } else { found_bw = true; - prev = sshash_to_graph_index(nb.backward[i].kmer_id); + prev = DBGSSHash::sshash_to_graph_index(nb.backward[i].kmer_id); if (nb.backward[i].kmer_orientation < 0) prev = graph.reverse_complement(prev); } @@ -175,8 +168,12 @@ generate_succ_pred(const DBGSSHash& graph, uint64_t num_nodes) { has_duplicates.load()); } -DBGSSHash::DBGSSHash(const std::string& input_filename, size_t k, Mode mode, size_t num_chars) - : DBGSSHash(k, mode) { +DBGSSHash::DBGSSHash(const std::string& input_filename, + size_t k, + Mode mode, + size_t num_chars, + const std::filesystem::path &swap_dir) + : DBGSSHash(k, mode) { if (k <= 1) throw std::domain_error("k must be at least 2"); @@ -206,6 +203,8 @@ DBGSSHash::DBGSSHash(const std::string& input_filename, size_t k, Mode mode, siz build_config.verbose = common::get_verbose(); build_config.num_threads = get_num_threads(); + if (!swap_dir.empty()) + build_config.tmp_dirname = swap_dir; build_config.allow_reverse_complement_duplicates = mode == BASIC; // put SSHash's temp files under --disk-swap (default: ./) and clean them up @@ -240,6 +239,19 @@ void DBGSSHash::add_sequence(std::string_view sequence, throw std::logic_error("adding sequences not supported"); } +sdsl::bit_vector DBGSSHash::get_invalid_char_mask(std::string_view sequence) const { + sdsl::bit_vector invalid_char(sequence.size()); + + std::visit([&](const auto& dict) { + using kmer_t = typename std::decay_t::kmer_type; + for (size_t i = 0; i < sequence.size(); ++i) { + invalid_char[i] = !kmer_t::is_valid(sequence[i]); + } + }, dict_); + + return invalid_char; +} + template void map_to_nodes_with_rc_impl(const DBGSSHash& graph, const Dict& dict, @@ -272,11 +284,7 @@ void map_to_nodes_with_rc_impl(const DBGSSHash& graph, callback(ret_val); } } else { - std::vector invalid_char(n); - for (size_t i = 0; i < n; ++i) { - invalid_char[i] = !kmer_t::is_valid(sequence[i]); - } - + auto invalid_char = graph.get_invalid_char_mask(sequence); auto invalid_kmer = utils::drag_and_mark_segments(invalid_char, true, k); bool begin = true; @@ -362,6 +370,15 @@ DBGSSHash::node_index DBGSSHash::reverse_complement(node_index node) const { return str_kmer != rc_str_kmer ? node + dict_size() : node; } +DBGSSHash::node_index DBGSSHash::get_base_node(node_index node) const { + if (node == npos) + return npos; + + assert(in_graph(node)); + + return node > dict_size() ? node - dict_size() : node; +} + uint64_t DBGSSHash::num_nodes() const { return mode_ != CANONICAL ? dict_size() : dict_size() * 2; } @@ -747,6 +764,34 @@ std::string DBGSSHash::get_node_sequence(node_index node) const { return str_kmer; } +uint64_t DBGSSHash::get_sequence_id(node_index node) const { + assert(node <= max_index()); + if (node == npos) + return npos; + + if (node > dict_size()) { + assert(mode_ == CANONICAL); + size_t num_base_contigs = succ_is_next_.size() + pred_is_prev_.size() - succ_is_next_.num_set_bits() - pred_is_prev_.num_set_bits(); + return get_sequence_id(reverse_complement(node)) + num_base_contigs; + } + + // treat each unitig as a sequence + size_t srank = succ_is_next_.rank0(node); + size_t prank = pred_is_prev_.rank0(node); + + uint64_t contig_id = srank + prank; + + bool s_unitig_back = !succ_is_next_[node]; + bool s_unitig_front = !pred_is_prev_[node]; + + if (s_unitig_back && !s_unitig_front) { + assert(contig_id); + --contig_id; + } + + return contig_id + 1; +} + void DBGSSHash::serialize(std::ostream& out) const { essentials::generic_saver saver(out); @@ -764,6 +809,12 @@ void DBGSSHash::serialize(std::ostream& out) const { void DBGSSHash::serialize(const std::string& filename) const { std::string suffixed_filename = utils::make_suffix(filename, kExtension); std::ofstream fout(suffixed_filename, std::ios::binary); + + if (!fout.good()) { + common::logger->error("Failed to create file {}", suffixed_filename); + throw std::runtime_error("Invalid out file"); + } + serialize(fout); } diff --git a/metagraph/src/graph/representation/hash/dbg_sshash.hpp b/metagraph/src/graph/representation/hash/dbg_sshash.hpp index c07b760c0b..49e29823d8 100644 --- a/metagraph/src/graph/representation/hash/dbg_sshash.hpp +++ b/metagraph/src/graph/representation/hash/dbg_sshash.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -36,7 +37,8 @@ class DBGSSHash : public DeBruijnGraph { DBGSSHash(const std::string& input_filename, size_t k, Mode mode = BASIC, - size_t num_chars = 0); + size_t num_chars = 0, + const std::filesystem::path &swap_dir = {}); // SequenceGraph overrides void add_sequence( @@ -111,12 +113,18 @@ class DBGSSHash : public DeBruijnGraph { node_index node, const std::function& callback) const; + // Return the index of the unitig containing the node, or npos for npos. + // In canonical mode, unitigs of reverse-complement nodes follow those of + // the forward nodes. + uint64_t get_sequence_id(node_index node) const; + const std::string& alphabet() const override final { return alphabet_; } const dict_t& data() const { return dict_; } node_index reverse_complement(node_index node) const; + node_index get_base_node(node_index node) const; void adjacent_outgoing_nodes(node_index node, const std::function& callback) const override; @@ -134,6 +142,15 @@ class DBGSSHash : public DeBruijnGraph { adjacent_incoming_nodes_with_rc(node_index node, const std::function& callback) const; + sdsl::bit_vector get_invalid_char_mask(std::string_view sequence) const; + + static constexpr node_index sshash_to_graph_index(uint64_t idx) { + return idx + 1; + } + + static constexpr uint64_t graph_index_to_sshash(node_index idx) { + return idx - 1; + } private: static const std::string alphabet_; diff --git a/metagraph/src/graph/representation/succinct/boss.hpp b/metagraph/src/graph/representation/succinct/boss.hpp index 7782a92d99..806ce4db7d 100644 --- a/metagraph/src/graph/representation/succinct/boss.hpp +++ b/metagraph/src/graph/representation/succinct/boss.hpp @@ -481,6 +481,8 @@ class BOSS { inline std::tuple index_range(RandomAccessIt begin, RandomAccessIt end) const; + inline void call_tighter_ranges(edge_index rl, edge_index ru, + const std::function &callback) const; inline bool tighten_range(edge_index *rl, edge_index *ru, TAlphabet s) const; /** @@ -695,6 +697,30 @@ inline bool BOSS::tighten_range(edge_index *rl, edge_index *ru, TAlphabet s) con return true; } +inline void BOSS +::call_tighter_ranges(edge_index rl, edge_index ru, + const std::function &callback) const { + if (ru - rl + 1 >= alphabet.size() - 1) { + for (TAlphabet s = 1; s < alph_size; ++s) { + edge_index next_rl = rl; + edge_index next_ru = ru; + if (tighten_range(&next_rl, &next_ru, s)) + callback(next_rl, next_ru, s); + } + } else { + for (edge_index e = rl; e <= ru; ++e) { + if (TAlphabet s = get_W(e) % alph_size) { + edge_index next_rl = rl; + edge_index next_ru = ru; + bool ret_val = tighten_range(&next_rl, &next_ru, s); + std::ignore = ret_val; + assert(ret_val); + callback(next_rl, next_ru, s); + } + } + } +} + template BOSS::edge_index BOSS::index(RandomAccessIt begin, RandomAccessIt end) const { static_assert(std::is_same_v, TAlphabet>, diff --git a/metagraph/src/graph/representation/succinct/dbg_succinct.cpp b/metagraph/src/graph/representation/succinct/dbg_succinct.cpp index d3ddf10ea9..1d5ff3bfce 100644 --- a/metagraph/src/graph/representation/succinct/dbg_succinct.cpp +++ b/metagraph/src/graph/representation/succinct/dbg_succinct.cpp @@ -393,12 +393,10 @@ ::call_nodes_with_suffix_matching_longest_prefix( } void DBGSuccinct::traverse(node_index start, - const char *begin, - const char *end, + std::string_view seq, const std::function &callback, const std::function &terminate) const { assert(in_graph(start)); - assert(end >= begin); if (terminate()) return; @@ -407,7 +405,7 @@ void DBGSuccinct::traverse(node_index start, assert(edge); BOSS::TAlphabet w; - for (; begin != end && !terminate() && (w = boss_graph_->get_W(edge)); ++begin) { + for (auto begin = seq.begin(); begin != seq.end() && !terminate() && (w = boss_graph_->get_W(edge)); ++begin) { // stop traversal if the character is invalid if (boss_graph_->encode(*begin) == boss_graph_->alph_size) return; diff --git a/metagraph/src/graph/representation/succinct/dbg_succinct.hpp b/metagraph/src/graph/representation/succinct/dbg_succinct.hpp index 2f243189f4..6f0f1e4edf 100644 --- a/metagraph/src/graph/representation/succinct/dbg_succinct.hpp +++ b/metagraph/src/graph/representation/succinct/dbg_succinct.hpp @@ -91,12 +91,11 @@ class DBGSuccinct : public DeBruijnGraph { size_t max_num_allowed_matches = std::numeric_limits::max()) const; // Given a starting node, traverse the graph forward following the edge - // sequence delimited by begin and end. Terminate the traversal if terminate() + // sequence seq. Terminate the traversal if terminate() // returns true, or if the sequence is exhausted. // In canonical mode, non-canonical k-mers are NOT mapped to canonical ones virtual void traverse(node_index start, - const char *begin, - const char *end, + std::string_view seq, const std::function &callback, const std::function &terminate = [](){ return false; }) const override final; diff --git a/metagraph/src/pybindings/pymetagraph_bindings.cpp b/metagraph/src/pybindings/pymetagraph_bindings.cpp new file mode 100644 index 0000000000..11a891710b --- /dev/null +++ b/metagraph/src/pybindings/pymetagraph_bindings.cpp @@ -0,0 +1,361 @@ +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "cli/align.hpp" +#include "cli/config/config.hpp" +#include "cli/load/load_annotated_graph.hpp" +#include "cli/load/load_graph.hpp" +#include "common/logger.hpp" +#include "common/threads/threading.hpp" +#include "graph/alignment_redone/aln_match.hpp" +#include "graph/alignment_redone/aln_query.hpp" +#include "graph/alignment_redone/aln_seeder.hpp" +#include "graph/alignment_redone/annotation_buffer.hpp" +#include "graph/annotated_dbg.hpp" +#include "graph/representation/canonical_dbg.hpp" + +namespace py = pybind11; + +namespace { + +// Converts Python kwargs into a fake argv so metagraph's existing CLI argument +// parser (mtg::cli::Config) can be reused unchanged instead of reimplementing +// config parsing/validation here. argv[0] is a dummy program name; argv[1] (the +// CLI "identity") is picked based on "method", since Config's validation rules +// differ: +// - "query" -> "server_query": the one identity whose validation allows +// loading a graph/annotation from -i/-a with no positional FASTA files +// (plain "query" always demands input files), but it also *requires* +// exactly one annotator -- fine here, since label-presence queries can't +// work without one anyway. +// - "align" -> "align": also allows -i (+ no positional files) with no +// validation-time exit(), but unlike "server_query" it does not require an +// annotator at all, matching how `metagraph align` itself only needs a +// graph. This is what lets Index(method="align") skip annotation entirely. +// Neither identity's *parsing* (as opposed to validation) is gated on which +// one is picked -- e.g. -a/--num-top-labels/etc. are recognized the same way +// under both -- so this only changes which validation rules apply. +// +// "method", "connect_anchors" and "extend_chains" are consumed here rather +// than forwarded to Config: "method" selects which pymetagraph binding code +// path runs (label-presence vs. real alignment), and the other two only +// configure the alignment path, which no CLI command uses. +struct BindingOptions { + std::string method = "query"; + bool connect_anchors = true; + bool extend_chains = true; +}; + +std::vector kwargs_to_argv(const py::kwargs &kwargs, BindingOptions *options) { + static const std::unordered_map kFlagAliases = { + { "input", "-i" }, + { "annotator", "-a" }, + { "threads", "-p" }, + { "num_top_labels", "--num-top-labels" }, + { "discovery_fraction", "--min-kmers-fraction-label" }, + { "presence_fraction", "--min-kmers-fraction-graph" }, + { "seed_length", "--align-min-seed-length" }, + { "max_num_nodes_per_seq_char", "--align-max-nodes-per-seq-char" }, + { "min_exact_match", "--align-min-exact-match" }, + { "max_num_seeds_per_locus", "--align-max-num-seeds-per-locus" }, + { "xdrop", "--align-xdrop" }, + { "rel_score_cutoff", "--align-rel-score-cutoff" }, + }; + + // args[1] (the identity) is filled in below once "method" has been seen, + // since kwargs iteration order isn't guaranteed to visit it first. + std::vector args = { "pymetagraph", "" }; + + for (auto item : kwargs) { + auto key = item.first.cast(); + + if (key == "method") { + options->method = py::str(item.second).cast(); + continue; + } + if (key == "connect_anchors") { + options->connect_anchors = item.second.cast(); + continue; + } + if (key == "extend_chains") { + options->extend_chains = item.second.cast(); + continue; + } + + auto it = kFlagAliases.find(key); + if (it == kFlagAliases.end()) { + throw std::invalid_argument("pymetagraph: unknown keyword argument '" + key + "'"); + } + args.push_back(it->second); + args.push_back(py::str(item.second).cast()); + } + + if (options->method != "query" && options->method != "align") { + throw std::invalid_argument( + "pymetagraph: unknown method '" + options->method + "' (expected 'query' or 'align')"); + } + + args[1] = (options->method == "align") ? "align" : "server_query"; + + // Unlike "server_query", the "align" identity's validation unconditionally + // requires at least one positional input file (config.cpp: "No input + // file(s) passed"), since it's normally used as `align -i graph.dbg + // reads.fq`. We never read config.fnames ourselves -- queries come in via + // Index::query/query_batch instead -- so a placeholder that's never opened + // is enough to satisfy that check. + if (options->method == "align") + args.push_back("__pymetagraph_placeholder__"); + + return args; +} + +} // namespace + +// An alignment-like result. Only the fields Readfish's decision logic reads +// (`ctg`, `r_st`, `r_en`, `strand`) are populated; r_st/r_en/strand are dummy +// placeholders since pymetagraph reports label presence, not real coordinates. +struct Alignment { + std::string ctg; + int64_t r_st = 0; + int64_t r_en = 0; + int strand = 1; + + Alignment() = default; + Alignment(std::string ctg, int64_t r_st, int64_t r_en, int strand) + : ctg(std::move(ctg)), r_st(r_st), r_en(r_en), strand(strand) {} +}; + +// Loads a prebuilt metagraph graph, plus an optional single annotation file, +// and answers queries against it, using one of two unrelated underlying +// methods picked at construction time via the "method" kwarg: +// - "query" (default): AnnotatedDBG::get_top_labels, the same label-presence +// primitive backing the CLI's `query` command. Requires an annotation. +// Read-only, no mutable per-call state, safe to call concurrently from +// multiple threads. +// - "align": the alignment_redone seed-and-extend pipeline. Neither the CLI's +// `align` command nor the server's `/align` endpoint uses this engine. +// Reports real r_st/r_en/strand instead of dummy placeholders. Each call +// builds its own AnnotationBuffer, since that type caches mutable +// per-query state and is not meant to be shared across concurrent queries. +// Annotation is optional here, mirroring align.cpp's own behavior when no +// "-a" is given: without one, plain ExactSeeder/Extender is used (no label +// lookups at all, so the O(num_labels) annotation-fetch cost is avoided +// entirely), and the reported "ctg" degrades to a generic "mapped"/"*" +// hit-or-miss marker instead of an actual label. +class Index { + public: + explicit Index(const py::kwargs &kwargs) { + BindingOptions options; + auto arg_strings = kwargs_to_argv(kwargs, &options); + method_ = (options.method == "align") ? Method::kAlign : Method::kQuery; + + std::vector argv; + argv.reserve(arg_strings.size()); + for (auto &s : arg_strings) + argv.push_back(s.data()); + + mtg::cli::Config config(static_cast(argv.size()), argv.data()); + + if (config.infbase.empty()) + throw std::invalid_argument("pymetagraph: missing required argument 'input'"); + if (config.infbase_annotators.size() > 1) { + throw std::invalid_argument( + "pymetagraph: at most one 'annotator' file is supported"); + } + if (method_ == Method::kQuery && config.infbase_annotators.empty()) { + throw std::invalid_argument( + "pymetagraph: method 'query' requires an 'annotator' file " + "('align' can be used without one)"); + } + + graph_ = mtg::cli::load_critical_dbg(config.infbase); + if (!config.infbase_annotators.empty()) + anno_dbg_ = mtg::cli::initialize_annotated_dbg(graph_, config); + + num_top_labels_ = config.num_top_labels; + discovery_fraction_ = config.discovery_fraction; + presence_fraction_ = config.presence_fraction; + + if (method_ == Method::kAlign) { + connect_anchors_ = options.connect_anchors; + extend_chains_ = options.extend_chains; + + if (anno_dbg_) { + aligner_config_ = mtg::cli::initialize_redone_aligner_config(config, anno_dbg_->get_graph()); + } else { + // Mirrors align.cpp's wrap_graph: a PRIMARY-mode graph must be + // wrapped into CanonicalDBG before it can be aligned to directly. + // When an annotation is loaded, initialize_annotated_dbg does this + // wrapping internally; here we have to do it ourselves. + align_graph_ = graph_; + if (align_graph_->get_mode() == mtg::graph::DeBruijnGraph::PRIMARY) + align_graph_ = std::make_shared(align_graph_); + + aligner_config_ = mtg::cli::initialize_redone_aligner_config(config, *align_graph_); + } + } + + thread_pool_ = std::make_unique(std::max(1u, get_num_threads())); + } + + Alignment query(const std::string &seq) const { + if (seq.empty()) + return Alignment("*", 0, 0, 1); + + return method_ == Method::kAlign ? query_by_alignment(seq) : query_by_label(seq); + } + + std::vector query_batch(const std::vector &seqs) const { + std::vector results(seqs.size()); + std::vector> futures; + futures.reserve(seqs.size()); + + for (size_t i = 0; i < seqs.size(); ++i) { + // A read that throws is reported as unmapped, so it neither + // terminates the process nor fails the rest of the batch. + futures.push_back(thread_pool_->enqueue([&, i]() { + try { + results[i] = query(seqs[i]); + } catch (const std::exception &e) { + mtg::common::logger->warn("pymetagraph: read {} failed: {}", i, e.what()); + results[i] = Alignment("*", 0, 0, 1); + } + })); + } + for (auto &f : futures) + f.get(); + + return results; + } + + private: + enum class Method { kQuery, kAlign }; + + Alignment query_by_label(const std::string &seq) const { + auto top_labels = anno_dbg_->get_top_labels(seq, num_top_labels_, + discovery_fraction_, + presence_fraction_); + if (top_labels.empty()) + return Alignment("*", 0, 0, 1); + + return Alignment(top_labels.front().first, 0, static_cast(seq.size()), 1); + } + + Alignment query_by_alignment(const std::string &seq) const { + namespace align_redone = mtg::graph::align_redone; + + const auto &graph = anno_dbg_ ? anno_dbg_->get_graph() : *align_graph_; + align_redone::Query aln_query(graph, seq); + + std::optional anno_buffer; + std::unique_ptr seeder; + if (anno_dbg_) { + anno_buffer.emplace(graph, anno_dbg_->get_annotator()); + seeder = std::make_unique(*anno_buffer, aln_query, *aligner_config_); + } else { + seeder = std::make_unique(aln_query, *aligner_config_); + } + + std::vector paths; + auto aln_callback = [&](align_redone::Alignment &&aln) { + paths.emplace_back(std::move(aln)); + }; + + if (extend_chains_) { + std::unique_ptr extender; + if (anno_dbg_) { + extender = std::make_unique(*anno_buffer, aln_query, *aligner_config_); + } else { + extender = std::make_unique(aln_query, *aligner_config_); + } + align_redone::align_query(aln_query, *seeder, *extender, aln_callback, connect_anchors_); + } else { + align_redone::align_query(aln_query, *seeder, aln_callback, connect_anchors_); + } + + if (paths.empty()) + return Alignment("*", 0, 0, 1); + + const auto &best = *std::max_element(paths.begin(), paths.end(), + [&](const align_redone::Alignment &a, const align_redone::Alignment &b) { + return align_redone::score_match(a, *aligner_config_) + < align_redone::score_match(b, *aligner_config_); + }); + + std::string label; + if (anno_dbg_) { + label = "*"; + if (!best.get_label_classes().empty()) { + auto label_class = best.get_label_classes().front(); + if (label_class != align_redone::Anchor::nannot) { + const auto &columns = anno_buffer->get_cached_column_set(label_class); + if (!columns.empty()) + label = anno_buffer->get_annotator().get_label_encoder().decode(columns.front()); + } + } + } else { + // Without an annotation there are no labels to decode; a non-"*" + // ctg here just signals "this read aligned to the graph". + label = "mapped"; + } + + // Query coordinates are always reported on the original (forward) read, + // mirroring the reverse-complement clip positions back when needed. + auto qlen = static_cast(seq.size()); + int64_t r_st, r_en; + if (!best.get_orientation()) { + r_st = static_cast(best.get_clipping()); + r_en = qlen - static_cast(best.get_end_clipping()); + } else { + r_st = static_cast(best.get_end_clipping()); + r_en = qlen - static_cast(best.get_clipping()); + } + int strand = best.get_orientation() ? -1 : 1; + + return Alignment(label, r_st, r_en, strand); + } + + Method method_ = Method::kQuery; + std::shared_ptr graph_; + std::unique_ptr anno_dbg_; + // Only populated when method_ == Method::kAlign and there's no annotation; + // may be a CanonicalDBG wrapper around graph_ (see constructor). + std::shared_ptr align_graph_; + std::unique_ptr thread_pool_; + size_t num_top_labels_ = 1; + double discovery_fraction_ = 0.7; + double presence_fraction_ = 0.0; + + // Only populated when method_ == Method::kAlign. + std::optional aligner_config_; + bool connect_anchors_ = true; + bool extend_chains_ = true; +}; + +PYBIND11_MODULE(_pymetagraph_core, m) { + py::class_(m, "Alignment") + .def(py::init<>()) + .def(py::init(), + py::arg("ctg"), py::arg("r_st"), py::arg("r_en"), py::arg("strand")) + .def_readonly("ctg", &Alignment::ctg) + .def_readonly("r_st", &Alignment::r_st) + .def_readonly("r_en", &Alignment::r_en) + .def_readonly("strand", &Alignment::strand); + + py::class_(m, "Index") + .def(py::init()) + .def("query", &Index::query, py::arg("seq"), + py::call_guard()) + .def("query_batch", &Index::query_batch, py::arg("seqs"), + py::call_guard()); +} diff --git a/metagraph/tests/annotation/test_aligner_labeled.cpp b/metagraph/tests/annotation/test_aligner_labeled.cpp index 2e9653ec9f..fee23d8716 100644 --- a/metagraph/tests/annotation/test_aligner_labeled.cpp +++ b/metagraph/tests/annotation/test_aligner_labeled.cpp @@ -806,4 +806,54 @@ TYPED_TEST(LabeledAlignerTest, CanonicalBackwardExtensionMismatch) { } #endif +#if ! _PROTEIN_GRAPH +// Seeds shorter than k on a succinct graph (k = 5, min seed length 3), with +// A = TTACGAAA and B = ACGTTGCA. ACGT, TGCA and GTTG occur only in B, TACGA +// only in A. ACGT is a suffix only of the dummy source k-mer $ACGT that starts B. +// The existing aligner matches sub-k seeds against node suffixes, so on the +// basic graph it finds no alignment for ACGT. Disabled until that is fixed. +TEST(LabeledAlignerTest, DISABLED_SubkSeedLabels) { + size_t k = 5; + const std::vector sequences { "TTACGAAA", "ACGTTGCA" }; + const std::vector labels { "A", "B" }; + const std::vector> queries { + { "ACGT", "B" }, { "TACGA", "A" }, { "CGTTG", "B" }, { "TGCA", "B" } + }; + + for (auto mode : { DeBruijnGraph::BASIC, DeBruijnGraph::CANONICAL, DeBruijnGraph::PRIMARY }) { + auto anno_graph = build_anno_graph>( + k, sequences, labels, mode); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + config.min_seed_length = 3; + LabeledAligner<> aligner(anno_graph->get_graph(), config, anno_graph->get_annotator()); + const auto &label_encoder = anno_graph->get_annotator().get_label_encoder(); + + for (const auto &[query, label] : queries) { + bool found = false; + for (const auto &alignment : aligner.align(query)) { + // every label of an alignment, including a partial one, must + // contain its sequence on either strand + std::string sequence(alignment.get_sequence()); + std::string sequence_rc = sequence; + reverse_complement(sequence_rc.begin(), sequence_rc.end()); + std::vector alignment_labels; + for (auto column : alignment.label_columns) { + alignment_labels.emplace_back(label_encoder.decode(column)); + const std::string &labeled = sequences[alignment_labels.back() == "A" ? 0 : 1]; + EXPECT_TRUE(labeled.find(sequence) != std::string::npos + || labeled.find(sequence_rc) != std::string::npos) + << mode << " " << query << " " << alignment; + } + found |= (sequence == query || sequence_rc == query) + && alignment_labels == std::vector{ label }; + } + EXPECT_TRUE(found) << mode << " " << query; + } + } +} + +#endif + } // namespace diff --git a/metagraph/tests/annotation/test_aligner_labeled_redone.cpp b/metagraph/tests/annotation/test_aligner_labeled_redone.cpp new file mode 100644 index 0000000000..35bc432059 --- /dev/null +++ b/metagraph/tests/annotation/test_aligner_labeled_redone.cpp @@ -0,0 +1,647 @@ +#include + +#include "../graph/all/test_dbg_helpers.hpp" +#include "../test_helpers.hpp" +#include "test_annotated_dbg_helpers.hpp" + +#include + +#include "annotation/representation/column_compressed/annotate_column_compressed.hpp" +#include "annotation/annotation_converters.hpp" +#include "common/seq_tools/reverse_complement.hpp" +#include "graph/alignment_redone/aln_query.hpp" +#include "graph/alignment_redone/aln_match.hpp" +#include "graph/alignment_redone/aligner_config.hpp" +#include "graph/alignment_redone/aln_seeder.hpp" +#include "graph/representation/canonical_dbg.hpp" +#include "graph/annotated_dbg.hpp" +#include "kmer/alphabets.hpp" +#include "seq_io/sequence_io.hpp" + +namespace { + +using namespace mtg; +using namespace mtg::graph; +using namespace mtg::graph::align_redone; +using namespace mtg::test; +using namespace mtg::kmer; + + +void run_alignment(const AnnotatedDBG &anno_graph, + DBGAlignerConfig config, + std::string_view query, + const std::vector> &mappings, + size_t end_trim = 0, + bool needs_extension = false, + bool needs_extension_long_seed = false) { + const auto &graph = anno_graph.get_graph(); + const auto &label_encoder = anno_graph.get_annotator().get_label_encoder(); + AnnotationBuffer anno_buffer(graph, anno_graph.get_annotator()); + + size_t k = graph.get_k(); + if (config.min_seed_length == 0) + config.min_seed_length = k; + + for (auto mx : { k, std::numeric_limits::max() }) { + config.max_seed_length = std::max(mx, config.min_seed_length); + bool check_chaining = mx == std::numeric_limits::max() + ? !needs_extension_long_seed + : !needs_extension; + + Query aln_query(graph, query); + LabeledSeeder seeder(anno_buffer, aln_query, config); + std::vector paths; + LabeledExtender extender(anno_buffer, aln_query, config); + std::vector paths_no_extend = seeder.get_inexact_anchors(); + std::vector paths_no_align = seeder.get_inexact_anchors(false); + for (const auto &base_path : paths_no_extend) { + extender.extend(base_path, [&](Alignment&& path) { + paths.emplace_back(std::move(path)); + }); + } + + auto aln_sort = [&config,&anno_buffer](const auto &a, const auto &b) { + return std::make_tuple(score_match(a, config), + anno_buffer.get_cached_column_set(b.get_label_classes()[0])[0]) + > std::make_tuple(score_match(b, config), + anno_buffer.get_cached_column_set(a.get_label_classes()[0])[0]); + }; + + std::sort(paths_no_extend.begin(), paths_no_extend.end(), aln_sort); + std::sort(paths_no_align.begin(), paths_no_align.end(), aln_sort); + std::sort(paths.begin(), paths.end(), aln_sort); + + ASSERT_LE(mappings.size(), paths.size()) << mx; + paths.resize(mappings.size()); + + tsl::hopscotch_set seen_labels; + + for (size_t i = 0; i < mappings.size(); ++i) { + const auto &[label, reference, cigar_str] = mappings[i]; + Cigar cigar(cigar_str); + + auto check_ref = [&](const Alignment &path, const std::string &type) { + const auto &label = std::get<0>(mappings[i]); + std::string reference = std::get<1>(mappings[i]); + std::string cigar_str = std::get<2>(mappings[i]); + + if (path.get_orientation() && !path.get_end_trim()) { + ::reverse_complement(reference.begin(), reference.end()); + cigar = Cigar(cigar_str); + std::reverse(cigar.data().begin(), cigar.data().end()); + cigar_str = cigar.to_string(); + } + + if (label.size()) { + // make sure the labels are correct before moving on + size_t target_column = label_encoder.encode(label); + + size_t old_num_column_sets = anno_buffer.num_column_sets(); + size_t target = anno_buffer.cache_column(target_column); + ASSERT_EQ(old_num_column_sets, anno_buffer.num_column_sets()); + ASSERT_LE(target, old_num_column_sets); + + const auto &path_targets = path.get_label_classes(); + std::string path_target_labels = anno_buffer.generate_column_set_str(path_targets[0], path.get_spelling().size()); + + if (path_targets.size() > 1) { + ASSERT_TRUE(std::all_of(path_targets.begin(), path_targets.end(), + [&](auto node_target) { return node_target == path_targets[0]; })) + << path << "\t" << path_target_labels << "\t" << fmt::format("{}", fmt::join(path_targets, ",")); + } + // for (auto node_target : path_targets) { + // for (auto column : anno_buffer.get_cached_column_set(node_target)) { + // path_target_labels += label_encoder.decode(column) + ","; + // } + // path_target_labels += ";"; + // } + ASSERT_TRUE(std::all_of(path_targets.begin(), path_targets.end(), + [&](auto node_target) { return node_target == target; })) + << path << "\t" << path_target_labels << " vs. " << label << "\t" + << fmt::format("{}", fmt::join(path_targets, ",")) << " vs. " << target; + } + + ASSERT_EQ(cigar.to_string(), path.get_cigar().to_string()) << label << "\t" << mx << "\t" << type; + if (reference.size()) { + ASSERT_TRUE(cigar.is_valid(reference, path.get_query())); + EXPECT_EQ(reference.size() - k + 1 + end_trim, path.get_path().size()) << label << "\t" << mx << "\t" << type; + EXPECT_EQ(reference, path.get_spelling()) << label << "\t" << mx << "\t" << type; + } + + EXPECT_EQ(end_trim, path.get_end_trim()) << label << "\t" << mx << "\t" << type; + if (reference.size()) { + EXPECT_EQ(config.score_cigar(reference, path.get_query(), cigar), + score_match(path, config)) << label << "\t" << mx << "\t" << type; + } + + EXPECT_EQ(cigar.get_clipping(), path.get_clipping()) << label << "\t" << mx << "\t" << type; + EXPECT_EQ(cigar.get_end_clipping(), path.get_end_clipping()) << label << "\t" << mx << "\t" << type; + }; + + check_ref(paths[i], "extend"); + + if (check_chaining) { + // this alignment should work with chaining alone + ASSERT_LT(i, paths_no_extend.size()); + check_ref(paths_no_extend[i], "chain"); + + ASSERT_LT(i, paths_no_align.size()); + EXPECT_EQ(cigar.to_string(), paths_no_align[i].get_cigar().to_string()) << label << "\t" << mx << "\tnoalign"; + } + } + } +} + +template +class LabeledAlignerRedoneTest : public ::testing::Test {}; + +typedef ::testing::Types>, + std::pair>, + std::pair>, + std::pair, + std::pair, + std::pair> FewGraphAnnotationPairTypes; + +TYPED_TEST_SUITE(LabeledAlignerRedoneTest, FewGraphAnnotationPairTypes); + +TYPED_TEST(LabeledAlignerRedoneTest, SimpleLinearGraph) { + size_t k = 4; + /* + A A B B B B + GCAA-CAAT-AATG-ATGC-TGCT-GCTT + */ + const std::vector sequences { + "GCAAT", + "AATGCTT" + }; + const std::vector labels { "A", "B" }; + + auto anno_graph = build_anno_graph(k, sequences, labels); + + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("A"), + anno_graph->get_annotator().get_label_encoder().encode("B")); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + + std::vector>>> exp_alignments { + { std::string("GCAATGCTT"), { + std::make_tuple(std::string("B"), std::string("AATGCTT"), std::string("2S7=")), + std::make_tuple(std::string("A"), std::string("GCAAT"), std::string("5=4S")), + } } + }; + + for (const auto &[query, mappings] : exp_alignments) { + run_alignment(*anno_graph, config, query, mappings); + } +} + +TYPED_TEST(LabeledAlignerRedoneTest, SimpleTangleGraph) { + size_t k = 3; + /* B AB AB + CGA GCC-CCT + \ BC BC ABC ABC/ + GAA-AAT-ATG-TGC + C/ \ C C + GGA GCA-CAT + */ + const std::vector sequences { + "ATGCCT", + "CGAATGCCT", + "GGAATGCAT" + }; + const std::vector labels { "A", "B", "C" }; + + auto anno_graph = build_anno_graph(k, sequences, labels); + + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("A"), + anno_graph->get_annotator().get_label_encoder().encode("B")); + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("B"), + anno_graph->get_annotator().get_label_encoder().encode("C")); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -2, -2); + config.gap_opening_penalty = -4; + config.gap_extension_penalty = -4; + config.left_end_bonus = 6; + config.right_end_bonus = 6; + + std::vector>>> exp_alignments { + { std::string("CGAATGCAT"), { // RC: ATGCATTCG + std::make_tuple(std::string("B"), std::string("CGAATGCCT"), std::string("7=1X1=")), + std::make_tuple(std::string("C"), std::string("GGAATGCAT"), std::string("1X8=")), + } } + }; + + for (const auto &[query, mappings] : exp_alignments) { + run_alignment(*anno_graph, config, query, mappings, 0, true, true); + } +} + +TYPED_TEST(LabeledAlignerRedoneTest, SimpleTangleGraphSuffixSeed) { + if constexpr(!std::is_base_of_v + && !std::is_base_of_v) { + return; + } + + size_t k = 7; + /* + B B B B B B B B AB AB + TATTGCC-ATTGCCG-TTGCCGA-TGCCGAT-GCCGATT-CCGATTG GATTGCC-ATTGCCT-TTGCCTT-TGCCTTG + \ BC / + CGATTGC + C C C C C C / \ C C C C + TATTGGC-ATTGGCG-TTGGCGA-TGGCGAT-GGCGATT-GCGATTG GATTGCA-ATTGCAT-TTGCATT-TGCATTG + */ + const std::vector sequences { + "TTGCCTTG", + "TATTGCCGATTGCCTTG", + "TATTGGCGATTGCATTG" + }; + const std::vector labels { "A", "B", "C" }; + + auto anno_graph = build_anno_graph(k, sequences, labels); + + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("A"), + anno_graph->get_annotator().get_label_encoder().encode("B")); + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("B"), + anno_graph->get_annotator().get_label_encoder().encode("C")); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + config.min_seed_length = 5; + config.left_end_bonus = 5; + config.right_end_bonus = 5; + + { + std::vector>>> exp_alignments { + { std::string("TATTGACGATTGCATTG"), { + std::make_tuple(std::string("C"), std::string("TATTGGCGATTGCATTG"), std::string("5=1X11=")), + } } + }; + + for (size_t i = 0; i < exp_alignments.size(); ++i) { + bool first_seed_matches = true; + if constexpr(std::is_same_v) { + const auto &sshash = static_cast(anno_graph->get_graph()); + std::visit([&](const auto &dict) { + ASSERT_GE(config.min_seed_length, dict.m()); + const auto &[query, mappings] = exp_alignments[i]; + ASSERT_LT(0u, mappings.size()); + using kmer_t = typename std::decay_t::kmer_type; + auto first_kmer_ref = sshash::util::string_to_uint_kmer(std::get<1>(mappings[0]).data(), k); + auto first_kmer_query = sshash::util::string_to_uint_kmer(query.data(), k); + // we can only find the first seed of it's also the minimizer + if (sshash::util::compute_minimizer(first_kmer_ref, k, dict.m(), dict.hasher()).pos_in_seq + || sshash::util::compute_minimizer(first_kmer_query, k, dict.m(), dict.hasher()).pos_in_seq) { + first_seed_matches = false; + } + }, sshash.data()); + } + + const auto &[query, mappings] = exp_alignments[i]; + if (first_seed_matches) { + run_alignment(*anno_graph, config, query, mappings); + } else { + run_alignment(*anno_graph, config, query, mappings, 0, true, true); + } + } + } + { + std::vector>>> exp_alignments { + { std::string("TATTGCCGATTGCATTG"), { + std::make_tuple(std::string("B"), std::string("TATTGCCGATTGCCTTG"), std::string("13=1X3=")), + std::make_tuple(std::string("C"), std::string("TATTGGCGATTGCATTG"), std::string("5=1X11=")), + } } + }; + + for (const auto &[query, mappings] : exp_alignments) { + run_alignment(*anno_graph, config, query, mappings, 0, true, true); + } + } +} + +#if ! _PROTEIN_GRAPH +TYPED_TEST(LabeledAlignerRedoneTest, CanonicalTangleGraph) { + if constexpr(!std::is_base_of_v + && !std::is_base_of_v) { + return; + } + + size_t k = 7; + const std::vector sequences { + "CAAGGCAA", + "CAAGGCAATCGGCAATA", + "CAATGCAATCGCCAATA" + }; + const std::vector labels { "A", "B", "C" }; + + for (DeBruijnGraph::Mode mode : { DeBruijnGraph::CANONICAL, DeBruijnGraph::PRIMARY }) { + auto anno_graph = build_anno_graph( + k, sequences, labels, mode + ); + + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("A"), + anno_graph->get_annotator().get_label_encoder().encode("B")); + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("B"), + anno_graph->get_annotator().get_label_encoder().encode("C")); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + config.min_seed_length = 5; + config.left_end_bonus = 5; + config.right_end_bonus = 5; + + { + std::vector>>> exp_alignments { + { std::string("TATTGACGATTGCATTG"), { + std::make_tuple(std::string("C"), std::string("TATTGGCGATTGCATTG"), std::string("5=1X11=")), + } } + }; + + for (size_t i = 0; i < exp_alignments.size(); ++i) { + bool first_seed_matches = true; + if constexpr(std::is_same_v) { + const auto &sshash = static_cast(anno_graph->get_graph()); + std::visit([&](const auto &dict) { + ASSERT_GE(config.min_seed_length, dict.m()); + const auto &[query, mappings] = exp_alignments[i]; + ASSERT_LT(0u, mappings.size()); + using kmer_t = typename std::decay_t::kmer_type; + auto first_kmer_ref = sshash::util::string_to_uint_kmer(std::get<1>(mappings[0]).data(), k); + auto first_kmer_query = sshash::util::string_to_uint_kmer(query.data(), k); + // we can only find the first seed of it's also the minimizer + if (sshash::util::compute_minimizer(first_kmer_ref, k, dict.m(), dict.hasher()).pos_in_seq + || sshash::util::compute_minimizer(first_kmer_query, k, dict.m(), dict.hasher()).pos_in_seq) { + first_seed_matches = false; + } + }, sshash.data()); + } + + const auto &[query, mappings] = exp_alignments[i]; + if (first_seed_matches) { + run_alignment(*anno_graph, config, query, mappings); + } else { + run_alignment(*anno_graph, config, query, mappings, 0, true, true); + } + } + } + { + std::vector>>> exp_alignments { + { std::string("TATTGCCGATTGCATTG"), { + std::make_tuple(std::string("B"), std::string("TATTGCCGATTGCCTTG"), std::string("13=1X3=")), + std::make_tuple(std::string("C"), std::string("TATTGGCGATTGCATTG"), std::string("5=1X11=")), + } } + }; + + for (const auto &[query, mappings] : exp_alignments) { + run_alignment(*anno_graph, config, query, mappings, 0, true, true); + } + } + } +} +#endif + +TYPED_TEST(LabeledAlignerRedoneTest, SimpleTangleGraphCoords) { + // TODO: for now, not implemented for other annotators + if constexpr(!std::is_same_v> + && !std::is_same_v) { + return; + } + + size_t k = 3; + /* B AB AB + CGA GCC-CCT + \ BC BC BC ABC/ + GAA-AAT-ATG-TGC + C/ \ C C + GGA GCA-CAT + */ + const std::vector sequences { + "TGCCT", + "CGAATGCCT", + "GGAATGCAT" + }; + const std::vector labels { "A", "B", "C" }; + + auto anno_graph = build_anno_graph( + k, sequences, labels, DeBruijnGraph::BASIC, true + ); + + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("A"), + anno_graph->get_annotator().get_label_encoder().encode("B")); + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("B"), + anno_graph->get_annotator().get_label_encoder().encode("C")); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + config.left_end_bonus = 5; + config.right_end_bonus = 5; + + std::vector>>> exp_alignments { + { std::string("CGAATGCAT"), { // RC: ATGCATTCG + std::make_tuple(std::string("B"), std::string("CGAATGCCT"), std::string("7=1X1=")), + std::make_tuple(std::string("C"), std::string("GGAATGCAT"), std::string("1X8=")), + } } + }; + + for (const auto &[query, mappings] : exp_alignments) { + run_alignment(*anno_graph, config, query, mappings, 0, true, true); + } +} + +TYPED_TEST(LabeledAlignerRedoneTest, SimpleTangleGraphCoordsMiddle) { + // TODO: for now, not implemented for other annotators + if constexpr(!std::is_same_v> + && !std::is_same_v) { + return; + } + + size_t k = 3; + /* B AB AB + CGA GCC-CCT + \ BC BC BC ABC/ + GAA-AAT-ATG-TGC + C/ \ C C C C + GGA GCA-CAT-ATT-TTT + */ + const std::vector sequences { + "TGCCT", + "CGAATGCCT", + "GGAATGCATTTT" + }; + const std::vector labels { "A", "B", "C" }; + + auto anno_graph = build_anno_graph( + k, sequences, labels, DeBruijnGraph::BASIC, true + ); + + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("A"), + anno_graph->get_annotator().get_label_encoder().encode("B")); + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("B"), + anno_graph->get_annotator().get_label_encoder().encode("C")); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + + std::vector>>> exp_alignments { + { std::string("CGAAAGCCTTTT"), { // RC: AAAGGCTTTCG + std::make_tuple(std::string("C"), std::string("GAATGCATTTT"), std::string("1S3=1X2=1X4=")), + std::make_tuple(std::string("B"), std::string("CGAATGCCT"), std::string("4=1X4=3S")), + std::make_tuple(std::string("A"), std::string("GCCT"), std::string("5S4=3S")), + } } + }; + + for (const auto &[query, mappings] : exp_alignments) { + run_alignment(*anno_graph, config, query, mappings, 0, true, true); + } +} + +TYPED_TEST(LabeledAlignerRedoneTest, SimpleTangleGraphCoordsCycle) { + // TODO: for now, not implemented for other annotators + if constexpr(!std::is_same_v> + && !std::is_same_v) { + return; + } + + size_t k = 4; + /* + A A B B B B B + GCAA-CAAT-AATG-ATGC-TGCG-GCGC-CGCA + */ + const std::vector sequences { + "GCAAT", + "AATGCGCA" + }; + const std::vector labels { "A", "B" }; + + auto anno_graph = build_anno_graph( + k, sequences, labels, DeBruijnGraph::BASIC, true + ); + + ASSERT_LT(anno_graph->get_annotator().get_label_encoder().encode("A"), + anno_graph->get_annotator().get_label_encoder().encode("B")); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + + std::vector>>> exp_alignments { + { std::string("ATGCGCAATGCG"), { // RC: CGCATTGCGCAT + std::make_tuple(std::string("B"), std::string("ATGCGCA"), std::string("7=5S")), + } } + }; + + for (const auto &[query, mappings] : exp_alignments) { + run_alignment(*anno_graph, config, query, mappings, 0, true, true); + } +} + +#if ! _PROTEIN_GRAPH +// Seeds shorter than k on a succinct graph (k = 5, min seed length 3), with +// A = TTACGAAA and B = ACGTTGCA. ACGT and GTTG occur only in B, TACGA only in +// A. ACGT is a suffix only of the dummy source k-mer $ACGT that starts B. ACGC +// occurs nowhere and is shorter than k, so it is seeded past the last full +// k-mer; any partial alignment must still spell part of its label's sequence. +TEST(LabeledAlignerRedoneTest, SubkSeedLabels) { + size_t k = 5; + const std::vector sequences { "TTACGAAA", "ACGTTGCA" }; + const std::vector labels { "A", "B" }; + + for (auto mode : { DeBruijnGraph::BASIC, DeBruijnGraph::CANONICAL, DeBruijnGraph::PRIMARY }) { + auto anno_graph = build_anno_graph>( + k, sequences, labels, mode); + const auto &graph = anno_graph->get_graph(); + const auto &label_encoder = anno_graph->get_annotator().get_label_encoder(); + AnnotationBuffer anno_buffer(graph, anno_graph->get_annotator()); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + config.min_seed_length = 3; + config.max_seed_length = k; + + for (const auto &[query, label] : std::vector> { + { "ACGT", "B" }, { "TACGA", "A" }, { "CGTTG", "B" }, { "ACGC", "" } + }) { + Query aln_query(graph, query); + LabeledSeeder seeder(anno_buffer, aln_query, config); + LabeledExtender extender(anno_buffer, aln_query, config); + std::vector paths; + align_query(aln_query, seeder, extender, + [&](Alignment&& path) { paths.emplace_back(std::move(path)); }); + + if (label.empty()) { + EXPECT_TRUE(paths.empty()) << mode << " " << query; + continue; + } + + ASSERT_FALSE(paths.empty()) << mode << " " << query; + for (const auto &path : paths) { + EXPECT_EQ(query, path.get_spelling()) << mode; + for (auto label_class : path.get_label_classes()) { + ASSERT_NE(Anchor::nannot, label_class) << mode << " " << query; + for (auto column : anno_buffer.get_cached_column_set(label_class)) { + EXPECT_EQ(label, label_encoder.decode(column)) << mode << " " << query; + } + } + } + } + } +} + +// TGCA ends B. Sub-k seeds are extended forward to full k-mers, and on basic +// graphs the reverse complement is not seeded, so the seed is lost. Disabled +// until that is fixed. +TEST(LabeledAlignerRedoneTest, DISABLED_SubkSeedAtSequenceEnd) { + size_t k = 5; + const std::vector sequences { "TTACGAAA", "ACGTTGCA" }; + const std::vector labels { "A", "B" }; + + for (auto mode : { DeBruijnGraph::BASIC }) { + auto anno_graph = build_anno_graph>( + k, sequences, labels, mode); + const auto &graph = anno_graph->get_graph(); + const auto &label_encoder = anno_graph->get_annotator().get_label_encoder(); + AnnotationBuffer anno_buffer(graph, anno_graph->get_annotator()); + + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + config.min_seed_length = 3; + config.max_seed_length = k; + + for (const auto &[query, label] : std::vector> { { "TGCA", "B" } }) { + Query aln_query(graph, query); + LabeledSeeder seeder(anno_buffer, aln_query, config); + LabeledExtender extender(anno_buffer, aln_query, config); + std::vector paths; + align_query(aln_query, seeder, extender, + [&](Alignment&& path) { paths.emplace_back(std::move(path)); }); + + if (label.empty()) { + EXPECT_TRUE(paths.empty()) << mode << " " << query; + continue; + } + + ASSERT_FALSE(paths.empty()) << mode << " " << query; + for (const auto &path : paths) { + EXPECT_EQ(query, path.get_spelling()) << mode; + for (auto label_class : path.get_label_classes()) { + ASSERT_NE(Anchor::nannot, label_class) << mode << " " << query; + for (auto column : anno_buffer.get_cached_column_set(label_class)) { + EXPECT_EQ(label, label_encoder.decode(column)) << mode << " " << query; + } + } + } + } + } +} + +#endif + +} // namespace \ No newline at end of file diff --git a/metagraph/tests/annotation/test_annotated_dbg_helpers.cpp b/metagraph/tests/annotation/test_annotated_dbg_helpers.cpp index 71c22d7e44..ee17300c65 100644 --- a/metagraph/tests/annotation/test_annotated_dbg_helpers.cpp +++ b/metagraph/tests/annotation/test_annotated_dbg_helpers.cpp @@ -237,6 +237,7 @@ template std::unique_ptr build_anno_graph build_anno_graph(uint64_t, const std::vector &, const std::vector&, DeBruijnGraph::Mode, bool, const std::vector&); template std::unique_ptr build_anno_graph(uint64_t, const std::vector &, const std::vector&, DeBruijnGraph::Mode, bool, const std::vector&); template std::unique_ptr build_anno_graph(uint64_t, const std::vector &, const std::vector&, DeBruijnGraph::Mode, bool, const std::vector&); +template std::unique_ptr build_anno_graph(uint64_t, const std::vector &, const std::vector&, DeBruijnGraph::Mode, bool, const std::vector&); template std::unique_ptr build_anno_graph(uint64_t, const std::vector &, const std::vector&, DeBruijnGraph::Mode, bool, const std::vector&); template std::unique_ptr build_anno_graph(uint64_t, const std::vector &, const std::vector&, DeBruijnGraph::Mode, bool, const std::vector&); diff --git a/metagraph/tests/graph/all/test_dbg_traverse.cpp b/metagraph/tests/graph/all/test_dbg_traverse.cpp index b0cc29605f..660f7c083f 100644 --- a/metagraph/tests/graph/all/test_dbg_traverse.cpp +++ b/metagraph/tests/graph/all/test_dbg_traverse.cpp @@ -1,5 +1,6 @@ #include "gtest/gtest.h" +#include #include #include "../../test_helpers.hpp" @@ -26,12 +27,23 @@ TYPED_TEST(DeBruijnGraphTest, traverse_string) { std::string path; graph->traverse( cur_node, - sequence.data() + i + k, sequence.data() + sequence.size(), + sequence.substr(i + k), [&](auto node) { path += graph->get_node_sequence(node).back(); } ); EXPECT_EQ(std::string(sequence.begin() + i + k, sequence.end()), path); + + std::string path_back; + graph->traverse_back( + cur_node, + sequence.substr(0, i), + [&](auto node) { + path_back += graph->get_node_sequence(node)[0]; + } + ); + std::reverse(path_back.begin(), path_back.end()); + EXPECT_EQ(std::string(sequence.begin(), sequence.begin() + i), path_back); } } } @@ -44,8 +56,16 @@ TYPED_TEST(DeBruijnGraphTest, traverse_string_stop_when_no_edge) { std::string query = "CCCTGTTTG"; graph->traverse( graph->kmer_to_node("AGGC"), - query.data() + 4, - query.data() + query.size(), + query.substr(4), + [&](auto node) { + EXPECT_FALSE(true) << node << " " << graph->get_node_sequence(node); + }, + []() { return false; } + ); + + graph->traverse_back( + graph->kmer_to_node("TTTG"), + query.substr(0, query.size() - 1), [&](auto node) { EXPECT_FALSE(true) << node << " " << graph->get_node_sequence(node); }, diff --git a/metagraph/tests/graph/test_aligner_redone.cpp b/metagraph/tests/graph/test_aligner_redone.cpp new file mode 100644 index 0000000000..d798f5b67f --- /dev/null +++ b/metagraph/tests/graph/test_aligner_redone.cpp @@ -0,0 +1,1053 @@ +#include + +#include "all/test_dbg_helpers.hpp" +#include "../test_helpers.hpp" + +#include "graph/representation/canonical_dbg.hpp" +#include "seq_io/sequence_io.hpp" + +#include "common/seq_tools/reverse_complement.hpp" +#include "kmer/alphabets.hpp" + +#include "graph/alignment_redone/aln_query.hpp" +#include "graph/alignment_redone/aln_match.hpp" +#include "graph/alignment_redone/aligner_config.hpp" +#include "graph/alignment_redone/aln_seeder.hpp" + +namespace { + +using namespace mtg; +using namespace mtg::graph; +using namespace mtg::graph::align_redone; +using namespace mtg::test; +using namespace mtg::kmer; + +const std::string test_data_dir = "../tests/data"; + +template +class DBGAlignerRedoneTest : public DeBruijnGraphTest {}; + +typedef ::testing::Types AlignGraphTypes; +TYPED_TEST_SUITE(DBGAlignerRedoneTest, AlignGraphTypes); + +void run_alignment(const DeBruijnGraph &graph, + DBGAlignerConfig config, + std::string_view query, + const std::vector &reference, + const std::vector &cigar_str, + size_t end_trim = 0, + bool needs_extension = false, + bool needs_extension_long_seed = false) { + size_t k = graph.get_k(); + if (config.min_seed_length == 0) + config.min_seed_length = k; + + for (auto mx : { k, std::numeric_limits::max() }) { + config.max_seed_length = std::max(mx, config.min_seed_length); + bool check_chaining = mx == std::numeric_limits::max() + ? !needs_extension_long_seed + : !needs_extension; + + Query aln_query(graph, query); + ExactSeeder seeder(aln_query, config); + std::vector paths; + Extender extender(aln_query, config); + std::vector paths_no_extend = seeder.get_inexact_anchors(); + std::vector paths_no_align = seeder.get_inexact_anchors(false); + for (const auto &base_path : paths_no_extend) { + extender.extend(base_path, [&](Alignment&& path) { + paths.emplace_back(std::move(path)); + }); + } + + auto aln_sort = [&config](const auto &a, const auto &b) { + return std::make_pair(score_match(a, config), b.get_orientation()) + > std::make_pair(score_match(b, config), a.get_orientation()); + }; + + std::sort(paths_no_extend.begin(), paths_no_extend.end(), aln_sort); + std::sort(paths_no_align.begin(), paths_no_align.end(), aln_sort); + std::sort(paths.begin(), paths.end(), aln_sort); + + ASSERT_LE(reference.size(), paths.size()) << mx; + paths.resize(reference.size()); + + for (size_t i = 0; i < reference.size(); ++i) { + auto check_ref = [&](const Alignment &path, const std::string &reference, const std::string &type) { + if (reference.size()) { + EXPECT_EQ(reference.size() - k + 1 + end_trim, path.get_path().size()) << mx << "\t" << type; + EXPECT_EQ(reference, path.get_spelling()) << mx << "\t" << type; + } + }; + + auto check_aln = [&](const Alignment &path, const Cigar &cigar, const std::string &reference, const std::string &type) { + ASSERT_EQ(cigar.to_string(), path.get_cigar().to_string()) << mx << "\t" << type; + EXPECT_EQ(end_trim, path.get_end_trim()) << mx << "\t" << type; + if (reference.size()) { + ASSERT_TRUE(cigar.is_valid(reference, path.get_query())); + EXPECT_EQ(config.score_cigar(reference, path.get_query(), cigar), + score_match(path, config)) << mx << "\t" << type; + } + + EXPECT_EQ(cigar.get_clipping(), path.get_clipping()) << mx << "\t" << type; + EXPECT_EQ(cigar.get_end_clipping(), path.get_end_clipping()) << mx << "\t" << type; + }; + + std::string cur_reference = reference[i]; + if (paths[i].get_orientation() && !paths[i].get_end_trim()) + ::reverse_complement(cur_reference.begin(), cur_reference.end()); + + check_ref(paths[i], cur_reference, "extend"); + + if (i < cigar_str.size() && cigar_str[i].size()) { + std::string cur_cigar = cigar_str[i]; + Cigar cigar(cur_cigar); + if (paths[i].get_orientation() && !paths[i].get_end_trim()) { + std::reverse(cigar.data().begin(), cigar.data().end()); + cur_cigar = cigar.to_string(); + } + check_aln(paths[i], cigar, cur_reference, "extend"); + + if (check_chaining) { + // this alignment should work with chaining alone + ASSERT_LT(i, paths_no_extend.size()); + check_ref(paths_no_extend[i], cur_reference, "chain"); + check_aln(paths_no_extend[i], cigar, cur_reference, "chain"); + + ASSERT_LT(i, paths_no_align.size()); + EXPECT_EQ(cigar.to_string(), paths_no_align[i].get_cigar().to_string()) << mx << "\tnoalign"; + } + } + } + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_empty) { + size_t k = 4; + std::string reference = "CATTT"; + std::string query; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, {}, {}); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_sequence_much_too_short) { + size_t k = 4; + std::string reference = "CATTT"; + std::string query = "CA"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, {}, {}); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_sequence_too_short) { + size_t k = 4; + std::string reference = "CATTT"; + std::string query = "CAT"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, {}, {}); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_big_self_loop) { + size_t k = 3; + std::string reference = "AAAA"; + std::string query = "AAAAAAAAA"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { query }, { "9=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_single_node) { + size_t k = 3; + std::string reference = "CAT"; + std::string query = "CAT"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { query }, { "3=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_straight) { + size_t k = 4; + std::string reference = "AGCTTCGAGGCCAA"; + std::string query = reference; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { query }, { "14=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_straight_with_N) { + size_t k = 5; + std::string reference = "TAGCTTCGAGGCCAA"; + std::string query = "TAGCTNCGAGGCCAA"; + // X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + + run_alignment(*graph, config, query, { reference }, { "5=1X9=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 3; + run_alignment(*graph, config, query, { reference }, { "5=1X9=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_ending_branch) { + size_t k = 4; + std::string reference_1 = "AGCTTCGAA"; + std::string reference_2 = "AGCTTCGAC"; + // X + + std::string query = reference_2; + + auto graph = build_graph_batch(k, { reference_1, reference_2 }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { query }, { "9=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_branch) { + size_t k = 6; + std::string reference_1 = "AGCTTCGAATATTTGTT"; + std::string reference_2 = "AGCTTCGACGATTTGTT"; + // XX + + std::string query = reference_2; + + auto graph = build_graph_batch(k, { reference_1, reference_2 }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { query }, { "17=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_branch_with_cycle) { + size_t k = 4; + std::string reference_1 = "AGCTTCGAATATTTGTT"; + std::string reference_2 = "AGCTTCGACGATTTGTT"; + // XX + + std::string query = reference_2; + + auto graph = build_graph_batch(k, { reference_1, reference_2 }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { query }, { "17=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, repetitive_sequence_alignment) { + size_t k = 3; + std::string reference = "AGGGGGGGGGAAAAGGGGGGG"; + std::string query = "AGGGGG"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { query }, { "6=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, variation) { + size_t k = 5; + std::string reference = "AGCAACTCGAAA"; + std::string query = "AGCAATTCGAAA"; + // X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference }, { "5=1X6=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 3; + run_alignment(*graph, config, query, { reference }, { "5=1X6=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, variation_in_branching_point) { + size_t k = 4; + std::string reference_1 = "TTAAGCAACTCGAAA"; + std::string reference_2 = "TTAAGCAAGTCGAAA"; + std::string query = "TTAAGCAATGGGAAA"; + // XXX + + auto graph = build_graph_batch(k, { reference_1, reference_2 }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -2, -2); + config.gap_opening_penalty = -4; + config.gap_extension_penalty = -2; + + run_alignment(*graph, config, query, { "" }, { "8=3X4=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 2; + run_alignment(*graph, config, query, { "" }, { "8=3X4=" }, 0, true, true); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, multiple_variations) { + size_t k = 6; + std::string reference = "TACGCAACTCTCTGAACTTGT"; + std::string query = "TACGCAATTCTCTGTATTTGT"; + // X X X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference }, { "7=1X6=1X1=1X4=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 4; + run_alignment(*graph, config, query, { reference }, { "7=1X6=1X1=1X4=" }, 0, true, true); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_noise_in_branching_point) { + size_t k = 5; + std::string reference_1 = "AAAAACTTTTTT"; + std::string reference_2 = "AAAAATTGGGGG"; + std::string query = "AAAAATTTTTTT"; + // D + + auto graph = build_graph_batch(k, { reference_1, reference_2 }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -2, -2); + + config.gap_opening_penalty = -1; + config.gap_extension_penalty = -1; + run_alignment(*graph, config, query, { "AAAAACTTTTTTT" }, { "5=1D7=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 3; + run_alignment(*graph, config, query, { "AAAAACTTTTTTT" }, { "5=1D7=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, alternative_path_basic) { + size_t k = 5; + std::vector references = { + "AGACAATTTTTTAAA", + "AGACAATTTTTGAAA", + "AGACAAGTTTTTAAA", + "AGACAAGTTTTGAAA" }; + std::string query = "AGACAACTTTTCAAA"; + // X X + + auto graph = build_graph_batch(k, references); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + + run_alignment(*graph, config, query, { "" }, { "6=1X4=1X3=" }, 0, true, true); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_multiple_misalignment) { + size_t k = 5; + std::string reference = "AAAAGCGGACCCTTTCCGTTAT"; + std::string query = "AAAAGGGGACCCTTTTCGTTAT"; + // X X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference }, { "5=1X9=1X6=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 3; + run_alignment(*graph, config, query, { reference }, { "5=1X9=1X6=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_insert_non_existent) { + size_t k = 4; + std::string reference = "TTTCCTTGTT"; + std::string query = "TTTCCATTGTT"; + // I + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + config.gap_opening_penalty = -3; + config.gap_extension_penalty = -3; + run_alignment(*graph, config, query, { reference }, { "5=1I5=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 2; + run_alignment(*graph, config, query, { reference }, { "5=1I5=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_insert_multi) { + size_t k = 4; + std::string reference = "TTTCCTTGTT"; + std::string query = "TTTCCAATTGTT"; + // II + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + config.gap_opening_penalty = -3; + config.gap_extension_penalty = -3; + run_alignment(*graph, config, query, { reference }, { "5=2I5=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 2; + run_alignment(*graph, config, query, { reference }, { "5=2I5=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_insert_long) { + size_t k = 4; + std::string reference = "TTTCCTTGTT"; + std::string query = "TTTCCAAAAAAAAATTGTT"; + // IIIIIIIII + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + config.gap_opening_penalty = -1; + config.gap_extension_penalty = -1; + run_alignment(*graph, config, query, { reference }, { "5=9I5=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 2; + run_alignment(*graph, config, query, { reference }, { "5=9I5=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_insert_long_offset) { + size_t k = 5; + std::string reference = "TTTCCGGTTGTTA"; + std::string query = "TTTCCGCAAAAAAAAATTGTTA"; + // XIIIIIIIII + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + config.gap_opening_penalty = -1; + config.gap_extension_penalty = -1; + run_alignment(*graph, config, query, { reference }, { "6=1X9I6=" }, 0, true, true); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_delete) { + size_t k = 5; + // D + std::string reference = "TTCGAT""A""GGCCT"; + std::string query = "TTCGAT" "GGCCT"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + config.gap_opening_penalty = -3; + config.gap_extension_penalty = -3; + run_alignment(*graph, config, query, { reference }, { "6=1D5=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 3; + run_alignment(*graph, config, query, { reference }, { "6=1D5=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_gap) { + size_t k = 10; + // DDDD + std::string reference = "TTTTCTGTATACCTTGGCGCTCTCC"; + std::string query = "TTTTCTGTATA" "GGCGCTCTCC"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + config.gap_opening_penalty = -3; + config.gap_extension_penalty = -3; + run_alignment(*graph, config, query, { reference }, { "11=4D10=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 5; + run_alignment(*graph, config, query, { reference }, { "11=4D10=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_gap_after_seed) { + size_t k = 4; + // DDDD + std::string reference = "TTTCCCAAGGCGCTCTC"; + std::string query = "TTTC" "GGCGCTCTC"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -2, -2); + config.gap_opening_penalty = -4; + config.gap_extension_penalty = -2; + config.left_end_bonus = 4; + config.right_end_bonus = 4; + run_alignment(*graph, config, query, { reference }, { "4=4D9=" }, 0, true, true); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_loop_deletion) { + size_t k = 4; + std::string reference = "AAAATTTTCGAGGCCAA"; + std::string query = "AAAACGAGGCCAA"; + // DDDD + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::unit_scoring_matrix( + 2, alphabets::kAlphabetDNA, alphabets::kCharToDNA + ); + config.gap_opening_penalty = -2; + config.gap_extension_penalty = -2; + run_alignment(*graph, config, query, { "AAAATTTCGAGGCCAA" }, { "4=3D9=" }, 0, true, true); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_straight_long_xdrop) { + size_t k = 4; + std::string reference_1 = "AGCTTCGAGGCCAAGCCTGACTGATCGATGCATGCTAGCTAGTCAGTCAGCGTGAGCTAGCAT"; + std::string reference_2 = "AGCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT"; + std::string query = reference_1; + + auto graph = build_graph_batch(k, { reference_1, reference_2 }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + config.xdrop = 30; + config.rel_score_cutoff = 0.8; + run_alignment(*graph, config, query, { query }, { "63=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_drop_seed) { + size_t k = 4; + std::string reference = "TTTCCCTGGCGCTCTC"; + std::string query = "TTTCCGGGGCGCTCTC"; + // SSSSSSS + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -6, -6); + config.gap_opening_penalty = -10; + config.gap_extension_penalty = -4; + config.xdrop = 6; + run_alignment(*graph, config, query, { reference.substr(7) }, { "7S9=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_long_gap_after_seed) { + size_t k = 4; + std::string reference = "TTTCCCTTAAGGCGCTCTC"; + std::string query = "TTTC" "GGCGCTCTC"; + // SSSS + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + config.gap_opening_penalty = -5; + config.gap_extension_penalty = -1; + run_alignment(*graph, config, query, { reference.substr(10) }, { "4S9=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_repeat_sequence_no_delete_after_insert) { + size_t k = 27; + std::string reference = "TTTGTGGCTAGAGCTCGAGATCGCGCGGCCACAATTGACAAATGA" "GATCTAATTAAACTAAAGAGCTTCTGCACAGCAAAAGAAACTGTCATC"; + std::string query = "TTTGTGGCTAGAGCTCGAGATCGCGCGGCCACAATTGACAAATGACAAATGTGATCTAATGAAACTAAAGAGCTTCTGCACAGCAAAAGAAACTGTCATC"; + // IIIIIII X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + config.gap_opening_penalty = -3; + config.gap_extension_penalty = -3; + run_alignment(*graph, config, query, { reference }, { "45=7I8=1X39=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 25; + run_alignment(*graph, config, query, { reference }, { "45=7I8=1X39=" }, 0, true, true); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_clipping1) { + size_t k = 4; + std::string reference = "GGCCTGTTTG"; + std::string query = "ACCCTGTTTG"; + // SS + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference.substr(2) }, { "2S8=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_clipping2) { + size_t k = 4; + std::string reference = "AAAAGCTTCGAGGCCAA"; + std::string query = "TTAGCTTCGAGGCCAA"; + // SS + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference.substr(3) }, { "2S14=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_long_clipping) { + size_t k = 4; + std::string reference = "TTTTTTTAAAAGCTTCGAGGCCAA"; + std::string query = "CCCCCCCAAAAGCTTCGAGGCCAA"; + // SSSSSSS + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference.substr(7) }, { "7S17=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_end_clipping) { + size_t k = 4; + std::string reference = "AAAAGCTTCGAGGCCAATTTTTTT"; + std::string query = "AAAAGCTTCGAGGCCAACCCCCCC"; + // SSSSSSS + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference.substr(0, 17) }, { "17=7S" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_low_similarity) { + size_t k = 27; + std::string reference = "CTAGAACTTAAAGTATAATAATACTAATAATAAAATAAAATACA"; + std::string query = "CTAGAACTTAAAGTATAATAATACTAATAA""AAGTACAATACA"; + // DD X X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + run_alignment(*graph, config, query, { reference }, { "30=2D2=1X2=1X6=" }, 0, true, true); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_low_similarity2) { + size_t k = 27; + std::string reference = "GCCACAATTGACAAATGA" "GATCTAATTAAACTAAAGAGCTTCTGCACAGCAAAAGAAACTGTCATC"; + std::string query = "GCCACAATTGACAAATGACAAATGTGATCTAATGAAACTAAAGAGCTTCTGCACAGCAAAAGAAACTGTCATC"; + // IIIIIII X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + + run_alignment(*graph, config, query, { reference.substr(3) }, { "10S4=1X9=1X8=1X39=" }, 0, true, true); + + config.gap_opening_penalty = -3; + run_alignment(*graph, config, query, { reference }, { "18=7I8=1X39=" }, 0, true, true); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_low_similarity2_del) { + size_t k = 27; + std::string query = "GCCACAATTGACAAATGA" "GATCTAATTAAACTAAAGAGCTTCTGCACAGCAAAAGAAACTGTCATC"; + std::string reference = "GCCACAATTGACAAATGACAAATGTGATCTAATGAAACTAAAGAGCTTCTGCACAGCAAAAGAAACTGTCATC"; + // DDDDDDD X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + + config.gap_opening_penalty = -5; + run_alignment(*graph, config, query, { reference.substr(10) }, { "3S4=1X9=1X8=1X39=" }, 0, true, true); + + config.gap_opening_penalty = -3; + config.gap_extension_penalty = -1; + run_alignment(*graph, config, query, { reference }, { "18=7D8=1X39=" }, 0, true, true); +} + +// A repetitive query in a graph of transcripts. The node budget keeps its +// alignment from running for minutes, and a sequence in the graph must still +// align exactly. +void check_low_similarity4(const DeBruijnGraph &graph) { + size_t k = graph.get_k(); + std::string query = "TCGATCGATCGATCGATCGATCGACGATCGATCGATCGATCGATCGACGATCGAT" + "CGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGA" + "TCGATCGATCGATCGACGATCGATCGATCGATCGATCGACGATCGATCGATCGAT" + "CGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGA" + "TCGATCGACGATCGATCGATCGATCGATCGACGATCGATCGATCGATCGATCGAT" + "CGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGA" + "CGATCGATCGATCGATCGATCGACGATCGATCGATCGATCGATCGATCGATCGAT" + "CGATCGATCGATCGATCGATCGA"; + std::string match = "TCGATCAATCGATCAATCGATCAACGATCAATCGATCAATCGATCAACGATCAAT" + "CGATCAATCGATCAATCGATCAATCGATCAATCGATCAATCGATCAATCGATCAA" + "TCGATCAATCGATCAACGATCAATCGATCAATCGATCAACGATCAATCGATCAAT" + "CGATCAATCGATCAATCGATCAATCGATCAATCGATCAATCGATCAATCGATCAA" + "TCGATCAACGATCAATCGATCAATCGATCAACGATCAATCGATCAATCGATCAAT" + "CGATCAATCGATCAATCGATCAATCGATCAATCGATCAATCGATCAATCGATCAA" + "CGATCAATCGATCAATCGATCAACGATCAATCGATCAATCGATCAATCGATCAAT" + "CGATCAATCGATCAATCGATC"; + + ASSERT_TRUE(graph.find(match, 1.0)); + + for (double nodes_per_seq_char : { 10.0, 50.0 }) { + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + config.gap_opening_penalty = -5; + config.gap_extension_penalty = -2; + config.max_nodes_per_seq_char = nodes_per_seq_char; + config.min_seed_length = k; + config.max_seed_length = k; + + Query aln_query(graph, query); + ExactSeeder seeder(aln_query, config); + Extender extender(aln_query, config); + align_query(aln_query, seeder, extender, [](Alignment&&) {}); + + run_alignment(graph, config, match, { match }, + { std::to_string(match.size()) + "=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_low_similarity4) { + size_t k = 6; + std::vector seqs; + mtg::seq_io::read_fasta_file_critical(test_data_dir + "/transcripts_100.fa", + [&](auto *seq) { seqs.emplace_back(seq->seq.s); }); + check_low_similarity4(*build_graph_batch(k, std::move(seqs))); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_low_similarity5) { + size_t k = 31; + std::string reference = "GTCGTCAGATCGGAAGAGCGTCGTGTAGGGAAAG" "GTCTTC""GCCTGTGTAGATCTC"; + std::string query = "GTCAGATCGGAAGAGCGTCGTGTAGGGAAAG""AGTGTTCCTGGTGGTGTAGATCTC"; + // I X II XXX + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -2, -2); + config.left_end_bonus = 5; + config.right_end_bonus = 5; + run_alignment(*graph, config, query, {}, {}); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 11; + run_alignment(*graph, config, query, { reference.substr(3, 52) }, { "31=1I2=1X3=2I1=3X11=" }, 0, true, true); + } +} + +#if ! _PROTEIN_GRAPH + +TYPED_TEST(DBGAlignerRedoneTest, align_low_similarity4_rep_primary) { + if constexpr(std::is_same_v) { + // primary SSHash graphs need an odd k + return; + } + + size_t k = 6; + std::vector seqs; + mtg::seq_io::read_fasta_file_critical(test_data_dir + "/transcripts_100.fa", + [&](auto *seq) { seqs.emplace_back(seq->seq.s); }); + check_low_similarity4(*build_graph_batch(k, std::move(seqs), DeBruijnGraph::PRIMARY)); +} +#endif + +#if ! _PROTEIN_GRAPH +TYPED_TEST(DBGAlignerRedoneTest, align_both_directions) { + size_t k = 7; + std::string reference = "AAAAGCTTTCGAGGCCAA"; + std::string query = "AAAAGTTTTCGAGGCCAA"; + // X + + auto graph = build_graph_batch(k, { reference }, DeBruijnGraph::CANONICAL); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference }, { "5=1X12=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 5; + run_alignment(*graph, config, query, { reference }, { "5=1X12=" }); + } +} +#endif + +TYPED_TEST(DBGAlignerRedoneTest, align_both_directions2) { + size_t k = 11; + std::string reference = "GTAGTGCTAGCTTGTAGTCGTGCTGATGCA"; + std::string query = "GTAGTGCTACCTTGTAGTCGTGGTGATGCA"; + // X X + + auto graph = build_graph_batch(k, { reference }, DeBruijnGraph::BASIC); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference.substr(0, query.size()) }, { "9=1X12=1X7=" }, 0, true, true); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_seed_to_end) { + size_t k = 5; + // SSSSSSSSSSSSSSSSSSSSSS + std::string reference = "ATCCCTTTTAAAA"; + std::string query = "ATCCCGGGGGGGGGGGGGGGGGTTTTAAAA"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + run_alignment(*graph, config, query, { reference.substr(5) }, { "22S8=" }); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_bfs_vs_dfs_xdrop) { + size_t k = 31; + std::string reference_1 = "TCGGGGCAAGAAACACACAGCCTTCTCATCCAAGGGCCTCAGTGATGAAGAGTACGATGAGTACAAGAGGATCAGAGAAGAAAGGAATGGCAAATACTCCATAGAAGAGTACCTTCAGGACAGGGACAGATACTATGAGGAGGTGGCCAT"; + std::string reference_2 = "TCGGGGCAAGAAACACACAGCCTTCTCATCCAAGGGCCTCAGTGATGAAGAGTACGATGAGTACAAGAGAATCAGAGAGGAGAGGAATGGCAAATACTCAATAGAGGAATACCTCCAAGATAGGGACAGATACTATGAAGAGCTTGCCAT"; + std::string query = "TCGGGGCAAGAAACACACAGCCTTCTCATCCAAGGGCCTCAGTGATGATGAGTACGATGAGTACAAGAGCATCAGAGAGGAGAGGAATGGCAAATACTCAATAGAGGAATACCTCCAAGATAGGGACAGATACTATGAAGAGCTTGCCAT"; + + auto graph = build_graph_batch(k, { reference_1, reference_2 }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + config.xdrop = 27; + run_alignment(*graph, config, query, { "" }, { "48=1X20=1X80=" }, 0, true, true); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 19; + run_alignment(*graph, config, query, { "" }, { "48=1X20=1X80=" }); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_dummy_short) { + size_t k = 7; + std::string reference = "GAAAAGCT"; + std::string query = "GAAAAGTT"; + // X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + + run_alignment(*graph, config, query, {}, {}); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 5; + run_alignment(*graph, config, query, { reference }, { "6=1X1=" }, 0, true, true); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_dummy_short2) { + size_t k = 31; + std::string reference = "GCTAGTCTATAGCGTACGATCTAGCTAGTCT"; + std::string query = "GCTAGTCTATAGCGTACGATCTTTTTAGTC"; + // XXX + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -1); + + run_alignment(*graph, config, query, {}, {}); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 19; + run_alignment(*graph, config, query, { reference.substr(0, 30) }, { "22=3X5=" }, 1, true, true); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_dummy) { + size_t k = 7; + std::string reference = "AAAAGCTTTCGATT"; + std::string query = "AAAAGTTTTCGAA"; + // X X + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + config.left_end_bonus = 2; + config.right_end_bonus = 2; + run_alignment(*graph, config, query, {}, {}); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 5; + run_alignment(*graph, config, query, { reference.substr(0, query.size()) }, { "5=1X6=1X" }, 0, true, true); + } +} + +TYPED_TEST(DBGAlignerRedoneTest, align_extended_insert_after_match) { + size_t k = 27; + std::string reference_1 = "CGTGGCCCAGGCCCAGGCCCAG" "GCCCAGGCCCAGGC""CCA" "GGCCCAGGC""CC" "AGGCCCAGGCCCAGGCCCAAGCC"; + std::string reference_2 = "CGTGGCCCAGGCCCAGGCCCAG" "CCCCAGGCCCAGGC""CCA" "GGCCCAGGC""CC" "AGGCCCAGGCCCAGGCCCAAGCC"; + std::string query = "CGTGGCCCAGGCCCAGGCCCAGTGGGCGTTGGCCCAGGCGGCCA""CGG" "TGGCTGCG""CAGGCCCGCCTGGCACAAGCCACGCTG"; + // III XXX II I DDDX II X I + + auto graph = build_graph_batch(k, { reference_1, reference_2 }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + run_alignment(*graph, config, query, {}, {}); + + if constexpr(std::is_base_of_v) { + config.min_seed_length = 15; + run_alignment(*graph, config, query, + { "CGTGGCCCAGGCCCAGGCCCAGGCCCAGGCCCAGGCCCAGGCCCAGGCCCAGGCCCAGGCCCAAGCC" }, + { "" }, 0, true, true); + } +} + +#if ! _PROTEIN_GRAPH +TYPED_TEST(DBGAlignerRedoneTest, align_suffix_seed_no_full_seeds) { + size_t k = 31; + std::string reference = "CTGCTGCGCCATCGCAACCCACGGTTGCTTTTTGAGTCGCTGCTCACGTTAGCCATCACACTGACGTTAAGCTGGCTTTCGATGCTGTATCTTTTTTTT"; + std::string query = "CTTACTGCTGCGCTCTTCGCAAACCCCACGGTTTCTTGTTTTGAGCTCGCCTGCTCACGATACCCATACACACTGACGTTCAAGCTGGCTTTCGATGTTGTATC"; + // SSSS + + auto graph = build_graph_batch(k, { reference }, DeBruijnGraph::PRIMARY); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -3, -3); + run_alignment(*graph, config, query, {}, {}); + if constexpr(std::is_base_of_v) { + config.min_seed_length = 9; + run_alignment(*graph, config, query, + { reference.substr(0, 91) }, + { "4S9=1I1=1X6=2I9=1X2=2I7=1I3=1I9=1X2=1X4=1I12=1I16=1X6=" }, 0, true, true); + } +} +#endif + +// Records the connect flag that align_query passes to the seeder +class RecordingSeeder : public ExactSeeder { + public: + using ExactSeeder::ExactSeeder; + + std::vector get_inexact_anchors(bool align = true) const override { + flags.push_back(align); + return ExactSeeder::get_inexact_anchors(align); + } + + mutable std::vector flags; +}; + +TYPED_TEST(DBGAlignerRedoneTest, align_query_forwards_connect_anchors) { + size_t k = 10; + std::string reference = "TTTTCTGTATACCTTGGCGCTCTCC"; + std::string query = "TTTTCTGTATA" "GGCGCTCTCC"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + config.gap_opening_penalty = -3; + config.gap_extension_penalty = -3; + config.min_seed_length = k; + config.max_seed_length = k; + Query aln_query(*graph, query); + Extender extender(aln_query, config); + + for (bool connect : { true, false }) { + RecordingSeeder seeder(aln_query, config); + std::vector paths; + auto callback = [&](Alignment&& path) { paths.emplace_back(std::move(path)); }; + + align_query(aln_query, seeder, extender, callback, connect); + EXPECT_EQ(std::vector{ connect }, seeder.flags) << "extend"; + EXPECT_FALSE(paths.empty()) << connect; + + seeder.flags.clear(); + align_query(aln_query, seeder, callback, connect); + EXPECT_EQ(std::vector{ connect }, seeder.flags) << "no extend"; + } +} + +TYPED_TEST(DBGAlignerRedoneTest, extend_unconnected_chain) { + size_t k = 5; + std::string reference = "AGCAACTCGAAA"; + std::string query = "AGCAATTCGAAA"; + + auto graph = build_graph_batch(k, { reference }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + // Anchors are chained across the mismatch with sub-k seeds (DBGSuccinct) + // or by SSHash's own chaining + config.min_seed_length = std::is_base_of_v ? 3 : k; + config.max_seed_length = k; + Query aln_query(*graph, query); + ExactSeeder seeder(aln_query, config); + Extender extender(aln_query, config); + + bool has_placeholder = false; + for (const auto &chain : seeder.get_inexact_anchors(false)) { + const auto &path = chain.get_path(); + if (std::find(path.begin(), path.end(), DeBruijnGraph::npos) == path.end()) + continue; + + has_placeholder = true; + std::vector paths; + extender.extend(chain, [&](Alignment&& path) { paths.emplace_back(std::move(path)); }); + ASSERT_EQ(1u, paths.size()); + EXPECT_EQ(chain.get_cigar().to_string(), paths[0].get_cigar().to_string()); + EXPECT_EQ(chain.get_path(), paths[0].get_path()); + } + if constexpr(!std::is_same_v) { + EXPECT_TRUE(has_placeholder); + } +} + +#if ! _PROTEIN_GRAPH +TYPED_TEST(DBGAlignerRedoneTest, query_reverse_complement) { + auto graph = build_graph_batch(4, { "AGCTTCGAGGCCAA" }); + Query query(*graph, std::string("AACCG")); + EXPECT_EQ("AACCG", query.get_query(false)); + EXPECT_EQ("CGGTT", query.get_query(true)); +} + +TYPED_TEST(DBGAlignerRedoneTest, align_reverse_complement_basic) { + size_t k = 5; + std::string reference = "AGCTTCGAGGCCAATCGT"; + std::string query = reference; + ::reverse_complement(query.begin(), query.end()); + + auto graph = build_graph_batch(k, { reference }, DeBruijnGraph::BASIC); + ASSERT_EQ(DeBruijnGraph::BASIC, graph->get_mode()); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + config.min_seed_length = k; + config.max_seed_length = k; + Query aln_query(*graph, query); + + ExactSeeder seeder(aln_query, config); + Extender extender(aln_query, config); + std::vector paths; + align_query(aln_query, seeder, extender, + [&](Alignment&& path) { paths.emplace_back(std::move(path)); }); + ASSERT_FALSE(paths.empty()); + auto best = std::max_element(paths.begin(), paths.end(), + [&](const auto &a, const auto &b) { + return score_match(a, config) < score_match(b, config); + }); + EXPECT_TRUE(best->get_orientation()); + EXPECT_EQ(reference, best->get_spelling()); + EXPECT_EQ("18=", best->get_cigar().to_string()); +} +#endif + +TYPED_TEST(DBGAlignerRedoneTest, max_num_seeds_per_locus) { + if constexpr(!std::is_base_of_v) { + return; + } + + size_t k = 5; + // ACG continues in four ways, so the query ACGAT has four seeds shorter + // than k at its first position + auto graph = build_graph_batch(k, { "ACGAA", "ACGCC", "ACGTT", "ACGGG" }); + DBGAlignerConfig config; + config.score_matrix = DBGAlignerConfig::dna_scoring_matrix(2, -1, -2); + config.min_seed_length = 3; + config.max_seed_length = k; + Query aln_query(*graph, std::string("ACGAT")); + + auto num_short_seeds_at_start = [&]() { + auto anchors = ExactSeeder(aln_query, config).get_anchors(); + return std::count_if(anchors.begin(), anchors.end(), [&](const Anchor &a) { + return !a.get_orientation() && !a.get_clipping() && a.get_seed().size() < k; + }); + }; + + ASSERT_LT(1, num_short_seeds_at_start()); + config.max_num_seeds_per_locus = 1; + EXPECT_EQ(0, num_short_seeds_at_start()); +} + +} // namespace