Preserve mapping attribution and run NumPy correctness checks in CI

.github/workflows/CI.yml

@@ -19,6 +19,27 @@ permissions:
   contents: read
 
 jobs:
+  correctness:
+    name: Python correctness (${{ matrix.python-version }})
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        python-version: ["3.12", "3.13", "3.14"]
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+      - uses: dtolnay/rust-toolchain@stable
+      - name: Install development dependencies
+        run: python -m pip install maturin numpy pytest
+      - name: Build and install the local package
+        run: python -m pip install .
+      - name: Check the installed wheel's public NumPy results
+        working-directory: ${{ runner.temp }}
+        run: python -m pytest --import-mode=importlib -q "$GITHUB_WORKSPACE/tests"
+
   linux:
     runs-on: ${{ matrix.platform.runner }}
     strategy:
@@ -157,7 +178,7 @@ jobs:
     name: Release
     runs-on: ubuntu-latest
     if: ${{ startsWith(github.ref, 'refs/tags/') || github.event_name == 'workflow_dispatch' }}
-    needs: [linux, musllinux, windows, macos, sdist]
+    needs: [correctness, linux, musllinux, windows, macos, sdist]
     permissions:
       # Use to sign the release artifacts
       id-token: write

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "ruranges-py"
 # Keep this synchronized with `project.version` in pyproject.toml.
-version = "0.2.8"
+version = "0.2.9"
 edition = "2021"
 
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
@@ -19,7 +19,7 @@ python-module = ["pyo3/extension-module"]
 
 [dependencies]
 # Core interval algorithms now live in the separate Rust crate/repo.
-ruranges-core = { version = "0.1.13" }
+ruranges-core = { version = "0.1.16" }
 
 # Python bindings.
 rustc-hash = "2.1.0"

README.md

@@ -280,3 +280,11 @@ Any function that needs strand expects a boolean array: True for the minus stran
 ## License
 
 Apache 2.0. See LICENSE for details.
+
+### Transcript projection
+
+`ruranges.numpy.map_to_global` accepts queries first and exon annotations in
+the `*2` arrays. Inputs can be unsorted: the core sorts them while preserving
+original row indices and clips queries to exon coverage. Pass
+`return_exon_indices=True` to receive a fifth result array naming the original
+annotation row for each projected fragment. The default still returns four arrays.

pyproject.toml

@@ -17,7 +17,7 @@ classifiers = [
 ]
 
 # Keep this synchronized with `package.version` in Cargo.toml.
-version = "0.2.8"
+version = "0.2.9"
 
 [tool.maturin]
 module-name = "ruranges.ruranges"

ruranges/numpy.py

@@ -143,13 +143,16 @@ def map_to_global(
     genome_end2:   NDArray[RangeInt],
     strand2:   NDArray[np.bool_],
     sort_output: bool = True,
-) -> tuple[NDArray, NDArray, NDArray, NDArray]:
+    return_exon_indices: bool = False,
+) -> tuple[NDArray, ...]:
     """
     Vectorised transcript-to-genome projection.
 
-    All arrays must be 1-D. Each table must be sorted by transcript ID and
-    local start. The first table contains queries; the ``*2`` arrays contain
-    the exon annotation. The *strand* arrays encode '+' as True and '−' as False.
+    All arrays must be 1-D. Rust sorts either table when needed. The first
+    table contains queries; the ``*2`` arrays contain non-overlapping local
+    exon spans. Only query/exon intersections are projected. The *strand*
+    arrays encode '+' as True and '−' as False. With ``sort_output=True``,
+    fragments follow original query and annotation row order.
 
     Returns
     -------
@@ -159,6 +162,9 @@ def map_to_global(
         Genomic coordinates of the mapped pieces.
     g_strand_bool : bool
         True ⇒ '+', False ⇒ '−'.
+    exon_idx : uint32, optional
+        When ``return_exon_indices=True``, a fifth array identifies the
+        original annotation row producing each fragment.
     """
 
     return _dispatch_map_global_binary(
@@ -167,6 +173,7 @@ def map_to_global(
         groups,  starts,  ends,  strand,
         chr_code2, genome_start2, genome_end2,
         sort_output=sort_output,
+        return_exon_indices=return_exon_indices,
     )
 
 
@@ -1594,6 +1601,8 @@ def _dispatch_map_global_binary(
     **extra_kw
 ):
     roles = RETURN_SIGNATURES[prefix]
+    if extra_kw.get("return_exon_indices", False):
+        roles += ("index",)
     ex_start = _as_1d_array(ex_local_start, name="starts")
     ex_end = _as_1d_array(ex_local_end, name="ends")
     q_start_arr = _as_1d_array(q_start, name="starts2")

src/bindings/numpy_bindings/map_to_global_numpy.rs

@@ -1,7 +1,7 @@
-use numpy::{IntoPyArray, PyArray1, PyReadonlyArray1};
-use pyo3::prelude::*;
+use numpy::{IntoPyArray, PyReadonlyArray1};
+use pyo3::{prelude::*, IntoPyObjectExt};
 
-use ruranges_core::map_to_global::map_to_global; // core algorithm
+use ruranges_core::map_to_global::map_to_global_with_exons; // core algorithm
 
 /* =======================================================================
    Macro:  expose map_to_global_<suffix>() functions to Python/NumPy
@@ -17,20 +17,7 @@ use ruranges_core::map_to_global::map_to_global; // core algorithm
 macro_rules! define_map_to_global_numpy {
     ($fname:ident, $code_ty:ty, $pos_ty:ty) => {
         #[pyfunction]
-        #[pyo3(signature = (
-            ex_tx,
-            ex_local_start,
-            ex_local_end,
-            q_tx,
-            q_start,
-            q_end,
-            ex_chr_code,
-            ex_genome_start,
-            ex_genome_end,
-            ex_fwd,
-            q_fwd,
-            sort_output = true,
-        ))]
+        #[pyo3(signature = (ex_tx, ex_local_start, ex_local_end, q_tx, q_start, q_end, ex_chr_code, ex_genome_start, ex_genome_end, ex_fwd, q_fwd, sort_output = true, return_exon_indices = false,))]
         #[allow(non_snake_case)]
         pub fn $fname<'py>(
             py: Python<'py>,
@@ -49,13 +36,9 @@ macro_rules! define_map_to_global_numpy {
             ex_fwd: PyReadonlyArray1<bool>,
             q_fwd: PyReadonlyArray1<bool>,
             sort_output: bool,
-        ) -> PyResult<(
-            Py<PyArray1<u32>>,     // indices back into query table
-            Py<PyArray1<$pos_ty>>, // genomic start
-            Py<PyArray1<$pos_ty>>, // genomic end
-            Py<PyArray1<bool>>,    // strand (+ = True)
-        )> {
-            let (idx, g_start, g_end, strand) = map_to_global(
+            return_exon_indices: bool,
+        ) -> PyResult<Py<PyAny>> {
+            let (idx, g_start, g_end, strand, exons) = map_to_global_with_exons(
                 /*  exons first (left triple)  */
                 ex_tx.as_slice()?,
                 ex_local_start.as_slice()?,
@@ -73,12 +56,15 @@ macro_rules! define_map_to_global_numpy {
                 sort_output,
             );
 
-            Ok((
-                idx.into_pyarray(py).to_owned().into(),
-                g_start.into_pyarray(py).to_owned().into(),
-                g_end.into_pyarray(py).to_owned().into(),
-                strand.into_pyarray(py).to_owned().into(),
-            ))
+            let idx = idx.into_pyarray(py);
+            let g_start = g_start.into_pyarray(py);
+            let g_end = g_end.into_pyarray(py);
+            let strand = strand.into_pyarray(py);
+            if return_exon_indices {
+                (idx, g_start, g_end, strand, exons.into_pyarray(py)).into_py_any(py)
+            } else {
+                (idx, g_start, g_end, strand).into_py_any(py)
+            }
         }
     };
 }

tests/test_map_to_global.py

@@ -65,3 +65,49 @@ def view(values, dtype=np.int64):
     np.testing.assert_array_equal(starts, np.array([102, 104, 204]) + shift)
     np.testing.assert_array_equal(ends, np.array([105, 107, 207]) + shift)
     np.testing.assert_array_equal(strand, [True, False, False])
+
+
+@pytest.mark.parametrize("dtype", [np.int32, np.int64])
+@pytest.mark.parametrize("sort_output", [False, True])
+def test_unsorted_projection_retains_original_query_and_exon_indices(dtype, sort_output):
+    out = map_to_global(
+        starts=np.array([5, 5, 0], dtype=dtype),
+        ends=np.array([15, 15, 2], dtype=dtype),
+        groups=np.array([7, 1, 7]), strand=np.array([True, False, True]),
+        starts2=np.array([10, 0, 0, 10], dtype=dtype),
+        ends2=np.array([20, 10, 10, 20], dtype=dtype),
+        groups2=np.array([7, 1, 7, 1]), chr_code2=np.array([1, 1, 1, 1]),
+        genome_start2=np.array([200, 1000, 100, 2000], dtype=dtype),
+        genome_end2=np.array([210, 1010, 110, 2010], dtype=dtype),
+        strand2=np.array([True, True, True, True]),
+        sort_output=sort_output, return_exon_indices=True,
+    )
+    expected = [(0, 200, 205, True, 0), (0, 105, 110, True, 2),
+                (1, 1005, 1010, False, 1), (1, 2000, 2005, False, 3),
+                (2, 100, 102, True, 2)]
+    actual = list(zip(*(a.tolist() for a in out)))
+    assert (actual if sort_output else sorted(actual)) == (expected if sort_output else sorted(expected))
+
+
+@pytest.mark.parametrize("forward", [False, True])
+@pytest.mark.parametrize("query_start,query_end,local_start,local_end", [
+    (-5, 5, 0, 10), (0, 5, 10, 20), (5, 15, 10, 20), (15, 25, 10, 20),
+])
+def test_projection_clips_both_exon_boundaries(forward, query_start, query_end, local_start, local_end):
+    out = map_to_global(
+        starts=np.array([query_start]), ends=np.array([query_end]),
+        groups=np.array([7]), strand=np.array([True]),
+        starts2=np.array([local_start]), ends2=np.array([local_end]),
+        groups2=np.array([7]), chr_code2=np.array([1]),
+        genome_start2=np.array([100]), genome_end2=np.array([110]),
+        strand2=np.array([forward]),
+    )
+    lo, hi = max(query_start, local_start), min(query_end, local_end)
+    if lo >= hi:
+        assert all(a.size == 0 for a in out)
+    else:
+        expected = (100 + lo - local_start, 100 + hi - local_start) if forward else (110 - hi + local_start, 110 - lo + local_start)
+        np.testing.assert_array_equal(out[0], [0])
+        np.testing.assert_array_equal(out[1], [expected[0]])
+        np.testing.assert_array_equal(out[2], [expected[1]])
+        np.testing.assert_array_equal(out[3], [forward])

tests/test_strided_input.py

@@ -76,19 +76,6 @@ def pair(lo):
     "genome_bounds": lambda lo: rn.genome_bounds(
         groups=groups(), starts=starts(lo), ends=ends(lo), chrom_length=strided([lo + 100, lo + 100])
     ),
-    "map_to_global": lambda lo: rn.map_to_global(
-        starts=starts(0),
-        ends=ends(0),
-        groups=groups(),
-        strand=negative_strand(),
-        starts2=starts(0),
-        ends2=ends(0),
-        groups2=groups(),
-        chr_code2=groups(),
-        genome_start2=starts(lo),
-        genome_end2=ends(lo),
-        strand2=negative_strand(),
-    ),
     "sort_groups": lambda lo: rn.sort_groups(groups()),
     "fold_ranks": lambda lo: rn.fold_ranks(groups(), 1, groups(), 1),
 }
@@ -98,3 +85,21 @@ def pair(lo):
 @pytest.mark.parametrize("name", CASES)
 def test_kernel_accepts_strided_input(name, lo):
     CASES[name](lo)
+
+
+@pytest.mark.parametrize("lo", [1, BIG])
+def test_strided_projection_checks_coordinates_and_query_attribution(lo):
+    # One query spans two exons. Exchanging the tables must fail this check.
+    idx, out_starts, out_ends, strand = rn.map_to_global(
+        starts=strided([1, 3]), ends=strided([5, 6]), groups=strided([7, 7], np.uint32),
+        strand=strided([True, True], np.bool_),
+        starts2=strided([0, 2, 6]), ends2=strided([2, 6, 8]),
+        groups2=strided([7, 7, 7], np.uint32), chr_code2=strided([1, 1, 1], np.uint32),
+        genome_start2=strided([lo, lo + 10, lo + 20]),
+        genome_end2=strided([lo + 2, lo + 14, lo + 22]),
+        strand2=strided([True, True, True], np.bool_),
+    )
+    np.testing.assert_array_equal(idx, [0, 0, 1])
+    np.testing.assert_array_equal(out_starts, [lo + 1, lo + 10, lo + 11])
+    np.testing.assert_array_equal(out_ends, [lo + 2, lo + 13, lo + 14])
+    np.testing.assert_array_equal(strand, [True, True, True])