Migrate alignment functions to Rust

Cargo.toml

@@ -4,14 +4,15 @@ version = "0.0.0"
 edition = "2018"
 
 [dependencies]
-numpy = "0.27"
+numpy = "0.29"
 ndarray = "0.17"
 num-traits = "0.2"
 smallvec = "1"
 itertools = "0.14"
+bitflags = "2"
 
 [dependencies.pyo3]
-version = "0.27"
+version = "0.29"
 features = ["extension-module"]
 
 [lib]
@@ -22,3 +23,6 @@ path = "src/rust/lib.rs"
 # Enable optimizations in dev mode for better performance during development
 [profile.dev]
 opt-level = 3
+
+[lints.rust]
+warnings = "deny"

benchmarks/sequence/align/benchmark_pairwise.py

@@ -31,6 +31,7 @@ def seed(seq_pair):
 
 
 @pytest.mark.benchmark
+@pytest.mark.parametrize("score_only", [False, True], ids=["traceback", "score_only"])
 @pytest.mark.parametrize(
     "method",
     [
@@ -40,9 +41,9 @@ def seed(seq_pair):
     ],
     ids=lambda x: x.func.__name__,
 )
-def benchmark_align_pairwise(seq_pair, matrix, seed, method):
+def benchmark_align_pairwise(seq_pair, matrix, seed, method, score_only):
     """
     Perform pairwise sequence alignment using different algorithms.
     """
     kwargs = {"seed": seed} if method.func is align.align_local_gapped else {}
-    method(seq_pair[0], seq_pair[1], matrix, **kwargs)
+    method(seq_pair[0], seq_pair[1], matrix, score_only=score_only, **kwargs)

doc/apidoc.json

@@ -151,6 +151,7 @@
         "Aligners" : [
             "align_ungapped",
             "align_optimal",
+            "SeedExtension",
             "align_local_ungapped",
             "align_local_gapped",
             "align_banded",

doc/examples/scripts/sequence/homology/genome_comparison.py

@@ -128,7 +128,7 @@
 
 ########################################################################
 # The next step is a local ungapped alignment at the positions of the
-# double hits using :func:`align_local_ungapped()`:
+# double hits using :class:`SeedExtension`:
 # For each hit, the alignment is extended into both directions from the
 # match until the similarity score drops more than a given threshold
 # below the maximum score found.
@@ -139,19 +139,17 @@
 # As a result, the ``score_only=True`` parameter increases the
 # performance drastically.
 
-
 X_DROP = 20
 ACCEPT_THRESHOLD = 100
 
 matrix = align.SubstitutionMatrix.std_nucleotide_matrix()
+extension = align.SeedExtension(matrix, threshold=X_DROP)
 ungapped_scores = np.array(
     [
-        align.align_local_ungapped(
+        extension.align(
             chloroplast_seq,
             bacterium_seqs[strand],
-            matrix,
             seed=(i, j),
-            threshold=X_DROP,
             score_only=True,
         )
         for i, strand, j in double_hits

doc/examples/scripts/structure/alphabet/structure_search.py

@@ -276,14 +276,13 @@ def filter_undefined_spans(sequence, min_length):
 # We only keep those hits that meet the given score threshold.
 
 substitution_matrix = align.SubstitutionMatrix.std_3di_matrix()
+extension = align.SeedExtension(substitution_matrix, threshold=X_DROP)
 ungapped_scores = np.array(
     [
-        align.align_local_ungapped(
+        extension.align(
             query_sequence,
             db_sequences[db_index],
-            substitution_matrix,
             seed=(i, j),
-            threshold=X_DROP,
             score_only=True,
         )
         for i, db_index, j in double_hits

doc/tutorial/sequence/align_heuristic.rst

@@ -160,12 +160,11 @@ In other words the alignment comprises only a region that is decently similar.
     THRESHOLD = 20
 
     matrix = align.SubstitutionMatrix.std_protein_matrix()
+    # The scoring scheme is preprocessed once and reused for each seed
+    extension = align.SeedExtension(matrix, threshold=THRESHOLD)
     alignments = []
     for query_pos, ref_id, ref_pos in matches:
-        alignment = align.align_local_ungapped(
-            reference, query, matrix,
-            seed=(ref_pos, query_pos), threshold=THRESHOLD
-        )
+        alignment = extension.align(reference, query, seed=(ref_pos, query_pos))
         alignments.append(alignment)
 
     for alignment in alignments:
@@ -190,7 +189,7 @@ The difference between the heuristic gapped sequence alignment methods and
 :meth:`align_optimal()` is that the former only ideally traverses through a
 small fraction of the possible alignment search space, allowing them to run
 much faster.
-However, like :meth:`align_local_ungapped()` they need to be informed with a
+However, like :class:`SeedExtension` they need to be informed with a
 match position to start from.
 Furthermore, in some cases they might not find the optimal alignment, when the
 assumption of the method does not hold for such alignment.

setup.py

@@ -2,8 +2,6 @@
 import shutil
 import sys
 import numpy as np
-from Cython.Build import cythonize
-from puccinialin import setup_rust
 from setuptools import setup
 from setuptools_rust import RustExtension
 
@@ -19,6 +17,9 @@ def _should_build_wheel():
     if not os.environ.get("BIOTITE_OMIT_RUST", False):
         if not shutil.which("cargo"):
             # Rust compiler is not installed -> Install it temporarily
+            # `puccinialin` is only required in this fallback case
+            from puccinialin import setup_rust
+
             extra_env = setup_rust()
             env = {**os.environ, **extra_env}
         else:
@@ -34,6 +35,8 @@ def _should_build_wheel():
         rust_extensions = None
 
     if not os.environ.get("BIOTITE_OMIT_CYTHON", False):
+        from Cython.Build import cythonize
+
         # Only build C files and compile them when building a wheel
         cython_extensions = cythonize(
             "src/**/*.pyx",