diff --git a/CHANGELOG.txt b/CHANGELOG.txt index 124757b..1826f8b 100644 --- a/CHANGELOG.txt +++ b/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) diff --git a/README.md b/README.md index b54b45e..42c51af 100644 --- a/README.md +++ b/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) diff --git a/docs/how_to_map.rst b/docs/how_to_map.rst index fcd30e7..e1d6f90 100644 --- a/docs/how_to_map.rst +++ b/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. diff --git a/pyproject.toml b/pyproject.toml index d6c748f..4339153 100644 --- a/pyproject.toml +++ b/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 diff --git a/pyranges1/_ruranges.py b/pyranges1/_ruranges.py index bc645af..7d25578 100644 --- a/pyranges1/_ruranges.py +++ b/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, ...]: diff --git a/pyranges1/core/pyranges_main.py b/pyranges1/core/pyranges_main.py index bad2943..d34edec 100644 --- a/pyranges1/core/pyranges_main.py +++ b/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 ) diff --git a/pyranges1/methods/map_to_global.py b/pyranges1/methods/map_to_global.py index 46308cc..badc851 100644 --- a/pyranges1/methods/map_to_global.py +++ b/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) diff --git a/tests/test_parallelism.py b/tests/test_parallelism.py deleted file mode 100644 index 7675ba6..0000000 --- a/tests/test_parallelism.py +++ /dev/null @@ -1,39 +0,0 @@ -import pandas.testing -import pytest - -import pyranges1 as pr - - -def test_split_chunks() -> None: - df = pr.example_data.ensembl_gtf.get_with_loc_columns(["gene_id"]) - frame1, frame2 = pr.parallelism.split_df_into_chunks_without_splitting_groups(df, by=["gene_id"], nb_splits=2) - assert len({*frame1["gene_id"]}) == 1 - assert len({*frame2["gene_id"]}) == 1 - assert len(frame1) + len(frame2) == len(df) - - -@pytest.mark.slow -def test_split_chunks_many() -> None: - gr = pr.PyRanges( - {"Chromosome": ["chr1"] * 1000, "Start": range(1000), "End": range(1, 1001), "GeneId": [*range(100)] * 10} - ) - res = pr.parallelism.split_df_into_chunks_without_splitting_groups(gr, by=["GeneId"], nb_splits=37) - assert sum(len(r) for r in res) == len(gr) - seen_ids = set() - for i, r in enumerate(res): - new_ids = {*r["GeneId"].drop_duplicates()} - assert new_ids.isdisjoint(seen_ids), i - seen_ids.update(new_ids) - - -@pytest.mark.slow -def test_run_in_parallel() -> None: - df = pr.example_data.ensembl_gtf.get_with_loc_columns(["gene_id"]) - dfs = [gdf for _, gdf in df.groupby("gene_id")] - res = pr.parallelism.run_in_parallel( - function=lambda x: x.merge_overlaps(by="gene_id", count_col="Counts"), - dfs=dfs, - nb_cpu=2, - ).reset_index(drop=True) - expected_result = df.merge_overlaps(by="gene_id", count_col="Counts") - pandas.testing.assert_frame_equal(res, expected_result)