Use Rust for transcript-to-genome projection

CHANGELOG.txt

@@ -1,4 +1,5 @@
 # 1.5.2 (27.10.26)
+- Project transcript queries through Rust, retaining original query/exon order, indices, strand composition and optional metadata. Negative offsets retain their transcript-end interpretation. References whose ids span multiple chromosomes or strands retain the pandas path. Requires ruranges >= 0.2.9.
 - Fix the `stats.fisher_exact` docstring example, which passed TN and FN in each other's slots and showed the odds ratio and p-values for the wrong table. Found and fixed by Mohit Arvind Khakharia (@Mohit-Ak) in pyranges/pyranges0#411
 
 # 1.5.1 (25.10.26)

README.md

@@ -23,6 +23,7 @@ available at https://pyranges1.readthedocs.io/
 
 ```
 # 1.5.2 (27.10.26)
+- Project transcript queries through Rust, retaining original query/exon order, indices, strand composition and optional metadata. Negative offsets retain their transcript-end interpretation. References whose ids span multiple chromosomes or strands retain the pandas path. Requires ruranges >= 0.2.9.
 - Fix the `stats.fisher_exact` docstring example, which passed TN and FN in each other's slots and showed the odds ratio and p-values for the wrong table. Found and fixed by Mohit Arvind Khakharia (@Mohit-Ak) in pyranges/pyranges0#411
 
 # 1.5.1 (25.10.26)

docs/how_to_map.rst

@@ -481,3 +481,13 @@ With the similar logic, we can easily map the start and stop codon positions to
   149      |    rna-DGYR_LOCUS12552-2  637      640      +         stop
   PyRanges with 56 rows, 5 columns, and 1 index columns.
   Contains 17 chromosomes and 1 strands.
+
+Implementation and coordinate bounds
+------------------------------------
+
+``map_to_global`` uses the Rust projection engine, which sorts query and exon
+coordinates internally and restores original input order. Negative local
+offsets count from the end of the transcript, and portions outside transcript
+coverage are omitted. Identifiers reused across chromosomes or strands retain
+the pandas implementation because each locus has a separate local coordinate
+system.

pyproject.toml

@@ -30,7 +30,7 @@ classifiers = [
     "Topic :: Scientific/Engineering",
 ]
 keywords = ["bioinformatics", "genomicranges", "genomics"]
-dependencies = ["pandas", "tabulate", "natsort", "ruranges>=0.2.8", "gtfreader>=0.3.1"]
+dependencies = ["pandas", "tabulate", "natsort", "ruranges>=0.2.9", "gtfreader>=0.3.1"]
 
 [project.optional-dependencies]
 # pyarrow parses on every core instead of one: read_bed is 10-18x faster on
@@ -134,7 +134,7 @@ deps =
     ty == 0.0.87
     joblib
     gtfreader >= 0.3.1
-    ruranges >= 0.2.8
+    ruranges >= 0.2.9
     # The readers run every test twice, once on each parser. Without pyarrow
     # here the fast half silently skips, which is how a broken fast path ships.
     pyarrow

pyranges1/_ruranges.py

@@ -4,7 +4,7 @@
 from importlib.metadata import PackageNotFoundError, version
 from typing import Any
 
-MIN_RURANGES_VERSION = "0.1.5"
+MIN_RURANGES_VERSION = "0.2.9"
 
 
 def _version_tuple(v: str) -> tuple[int, ...]:

pyranges1/core/pyranges_main.py

@@ -61,7 +61,7 @@
 from pyranges1.core.tostring import tostring
 from pyranges1.methods.complement import _complement
 from pyranges1.methods.interval_metrics import compute_interval_metrics
-from pyranges1.methods.map_to_global import _map_to_global_pandas
+from pyranges1.methods.map_to_global import _map_to_global_ruranges
 from pyranges1.methods.map_to_local import _map_to_local
 from pyranges1.methods.sort import resolve_sort_keys, sort_order
 from pyranges1.range_frame.range_frame import RangeFrame
@@ -2007,7 +2007,7 @@ def map_to_global(
         else:
             local = self
 
-        return _map_to_global_pandas(
+        return _map_to_global_ruranges(
             local_gr=local, global_gr=gr, global_on=global_on, local_on=local_on, keep_id=keep_id, keep_loc=keep_loc
         )
 

pyranges1/methods/map_to_global.py

@@ -17,8 +17,6 @@
 
 ruranges = require_ruranges()
 
-cumsum_start = "_local_start"
-cumsum_end = "_local_end"
 out_local_start = "Start_local"
 out_local_end = "End_local"
 out_local_strand = "Strand_local"
@@ -89,7 +87,26 @@ def _map_to_global_pandas(
     return ensure_pyranges(z[cols_output])
 
 
-# buggy:
+def _transcript_offsets(
+    local_gr: "PyRanges", q_tx_code: np.ndarray, ex_tx_code: np.ndarray, ex_local_end: np.ndarray
+) -> tuple[np.ndarray, np.ndarray]:
+    """Interpret negative query offsets relative to their transcript length."""
+    q_start = local_gr[START_COL].to_numpy(copy=False)
+    q_end = local_gr[END_COL].to_numpy(copy=False)
+
+    # slice_ranges historically interprets negative offsets from the transcript
+    # end. Normalize those offsets before the core clips to exon coverage.
+    if (q_start < 0).any() or (q_end < 0).any():
+        n_groups = max(int(q_tx_code.max(initial=0)), int(ex_tx_code.max(initial=0))) + 1
+        transcript_lengths = np.zeros(n_groups, dtype=np.int64)
+        np.maximum.at(transcript_lengths, ex_tx_code, ex_local_end)
+        lengths = transcript_lengths[q_tx_code]
+        q_start = np.where(q_start < 0, lengths + q_start, q_start)
+        q_end = np.where(q_end < 0, lengths + q_end, q_end)
+
+    return q_start, q_end
+
+
 def _map_to_global_ruranges(
     local_gr: "PyRanges",
     global_gr: "PyRanges",
@@ -100,100 +117,87 @@ def _map_to_global_ruranges(
     keep_loc: bool = False,
     preserve_input_order: bool = True,
 ) -> "PyRanges":
-    """Lift intervals from local_gr to genomic coordinates using global_gr exon annotations.
-
-    Lift intervals in *local_gr* (coordinates relative to a transcript) onto absolute genomic coordinates using the exon annotation in *global_gr*.
-    """
+    """Project queries in Rust and rebuild their original index and metadata."""
     local_has_strand = local_gr.has_strand
     global_has_strand = global_gr.has_strand
-
-    global_cum = global_gr.group_cumsum(
+    reference = global_gr.get_with_loc_columns(global_on)
+
+    # The core models one spliced coordinate system per feature. Preserve the
+    # existing pandas behavior for ids reused on different chromosomes/strands.
+    locus_columns = [CHROM_COL] + ([STRAND_COL] if global_has_strand else [])
+    locus_columns = [column for column in locus_columns if column != global_on]
+    if locus_columns and (
+        reference.groupby(global_on, observed=True, dropna=False)[locus_columns]
+        .nunique(dropna=False)
+        .gt(1)
+        .any(axis=None)
+    ):
+        return _map_to_global_pandas(
+            local_gr,
+            global_gr,
+            global_on,
+            local_on=local_on,
+            keep_id=keep_id,
+            keep_loc=keep_loc,
+        )
+
+    # Keep temporary cumulative-coordinate columns off the calling frame.
+    local_start_column, local_end_column = "__map_start", "__map_end"
+    while local_start_column in reference or local_start_column == global_on:
+        local_start_column += "_"
+    while local_end_column in reference or local_end_column == global_on:
+        local_end_column += "_"
+    ex_df = reference.group_cumsum(
         group_by=global_on,
         use_strand="auto",
-        cumsum_start_column=cumsum_start,
-        cumsum_end_column=cumsum_end,
-        keep_order=True,  # False,
-    )
-    ex_df = global_cum.copy()
-
-    ex_df["__global_start__"] = ex_df[START_COL]
-    ex_df["__global_end__"] = ex_df[END_COL]
-    ex_df[START_COL] = ex_df[cumsum_start]
-    ex_df[END_COL] = ex_df[cumsum_end]
-
-    tmp_tx_id = "__txid__"
-    local_df = local_gr.copy()
-
-    ex_df[tmp_tx_id] = ex_df[global_on]
-    local_df[tmp_tx_id] = local_df[local_on]
-
-    q_tx_code, ex_tx_code = factorize_binary(local_df, ex_df, by=tmp_tx_id)
-
-    ex_chr_code, chr_table = pd.factorize(ex_df[CHROM_COL], sort=False)
-    ex_chr_code = ex_chr_code.astype(ex_tx_code.dtype, copy=False)
-
-    ex_local_start = ex_df[START_COL].to_numpy(copy=False)
-    ex_local_end = ex_df[END_COL].to_numpy(copy=False)
-    ex_genome_start = ex_df["__global_start__"].to_numpy(copy=False)
-    ex_genome_end = ex_df["__global_end__"].to_numpy(copy=False)
-    ex_fwd = (
-        (ex_df[STRAND_COL] == "+").to_numpy(np.bool_, copy=False)
-        if global_has_strand
-        else np.ones(len(ex_df), dtype=bool)
+        cumsum_start_column=local_start_column,
+        cumsum_end_column=local_end_column,
+        keep_order=True,
     )
-
-    q_start = local_df[START_COL].to_numpy(copy=False)
-    q_end = local_df[END_COL].to_numpy(copy=False)
-    q_fwd = (
-        (local_df[STRAND_COL] == "+").to_numpy(np.bool_, copy=False)
+    query_ids = pd.DataFrame({"feature": local_gr[local_on].array})
+    exon_ids = pd.DataFrame({"feature": ex_df[global_on].array})
+    q_tx_code, ex_tx_code = factorize_binary(query_ids, exon_ids, by="feature")
+    ex_local_start = ex_df[local_start_column].to_numpy(copy=False)
+    ex_local_end = ex_df[local_end_column].to_numpy(copy=False)
+    q_start, q_end = _transcript_offsets(local_gr, q_tx_code, ex_tx_code, ex_local_end)
+
+    keep_idx, out_start, out_end, out_strand, exon_idx = ruranges.numpy.map_to_global(
+        groups=q_tx_code,
+        starts=q_start,
+        ends=q_end,
+        strand=(local_gr[STRAND_COL] == "+").to_numpy(np.bool_, copy=False)
         if local_has_strand
-        else np.ones(len(local_df), dtype=bool)
-    )
-
-    keep_idx, out_start, out_end, out_strand_bool = ruranges.numpy.map_to_global(
-        groups=ex_tx_code,
-        starts=ex_local_start,
-        ends=ex_local_end,
-        strand=ex_fwd,
-        groups2=q_tx_code,
-        starts2=q_start,
-        ends2=q_end,
-        chr_code2=ex_chr_code,
-        genome_start2=ex_genome_start,
-        genome_end2=ex_genome_end,
-        strand2=q_fwd,
+        else np.ones(len(local_gr), dtype=bool),
+        groups2=ex_tx_code,
+        starts2=ex_local_start,
+        ends2=ex_local_end,
+        chr_code2=np.zeros(len(ex_df), dtype=np.uint32),
+        genome_start2=ex_df[START_COL].to_numpy(copy=False),
+        genome_end2=ex_df[END_COL].to_numpy(copy=False),
+        strand2=(ex_df[STRAND_COL] == "+").to_numpy(np.bool_, copy=False)
+        if global_has_strand
+        else np.ones(len(ex_df), dtype=bool),
         sort_output=preserve_input_order,
+        return_exon_indices=True,
     )
-
-    if len(keep_idx) == 0:
-        return ensure_pyranges(local_df.head(0))
-
-    max_code = int(ex_tx_code.max()) + 1
-    tx2chr_code = np.empty(max_code, dtype=ex_chr_code.dtype)
-    tx2chr_code[ex_tx_code] = ex_chr_code
-    out_chr_code = tx2chr_code[q_tx_code[keep_idx]]
-    out_chr_str = chr_table[out_chr_code]
-
+    # An empty merge reproduces pandas' key/category dtype rules without
+    # materializing the query-by-exon table used by the previous implementation.
+    schema = local_gr.head(0).merge(reference.head(0), left_on=local_on, right_on=global_on, suffixes=("_local", ""))
     mapped_df = local_gr.take(keep_idx).copy()
-
     if keep_id:
-        mapped_df[global_on] = mapped_df[CHROM_COL].copy()
+        mapped_df[global_on] = ex_df[global_on].iloc[exon_idx].astype(schema[global_on].dtype).array
+    elif local_on == global_on and local_on != CHROM_COL:
+        mapped_df[local_on] = mapped_df[local_on].astype(schema[local_on].dtype)
     if keep_loc:
-        mapped_df[out_local_start] = mapped_df[START_COL].copy()
-        mapped_df[out_local_end] = mapped_df[END_COL].copy()
+        mapped_df[out_local_start] = mapped_df[START_COL].array.copy()
+        mapped_df[out_local_end] = mapped_df[END_COL].array.copy()
         if local_has_strand and global_has_strand:
-            mapped_df[out_local_strand] = mapped_df[STRAND_COL].copy()
-
-    mapped_df[CHROM_COL] = out_chr_str
-    mapped_df[START_COL] = out_start
-    mapped_df[END_COL] = out_end
-
-    if local_has_strand or global_has_strand:
-        mapped_df[STRAND_COL] = np.where(out_strand_bool, "+", "-")
-
-    if tmp_tx_id in mapped_df.columns:
-        mapped_df = mapped_df.drop(columns=[tmp_tx_id])
-        if mapped_df is None:
-            raise ValueError
-
+            mapped_df[out_local_strand] = mapped_df[STRAND_COL].array.copy()
+    mapped_df[CHROM_COL] = ex_df[CHROM_COL].iloc[exon_idx].astype(schema[CHROM_COL].dtype).array
+    mapped_df[START_COL] = pd.Series(out_start).astype(global_gr[START_COL].dtype).array
+    mapped_df[END_COL] = pd.Series(out_end).astype(global_gr[END_COL].dtype).array
+    if local_has_strand and global_has_strand:
+        mapped_df[STRAND_COL] = np.where(out_strand, "+", "-")
+    elif global_has_strand:
+        mapped_df[STRAND_COL] = ex_df[STRAND_COL].iloc[exon_idx].array
     return ensure_pyranges(mapped_df)