Add Harun's redone aligner as opt-in engine with pymetagraph bindings

.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*

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 \"$<TARGET_FILE:Boost::iostreams>\")
+      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

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.

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"]

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()