From 1c32281f5286624b340518176aca80d3c1e1ebc6 Mon Sep 17 00:00:00 2001 From: Sam Park Date: Fri, 3 Apr 2026 16:03:45 -0400 Subject: [PATCH 01/10] Add gtars genomicdist backend alongside R backend New files: - gtars_backend.py: GtarsStatBackend wrapping gtars genomicdist CLI - compress_distributions.py: fixed-size compression (histograms, KDE, dense arrays) - ref_utils.py: reference file resolution via refgenie + seqcol fallback Changes: - Register GtarsStatBackend in factory (replaces NotImplementedError) - Only create RServiceManager when backend is "r" - Skip R bedset plots when backend is "gtars" - Add --backend CLI override for run_all and run_stats - Add DEFAULT_PRECISION constant Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedboss.py | 5 +- bedboss/bedbuncher/bedbuncher.py | 19 +- bedboss/bedstat/backends/__init__.py | 5 +- bedboss/bedstat/backends/gtars_backend.py | 229 ++++++++++++++++++++ bedboss/bedstat/compress_distributions.py | 244 ++++++++++++++++++++++ bedboss/bedstat/ref_utils.py | 193 +++++++++++++++++ bedboss/cli.py | 12 ++ bedboss/const.py | 1 + 8 files changed, 700 insertions(+), 8 deletions(-) create mode 100644 bedboss/bedstat/backends/gtars_backend.py create mode 100644 bedboss/bedstat/compress_distributions.py create mode 100644 bedboss/bedstat/ref_utils.py diff --git a/bedboss/bedboss.py b/bedboss/bedboss.py index 0287e3b6..61af5132 100644 --- a/bedboss/bedboss.py +++ b/bedboss/bedboss.py @@ -420,7 +420,9 @@ def insert_pep( if rerun: skipper.reinitialize() - if not lite: + backend = bbagent.config.config.analysis.backend + + if not lite and backend == "r": r_service = RServiceManager() else: r_service = None @@ -504,6 +506,7 @@ def insert_pep( force_overwrite=force_overwrite, annotation=bedset_annotation, lite=lite, + backend=backend, ) else: _LOGGER.info( diff --git a/bedboss/bedbuncher/bedbuncher.py b/bedboss/bedbuncher/bedbuncher.py index 0a9c9ecd..1b2bf90d 100644 --- a/bedboss/bedbuncher/bedbuncher.py +++ b/bedboss/bedbuncher/bedbuncher.py @@ -102,6 +102,7 @@ def run_bedbuncher( no_fail: bool = False, force_overwrite: bool = False, lite: bool = False, + backend: str = "r", ) -> None: """ Add bedset to the database @@ -114,16 +115,17 @@ def run_bedbuncher( :param description: Bedset description :param annotation: Bedset annotation (author, source, summary, etc.) :param heavy: whether to use heavy processing (add all columns to the database). - if False -> R-script won't be executed, only basic statistics will be calculated + if False -> R-script won't be executed, only basic statistics will be calculated. + Ignored when backend is "gtars" (no R plots available). :param no_fail: whether to raise an error if bedset was not added to the database :param upload_pephub: whether to create a view in pephub :param upload_s3: whether to upload files to s3 :param force_overwrite: whether to overwrite the record in the database :param lite: whether to run the pipeline in lite mode - # TODO: force_overwrite is not working!!! Fix it! + :param backend: analysis backend ("r" or "gtars"). When "gtars", heavy is ignored. :return: """ - _LOGGER.info(f"Adding bedset { record_id} to the database") + _LOGGER.info(f"Adding bedset {record_id} to the database") if isinstance(bedbase_config, str): bbagent = BedBaseAgent(bedbase_config) @@ -140,7 +142,7 @@ def run_bedbuncher( "bedsets", ) - if heavy: + if heavy and backend != "gtars": _LOGGER.info("Heavy processing is True. Calculating plots...") plot_value = create_plots( bedset=bed_set, @@ -148,7 +150,12 @@ def run_bedbuncher( ) plots = BedSetPlots(region_commonality=FileModel(**plot_value)) else: - _LOGGER.info("Heavy processing is False. Plots won't be calculated") + if backend == "gtars" and heavy: + _LOGGER.info( + "Heavy processing ignored for gtars backend (no R plots available)" + ) + else: + _LOGGER.info("Heavy processing is False. Plots won't be calculated") plots = None bbagent.bedset.create( @@ -178,6 +185,7 @@ def run_bedbuncher_form_pep( upload_s3: bool = False, no_fail: bool = False, force_overwrite: bool = False, + backend: str = "r", ) -> str: """ Create bedset from pep and add it to the database @@ -224,6 +232,7 @@ def run_bedbuncher_form_pep( upload_s3=upload_s3, no_fail=no_fail, force_overwrite=force_overwrite, + backend=backend, ) return bedset_name diff --git a/bedboss/bedstat/backends/__init__.py b/bedboss/bedstat/backends/__init__.py index bcfff9c8..7578cb3c 100644 --- a/bedboss/bedstat/backends/__init__.py +++ b/bedboss/bedstat/backends/__init__.py @@ -1,8 +1,9 @@ from bedboss.bedstat.backends.base import StatBackend +from bedboss.bedstat.backends.gtars_backend import GtarsStatBackend from bedboss.bedstat.backends.r_backend import RStatBackend from bedboss.const import BACKEND_GTARS, BACKEND_R -__all__ = ["StatBackend", "RStatBackend", "create_backend"] +__all__ = ["StatBackend", "RStatBackend", "GtarsStatBackend", "create_backend"] def create_backend(name: str, **kwargs) -> StatBackend: @@ -15,6 +16,6 @@ def create_backend(name: str, **kwargs) -> StatBackend: if name == BACKEND_R: return RStatBackend(**kwargs) elif name == BACKEND_GTARS: - raise NotImplementedError("gtars backend not yet available. Install via PR 2a.") + return GtarsStatBackend(**kwargs) else: raise ValueError(f"Unknown analysis backend: {name!r}. Use 'r' or 'gtars'.") diff --git a/bedboss/bedstat/backends/gtars_backend.py b/bedboss/bedstat/backends/gtars_backend.py new file mode 100644 index 00000000..2968f989 --- /dev/null +++ b/bedboss/bedstat/backends/gtars_backend.py @@ -0,0 +1,229 @@ +import json +import logging +import os +import statistics +from pathlib import Path +from typing import Union + +import pypiper +from gtars.models import RegionSet + +from bedboss.bedstat.backends.base import StatBackend +from bedboss.bedstat.compress_distributions import ( + compress_distributions, + compress_to_kde, +) +from bedboss.bedstat.gc_content import calculate_gc_content +from bedboss.bedstat.ref_utils import ( + get_chrom_sizes_path, + get_gda_path, + get_osm_path_with_precompile, +) +from bedboss.const import BEDSTAT_OUTPUT, DEFAULT_PRECISION, OUTPUT_FOLDER_NAME +from bedboss.exceptions import BedBossException + +_LOGGER = logging.getLogger("bedboss") + +# Map gtars partition names to legacy DB column name prefixes +PARTITION_NAME_MAP = { + "promoterCore": "promotercore", + "promoterProx": "promoterprox", + "threeUTR": "threeutr", + "fiveUTR": "fiveutr", + "exon": "exon", + "intron": "intron", + "intergenic": "intergenic", +} + + +def round_floats(obj, precision: int = 4): + """Recursively round all floats in a nested dict/list structure.""" + if isinstance(obj, float): + return round(obj, precision) + elif isinstance(obj, dict): + return {k: round_floats(v, precision) for k, v in obj.items()} + elif isinstance(obj, list): + return [round_floats(v, precision) for v in obj] + return obj + + +class GtarsStatBackend(StatBackend): + """gtars genomicdist-based statistics backend.""" + + def __init__( + self, + region_dist_bins: int = 250, + promoter_upstream: int = 200, + promoter_downstream: int = 2000, + precision: int = DEFAULT_PRECISION, + **kwargs, + ): + self._region_dist_bins = region_dist_bins + self._promoter_upstream = promoter_upstream + self._promoter_downstream = promoter_downstream + self._precision = precision + + def compute( + self, + bedfile: str, + genome: str, + outfolder: str, + bed_digest: str = None, + ensdb: str = None, + open_signal_matrix: str = None, + just_db_commit: bool = False, + rfg_config: Union[str, Path] = None, + pm: pypiper.PipelineManager = None, + ) -> dict: + # Auto-fetch GDA/GTF annotation via refgenie if not provided + if not ensdb: + try: + ensdb = get_gda_path(genome, rfg_config=rfg_config) + except Exception: + _LOGGER.warning( + f"Could not fetch annotation for {genome}. " + "Partition and TSS analysis will be skipped." + ) + + # Pre-compile open signal matrix to .bin if available + if open_signal_matrix and os.path.exists(open_signal_matrix): + open_signal_matrix = get_osm_path_with_precompile(open_signal_matrix) + + # Auto-fetch chrom.sizes via refgenie if not provided + chrom_sizes = None + try: + chrom_sizes = get_chrom_sizes_path(genome, rfg_config=rfg_config) + except Exception: + _LOGGER.warning( + f"Could not fetch chrom.sizes for {genome}. " + "Region distribution will not be normalized." + ) + + outfolder_stats = os.path.join(outfolder, OUTPUT_FOLDER_NAME, BEDSTAT_OUTPUT) + os.makedirs(outfolder_stats, exist_ok=True) + + # Used to stop pipeline if bedstat is used independently + stop_pipeline = not pm + + bed_object = RegionSet(bedfile) + if not bed_digest: + bed_digest = bed_object.identifier + + outfolder_stats_results = os.path.abspath( + os.path.join(outfolder_stats, bed_digest) + ) + os.makedirs(outfolder_stats_results, exist_ok=True) + + json_file_path = os.path.abspath( + os.path.join(outfolder_stats_results, bed_digest + ".json") + ) + + if not just_db_commit: + if not pm: + pm_out_path = os.path.abspath( + os.path.join(outfolder_stats, "pypiper", bed_digest) + ) + os.makedirs(pm_out_path, exist_ok=True) + pm = pypiper.PipelineManager( + name="bedstat-pipeline", + outfolder=pm_out_path, + pipestat_sample_name=bed_digest, + ) + + # Build gtars genomicdist command + cmd_parts = [ + "gtars", + "genomicdist", + "--bed", + bedfile, + "--output", + json_file_path, + ] + if ensdb: + cmd_parts.extend(["--gtf", ensdb]) + if chrom_sizes: + cmd_parts.extend(["--chrom-sizes", chrom_sizes]) + if open_signal_matrix: + cmd_parts.extend(["--signal-matrix", open_signal_matrix]) + cmd_parts.extend( + [ + "--bins", + str(self._region_dist_bins), + "--promoter-upstream", + str(self._promoter_upstream), + "--promoter-downstream", + str(self._promoter_downstream), + "--compact", + ] + ) + + command = " ".join(cmd_parts) + try: + _LOGGER.info(f"Running gtars genomicdist: {command}") + pm.run(cmd=command, target=json_file_path) + except Exception as e: + _LOGGER.error(f"gtars genomicdist failed: {e}") + raise BedBossException(f"gtars genomicdist failed: {e}") + + # Read gtars JSON output + gtars_output = {} + if os.path.exists(json_file_path): + with open(json_file_path, "r", encoding="utf-8") as f: + gtars_output = json.load(f) + + # Extract scalars to flat dict keys + data = {} + scalars = gtars_output.get("scalars", {}) + data["number_of_regions"] = scalars.get("number_of_regions") + data["mean_region_width"] = scalars.get("mean_region_width") + data["median_tss_dist"] = scalars.get("median_tss_dist") + + # Populate legacy partition flat columns + partitions = gtars_output.get("partitions") + if partitions: + total = partitions.get("total", 0) + for name, count in partitions.get("counts", []): + db_name = PARTITION_NAME_MAP.get(name) + if db_name and total > 0: + data[f"{db_name}_frequency"] = count + data[f"{db_name}_percentage"] = round(count / total * 100, 4) + + # GC content: compute via Python bindings (requires refgenie FASTA) + try: + gc_contents = calculate_gc_content( + bedfile=bed_object, genome=genome, rfg_config=rfg_config + ) + except BaseException as e: + _LOGGER.warning( + f"GC content calculation skipped for {genome}: {e}. " + "Ensure refgenie is configured with a FASTA asset." + ) + gc_contents = None + + if gc_contents: + gc_mean = round(statistics.mean(gc_contents), 4) + data["gc_content"] = gc_mean + + # Compress per-region GC values to 512-pt KDE, inject into distributions + gc_kde = compress_to_kde(gc_contents, n_points=512, log_transform=False) + if gc_kde: + gc_kde["mean"] = gc_mean + if "distributions" not in gtars_output: + gtars_output["distributions"] = {} + gtars_output["distributions"]["gc_content"] = gc_kde + else: + data["gc_content"] = None + + # Compress distributions for DB storage + compress_distributions(gtars_output) + + # Store entire augmented gtars JSON as distributions blob + data["distributions"] = gtars_output + + if self._precision is not None: + data = round_floats(data, self._precision) + + if stop_pipeline and pm: + pm.stop_pipeline() + + return data diff --git a/bedboss/bedstat/compress_distributions.py b/bedboss/bedstat/compress_distributions.py new file mode 100644 index 00000000..b014155c --- /dev/null +++ b/bedboss/bedstat/compress_distributions.py @@ -0,0 +1,244 @@ +"""Compress raw gtars genomicdist distributions for DB storage. + +Produces fixed-size representations (~18KB total) regardless of region count. +Output formats match what bedbase-ui expects for client-side rendering. +""" + +from typing import List, Optional + +import numpy as np + + +def compress_to_histogram( + values: List[float], n_bins: int = 50, trim_percentile: float = 99.0 +) -> Optional[dict]: + """Compress a list of values to a fixed-bin histogram. + + Uses percentile trimming to remove outliers, matching the UI's + quantileTrimmedHistogram() behavior. + + Returns: {"x_min", "x_max", "bins", "counts", "total"} + """ + if not values: + return None + + arr = np.asarray(values, dtype=np.float64) + total = len(arr) + cutoff = np.percentile(arr, trim_percentile) + trimmed = arr[arr <= cutoff] + if len(trimmed) == 0: + return None + + x_min = float(trimmed.min()) + x_max = float(trimmed.max()) + if x_min == x_max: + return { + "x_min": x_min, + "x_max": x_max, + "bins": 1, + "counts": [int(len(trimmed))], + "total": total, + } + + counts, _ = np.histogram(trimmed, bins=n_bins, range=(x_min, x_max)) + overflow = total - len(trimmed) + return { + "x_min": x_min, + "x_max": x_max, + "bins": n_bins, + "counts": counts.tolist(), + "total": total, + "overflow": overflow, + } + + +def compress_to_kde( + values: List[float], + n_points: int = 512, + log_transform: bool = False, + trim_percentile: float = 99.0, + max_samples: int = 5000, +) -> Optional[dict]: + """Compress a list of values to a Gaussian KDE curve. + + Replicates the UI's exact math (genomicdist-plots.ts): + 1. Optional log10 transform (for neighbor_distances) + 2. Percentile trim + downsample + 3. Silverman bandwidth + 4. Evaluate Gaussian kernel over evenly-spaced points + + Returns: {"x_min", "x_max", "n", "densities"} + """ + if not values or len(values) < 2: + return None + + arr = np.asarray(values, dtype=np.float64) + + # Optional log10 transform (filter non-positive first) + if log_transform: + arr = arr[arr > 0] + if len(arr) < 2: + return None + arr = np.log10(arr) + + if len(arr) < 2: + return None + + # Percentile trim + cutoff = np.percentile(arr, trim_percentile) + trimmed = np.sort(arr[arr <= cutoff]) + if len(trimmed) < 2: + return None + + # Compute bandwidth stats from FULL trimmed data (matching UI behavior) + full_n = len(trimmed) + sd = float(np.std(trimmed, ddof=0)) + if sd == 0: + sd = 1e-10 + + # IQR + q1 = float(np.percentile(trimmed, 25)) + q3 = float(np.percentile(trimmed, 75)) + iqr = q3 - q1 + + # Silverman bandwidth (using full trimmed count, not downsampled) + h = 0.9 * min(sd, iqr / 1.34 if iqr > 0 else sd) * (full_n**-0.2) + + # Downsample AFTER bandwidth computation (only affects KDE evaluation speed) + if len(trimmed) > max_samples: + indices = np.linspace(0, len(trimmed) - 1, max_samples, dtype=int) + trimmed = trimmed[indices] + + n = len(trimmed) + if h <= 0: + h = sd * (n**-0.2) + if h <= 0: + return None + + # Evaluation range: min - 3h to max + 3h + x_min = float(trimmed[0]) - 3 * h + x_max = float(trimmed[-1]) + 3 * h + + # Vectorized Gaussian kernel evaluation + xs = np.linspace(x_min, x_max, n_points) + # Shape: (n_points, n_samples) — broadcast subtract + u = (xs[:, np.newaxis] - trimmed[np.newaxis, :]) / h + densities = np.sum(np.exp(-0.5 * u * u), axis=1) / (h * np.sqrt(2 * np.pi) * n) + + return { + "x_min": round(x_min, 6), + "x_max": round(x_max, 6), + "n": n_points, + "densities": np.round(densities, 8).tolist(), + } + + +def compress_tss_histogram( + values: List[float], n_bins: int = 100, max_distance: float = 100_000.0 +) -> Optional[dict]: + """Compress signed TSS distances to a fixed-range symmetric histogram. + + Expects signed distances (negative = upstream, positive = downstream) + from gtars calc_feature_distances. Bins into [-max_distance, +max_distance] + to match the UI's local TSS distance plot. + + Returns: {"x_min", "x_max", "bins", "counts", "total"} + """ + if not values: + return None + + arr = np.asarray(values, dtype=np.float64) + total = len(arr) + + # Clamp to symmetric range + clamped = arr[(arr >= -max_distance) & (arr <= max_distance)] + if len(clamped) == 0: + return None + + counts, _ = np.histogram(clamped, bins=n_bins, range=(-max_distance, max_distance)) + return { + "x_min": -max_distance, + "x_max": max_distance, + "bins": n_bins, + "counts": counts.tolist(), + "total": total, + } + + +def compress_region_distribution(raw: dict) -> Optional[dict]: + """Compress per-chromosome region distribution to dense count arrays. + + Input: gtars format {"chr1": [{"start": ..., "end": ..., "rid": ...}, ...], ...} + Output: {"chr1": [count_at_rid_0, count_at_rid_1, ...], ...} + + The array index is the rid (bin index) used by the UI's faceted chart. + """ + if not raw: + return None + + result = {} + for chrom, regions in raw.items(): + if not regions: + result[chrom] = [] + continue + rids = np.array( + [r.get("rid", 0) if isinstance(r, dict) else 0 for r in regions], + dtype=np.int32, + ) + counts_arr = np.array( + [r.get("n", 1) if isinstance(r, dict) else 1 for r in regions], + dtype=np.int32, + ) + bins = np.zeros(rids.max() + 1, dtype=np.int64) + np.add.at(bins, rids, counts_arr) + result[chrom] = bins.tolist() + + return result + + +def compress_distributions(gtars_output: dict) -> dict: + """Compress all distributions in a gtars genomicdist output. + + Modifies gtars_output["distributions"] in place, replacing raw arrays + with compressed formats. + + Returns the modified gtars_output. + """ + dists = gtars_output.get("distributions", {}) + + # Widths: histogram (50 bins) + if "widths" in dists and isinstance(dists["widths"], list): + dists["widths"] = compress_to_histogram(dists["widths"], n_bins=50) + + # TSS distances: fixed-range histogram (100 bins, 0–100kb) + if "tss_distances" in dists and isinstance(dists["tss_distances"], list): + dists["tss_distances"] = compress_tss_histogram( + dists["tss_distances"], n_bins=100 + ) + + # Neighbor distances: KDE with log10 transform + if "neighbor_distances" in dists and isinstance(dists["neighbor_distances"], list): + dists["neighbor_distances"] = compress_to_kde( + dists["neighbor_distances"], n_points=512, log_transform=True + ) + + # Drop nearest_neighbors (redundant with neighbor_distances) + dists.pop("nearest_neighbors", None) + + # Region distribution: dense count arrays + # gtars outputs a flat list of {chr, start, end, n, rid}; group by chr first + if "region_distribution" in dists: + rd = dists["region_distribution"] + if isinstance(rd, list): + grouped = {} + for entry in rd: + chrom = entry.get("chr", "unknown") + grouped.setdefault(chrom, []).append(entry) + rd = grouped + if isinstance(rd, dict): + dists["region_distribution"] = compress_region_distribution(rd) + + # chromosome_stats: unchanged (already compact) + + gtars_output["distributions"] = dists + return gtars_output diff --git a/bedboss/bedstat/ref_utils.py b/bedboss/bedstat/ref_utils.py new file mode 100644 index 00000000..355596fe --- /dev/null +++ b/bedboss/bedstat/ref_utils.py @@ -0,0 +1,193 @@ +"""Reference file resolution for gtars genomicdist. + +Handles auto-fetching and pre-compilation of GTF annotations, chrom.sizes, +and open signal matrices via refgenie with seqcol API fallback. +""" + +import logging +import os +import subprocess +from typing import Union + +from bedboss.const import HOME_PATH + +_LOGGER = logging.getLogger("bedboss") + + +def _get_chrom_sizes_seqcol(genome: str) -> Union[str, None]: + """Fallback: fetch chrom.sizes via seqcol API when refgenie doesn't have it. + + Resolves genome name to a seqcol digest via the refgenie /v4/genomes + endpoint, then fetches chromosome names + lengths directly. + """ + import requests + from refget.clients import SequenceCollectionClient + + refgenie_api = "https://api.refgenie.org" + cache_dir = os.path.join(HOME_PATH, "chrom_sizes") + + chrom_sizes_path = os.path.join(cache_dir, f"{genome}.chrom.sizes") + if os.path.exists(chrom_sizes_path): + _LOGGER.info(f"Chrom sizes (seqcol cache): {chrom_sizes_path}") + return chrom_sizes_path + + # Resolve genome name -> seqcol digest + try: + resp = requests.get( + f"{refgenie_api}/v4/genomes", + params={"limit": 1000}, + timeout=30, + ) + resp.raise_for_status() + except Exception as e: + _LOGGER.warning(f"seqcol fallback: failed to query genome list: {e}") + return None + + genome_lower = genome.lower() + digest = None + fallback_digest = None + for entry in resp.json().get("items", []): + for alias in entry.get("aliases", []): + alias_lower = alias.lower() + if alias_lower == f"{genome_lower}-refgenie": + digest = entry["digest"] + break + if not fallback_digest and alias_lower.startswith(genome_lower): + fallback_digest = entry["digest"] + if digest: + break + digest = digest or fallback_digest + + if not digest: + _LOGGER.warning(f"seqcol fallback: no digest found for '{genome}'") + return None + + try: + client = SequenceCollectionClient(urls=[f"{refgenie_api}/seqcol"]) + os.makedirs(cache_dir, exist_ok=True) + client.write_chrom_sizes(digest, chrom_sizes_path) + _LOGGER.info(f"Chrom sizes (seqcol): {chrom_sizes_path}") + return chrom_sizes_path + except Exception as e: + _LOGGER.warning(f"seqcol fallback: failed to fetch chrom.sizes: {e}") + return None + + +def get_chrom_sizes_path(genome: str, rfg_config=None) -> Union[str, None]: + """Get a chrom.sizes file for a genome. + + Tries refgenie first (rgc.seek/pull), falls back to the seqcol API. + """ + from refgenconf import RefgenconfError + from yacman.exceptions import UndefinedAliasError + + from bedboss.bedmaker.utils import get_rgc + + rgc = get_rgc(rfg_config=rfg_config) + try: + return rgc.seek( + genome_name=genome, + asset_name="fasta", + tag_name="default", + seek_key="chrom_sizes", + ) + except (UndefinedAliasError, RefgenconfError): + _LOGGER.info(f"chrom.sizes not local for {genome}, pulling from refgenie") + try: + rgc.pull(genome=genome, asset="fasta", tag="default") + return rgc.seek( + genome_name=genome, + asset_name="fasta", + tag_name="default", + seek_key="chrom_sizes", + ) + except Exception: + _LOGGER.info(f"refgenie pull failed for {genome}, trying seqcol API") + + return _get_chrom_sizes_seqcol(genome) + + +def get_gda_path(genome: str, rfg_config=None) -> Union[str, None]: + """Get a GDA (GenomicDist Annotation) binary for a genome. + + Pulls the Ensembl GTF from refgenie and pre-compiles it to a GDA .bin + using ``gtars prep``. Falls back to the raw .gtf.gz if compilation fails. + """ + from refgenconf import RefgenconfError + from yacman.exceptions import UndefinedAliasError + + from bedboss.bedmaker.utils import get_rgc + + _LOGGER.info(f"Getting GDA annotation for genome: {genome}") + rgc = get_rgc(rfg_config=rfg_config) + + try: + gtf_path = rgc.seek( + genome_name=genome, + asset_name="ensembl_gtf", + tag_name="default", + seek_key="ensembl_gtf", + ) + except (UndefinedAliasError, RefgenconfError): + _LOGGER.info(f"ensembl_gtf not local for {genome}, pulling from refgenie") + try: + rgc.pull(genome=genome, asset="ensembl_gtf", tag="default") + gtf_path = rgc.seek( + genome_name=genome, + asset_name="ensembl_gtf", + tag_name="default", + seek_key="ensembl_gtf", + ) + except Exception as e: + _LOGGER.warning(f"Could not fetch GTF for {genome}: {e}") + return None + + gda_bin_path = gtf_path + ".gda.bin" + + # Return pre-compiled GDA binary if it already exists + if os.path.exists(gda_bin_path): + _LOGGER.info(f"GDA annotation (pre-compiled): {gda_bin_path}") + return gda_bin_path + + # Pre-compile to GDA .bin + _LOGGER.info(f"Pre-compiling GDA: {gtf_path}") + result = subprocess.run( + ["gtars", "prep", "--gtf", gtf_path, "-o", gda_bin_path], + capture_output=True, + text=True, + ) + if result.returncode == 0 and os.path.exists(gda_bin_path): + _LOGGER.info(f"GDA annotation (pre-compiled): {gda_bin_path}") + return gda_bin_path + else: + _LOGGER.warning(f"gtars prep (GDA) failed, using raw GTF: {result.stderr}") + + return gtf_path + + +def get_osm_path_with_precompile( + osm_path: str, +) -> str: + """Pre-compile an open signal matrix to .bin if not already done. + + Returns the .bin path if compilation succeeds, otherwise the original path. + """ + osm_bin_path = osm_path + ".bin" + + if os.path.exists(osm_bin_path): + _LOGGER.info(f"Open Signal Matrix (pre-compiled): {osm_bin_path}") + return osm_bin_path + + _LOGGER.info(f"Pre-compiling signal matrix: {osm_path}") + result = subprocess.run( + ["gtars", "prep", "--signal-matrix", osm_path, "-o", osm_bin_path], + capture_output=True, + text=True, + ) + if result.returncode == 0 and os.path.exists(osm_bin_path): + _LOGGER.info(f"Open Signal Matrix (pre-compiled): {osm_bin_path}") + return osm_bin_path + else: + _LOGGER.warning(f"gtars prep failed, using raw file: {result.stderr}") + + return osm_path diff --git a/bedboss/cli.py b/bedboss/cli.py index e899ca41..1816bf76 100644 --- a/bedboss/cli.py +++ b/bedboss/cli.py @@ -99,6 +99,10 @@ def run_all( upload_qdrant: bool = typer.Option(False, help="Upload to Qdrant"), upload_s3: bool = typer.Option(False, help="Upload to S3"), upload_pephub: bool = typer.Option(False, help="Upload to PEPHub"), + backend: str = typer.Option( + None, + help="Override analysis backend ('r' or 'gtars'). If not set, uses config file value.", + ), # Universes universe: bool = typer.Option(False, help="Create a universe"), universe_method: str = typer.Option( @@ -122,6 +126,9 @@ def run_all( agent = BedBaseAgent(bedbase_config) + if backend: + agent.config.config.analysis.backend = backend + run_all_bedboss( input_file=input_file, input_type=input_type, @@ -394,6 +401,10 @@ def run_stats( None, help="Path to the open signal matrix file" ), just_db_commit: bool = typer.Option(False, help="Just commit to the database?"), + backend: str = typer.Option( + "r", + help="Analysis backend ('r' or 'gtars'). Default: 'r'.", + ), # PipelineManager multi: bool = typer.Option(False, help="Run multiple samples"), recover: bool = typer.Option(True, help="Recover from previous run"), @@ -408,6 +419,7 @@ def run_stats( ensdb=ensdb, open_signal_matrix=open_signal_matrix, just_db_commit=just_db_commit, + backend=backend, pm=create_pm(outfolder=outfolder, multi=multi, recover=recover, dirty=dirty), ) diff --git a/bedboss/const.py b/bedboss/const.py index 4cc6b1a0..9f7c2ecd 100644 --- a/bedboss/const.py +++ b/bedboss/const.py @@ -30,6 +30,7 @@ # bedstat BACKEND_R = "r" BACKEND_GTARS = "gtars" +DEFAULT_PRECISION = 3 # bedbuncher DEFAULT_BEDBASE_CACHE_PATH = "./bedabse_cache" From 9e7f90dc46e1b1542d7757978042a5e682f2c94f Mon Sep 17 00:00:00 2001 From: Sam Park Date: Fri, 3 Apr 2026 18:57:12 -0400 Subject: [PATCH 02/10] Prefer seqcol API over FASTA pull for chrom.sizes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Reorder chrom.sizes resolution: local refgenie cache → seqcol API (~50KB) → refgenie FASTA pull (~3GB). Avoids downloading a full genome FASTA just to get chromosome sizes on cold start. Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedstat/ref_utils.py | 37 +++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/bedboss/bedstat/ref_utils.py b/bedboss/bedstat/ref_utils.py index 355596fe..7c47cf71 100644 --- a/bedboss/bedstat/ref_utils.py +++ b/bedboss/bedstat/ref_utils.py @@ -76,13 +76,15 @@ def _get_chrom_sizes_seqcol(genome: str) -> Union[str, None]: def get_chrom_sizes_path(genome: str, rfg_config=None) -> Union[str, None]: """Get a chrom.sizes file for a genome. - Tries refgenie first (rgc.seek/pull), falls back to the seqcol API. + Tries refgenie local seek first, then seqcol API (lightweight ~50KB), + then refgenie pull as last resort (pulls full FASTA asset ~3GB). """ from refgenconf import RefgenconfError from yacman.exceptions import UndefinedAliasError from bedboss.bedmaker.utils import get_rgc + # 1. Check local refgenie cache (instant) rgc = get_rgc(rfg_config=rfg_config) try: return rgc.seek( @@ -92,19 +94,28 @@ def get_chrom_sizes_path(genome: str, rfg_config=None) -> Union[str, None]: seek_key="chrom_sizes", ) except (UndefinedAliasError, RefgenconfError): - _LOGGER.info(f"chrom.sizes not local for {genome}, pulling from refgenie") - try: - rgc.pull(genome=genome, asset="fasta", tag="default") - return rgc.seek( - genome_name=genome, - asset_name="fasta", - tag_name="default", - seek_key="chrom_sizes", - ) - except Exception: - _LOGGER.info(f"refgenie pull failed for {genome}, trying seqcol API") + pass + + # 2. Try seqcol API (lightweight — fetches only chrom.sizes, ~50KB) + _LOGGER.info(f"chrom.sizes not local for {genome}, trying seqcol API") + result = _get_chrom_sizes_seqcol(genome) + if result: + return result + + # 3. Last resort: pull full FASTA asset from refgenie (~3GB) + _LOGGER.info(f"seqcol failed for {genome}, pulling FASTA asset from refgenie") + try: + rgc.pull(genome=genome, asset="fasta", tag="default") + return rgc.seek( + genome_name=genome, + asset_name="fasta", + tag_name="default", + seek_key="chrom_sizes", + ) + except Exception: + _LOGGER.warning(f"Could not fetch chrom.sizes for {genome}") - return _get_chrom_sizes_seqcol(genome) + return None def get_gda_path(genome: str, rfg_config=None) -> Union[str, None]: From 30f08bc6f3ec769a9cd0481f6dd5c4b7c2881431 Mon Sep 17 00:00:00 2001 From: Sam Park Date: Fri, 3 Apr 2026 23:38:03 -0400 Subject: [PATCH 03/10] Fix partition percentage scale to match R (0-1, not 0-100) R stores partition percentages as fractions (0.0615), not percent (6.082). Remove the * 100 multiplier so gtars output matches the existing database convention. Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedstat/backends/gtars_backend.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bedboss/bedstat/backends/gtars_backend.py b/bedboss/bedstat/backends/gtars_backend.py index 2968f989..458d1fa1 100644 --- a/bedboss/bedstat/backends/gtars_backend.py +++ b/bedboss/bedstat/backends/gtars_backend.py @@ -186,7 +186,7 @@ def compute( db_name = PARTITION_NAME_MAP.get(name) if db_name and total > 0: data[f"{db_name}_frequency"] = count - data[f"{db_name}_percentage"] = round(count / total * 100, 4) + data[f"{db_name}_percentage"] = round(count / total, 4) # GC content: compute via Python bindings (requires refgenie FASTA) try: From a577e87b1a1fe803abf2ec55235759555eb38087 Mon Sep 17 00:00:00 2001 From: Sam Park Date: Sun, 5 Apr 2026 15:14:22 -0400 Subject: [PATCH 04/10] Compute median_neighbor_distance scalar in gtars backend MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Derives median absolute neighbor distance from the gtars CLI's raw neighbor_distances list, stores it in the data dict so it gets written to the new BedStats.median_neighbor_distance column (added in bbconf PR #113). Computed BEFORE compress_distributions replaces the flat list with a KDE. The per-file KDE is still stored in the distributions JSONB blob for single-file views; only the scalar is used for bedset aggregation (mean ± sd across files). Replaces what used to be an aggregated neighbor_distances KDE at the bedset level — per earlier discussion, per-file KDE variance is low within bedsets (assay type dominates), and a scalar median gives enough signal at collection level. Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedstat/backends/gtars_backend.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/bedboss/bedstat/backends/gtars_backend.py b/bedboss/bedstat/backends/gtars_backend.py index 458d1fa1..65552db0 100644 --- a/bedboss/bedstat/backends/gtars_backend.py +++ b/bedboss/bedstat/backends/gtars_backend.py @@ -178,6 +178,24 @@ def compute( data["mean_region_width"] = scalars.get("mean_region_width") data["median_tss_dist"] = scalars.get("median_tss_dist") + # Derive median_neighbor_distance from the raw neighbor_distances list. + # (The gtars CLI output includes the full list; we reduce it to one + # scalar here rather than aggregating the full distribution at bedset + # level — see bbconf aggregation decisions.) + neighbor_distances = gtars_output.get("distributions", {}).get( + "neighbor_distances" + ) + if neighbor_distances: + abs_vals = [abs(d) for d in neighbor_distances if d is not None] + if abs_vals: + data["median_neighbor_distance"] = round( + statistics.median(abs_vals), 4 + ) + else: + data["median_neighbor_distance"] = None + else: + data["median_neighbor_distance"] = None + # Populate legacy partition flat columns partitions = gtars_output.get("partitions") if partitions: From a1ef7e45f0a4acc1edcaf7f668f3b6d9eed2d266 Mon Sep 17 00:00:00 2001 From: Sam Park Date: Sun, 5 Apr 2026 15:22:39 -0400 Subject: [PATCH 05/10] Apply black formatting to gtars_backend.py Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedstat/backends/gtars_backend.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/bedboss/bedstat/backends/gtars_backend.py b/bedboss/bedstat/backends/gtars_backend.py index 65552db0..0914c790 100644 --- a/bedboss/bedstat/backends/gtars_backend.py +++ b/bedboss/bedstat/backends/gtars_backend.py @@ -188,9 +188,7 @@ def compute( if neighbor_distances: abs_vals = [abs(d) for d in neighbor_distances if d is not None] if abs_vals: - data["median_neighbor_distance"] = round( - statistics.median(abs_vals), 4 - ) + data["median_neighbor_distance"] = round(statistics.median(abs_vals), 4) else: data["median_neighbor_distance"] = None else: From 2c09caed419394c1e8e9d1285b4f18cf8c89ce4b Mon Sep 17 00:00:00 2001 From: Sam Park Date: Sun, 5 Apr 2026 15:43:17 -0400 Subject: [PATCH 06/10] Add GtarsPyStatBackend for side-by-side performance comparison New 'gtars-py' backend that uses gtars Python bindings directly, no subprocess. Exists alongside the existing 'gtars' CLI-subprocess backend so we can benchmark both on real data and keep only the better performer before merge. Architecture: - GtarsPyStatBackend in backends/gtars_py_backend.py - GenomeRefs dataclass caches reference data per genome on the backend instance (FASTA, chrom_sizes, GeneModel, PartitionList, TssIndex, SignalMatrix). Loaded lazily on first compute() call for a given genome, reused for all subsequent files. - Same compute() signature as GtarsStatBackend - Same output dict shape (compress_distributions + DB insertion unchanged) - Same pypiper integration for status tracking (pm.timestamp markers between Python steps, no subprocess wrapping) - median_neighbor_distance scalar derived same way Factory dispatch: create_backend("gtars-py") returns the new backend. Extended bbconf's analysis.backend Literal to accept "gtars-py" (bbconf commit on modular-backend-logic branch). Output shape normalization: the gtars Python binding returns calc_partitions as {partition: [names], count: [counts], total: n} but the CLI JSON schema has {counts: [[name, count], ...], total: n}. _normalize_partitions() converts between them so downstream code sees the same shape regardless of backend. Parity verified on a 126K-region hg38 ENCODE BED file: - Scalars identical (number_of_regions, mean_region_width, median_tss_dist, median_neighbor_distance, gc_content) - All 14 partition fields (frequency+percentage for 7 categories) match - widths, tss, neighbor_distances, region_distribution bins byte-match after compress_distributions Performance (5 ENCODE files, hg38): - gtars (CLI): total 11.15s, median per-file 2.55s - gtars-py: total 5.71s, median per-file 1.15s - Speedup: ~2x batch, ~2.2x per-file median - Biggest wins: FASTA loaded once (not per file), no subprocess startup, no JSON round-trip to disk Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedbuncher/bedbuncher.py | 9 +- bedboss/bedstat/backends/__init__.py | 19 +- bedboss/bedstat/backends/gtars_py_backend.py | 465 +++++++++++++++++++ bedboss/cli.py | 2 +- bedboss/const.py | 1 + 5 files changed, 487 insertions(+), 9 deletions(-) create mode 100644 bedboss/bedstat/backends/gtars_py_backend.py diff --git a/bedboss/bedbuncher/bedbuncher.py b/bedboss/bedbuncher/bedbuncher.py index 1b2bf90d..0a1c700f 100644 --- a/bedboss/bedbuncher/bedbuncher.py +++ b/bedboss/bedbuncher/bedbuncher.py @@ -122,7 +122,7 @@ def run_bedbuncher( :param upload_s3: whether to upload files to s3 :param force_overwrite: whether to overwrite the record in the database :param lite: whether to run the pipeline in lite mode - :param backend: analysis backend ("r" or "gtars"). When "gtars", heavy is ignored. + :param backend: analysis backend ("r", "gtars", or "gtars-py"). For gtars variants, heavy is ignored. :return: """ _LOGGER.info(f"Adding bedset {record_id} to the database") @@ -142,7 +142,8 @@ def run_bedbuncher( "bedsets", ) - if heavy and backend != "gtars": + gtars_like = backend in ("gtars", "gtars-py") + if heavy and not gtars_like: _LOGGER.info("Heavy processing is True. Calculating plots...") plot_value = create_plots( bedset=bed_set, @@ -150,9 +151,9 @@ def run_bedbuncher( ) plots = BedSetPlots(region_commonality=FileModel(**plot_value)) else: - if backend == "gtars" and heavy: + if gtars_like and heavy: _LOGGER.info( - "Heavy processing ignored for gtars backend (no R plots available)" + f"Heavy processing ignored for {backend} backend (no R plots available)" ) else: _LOGGER.info("Heavy processing is False. Plots won't be calculated") diff --git a/bedboss/bedstat/backends/__init__.py b/bedboss/bedstat/backends/__init__.py index 7578cb3c..c1c5028c 100644 --- a/bedboss/bedstat/backends/__init__.py +++ b/bedboss/bedstat/backends/__init__.py @@ -1,15 +1,22 @@ from bedboss.bedstat.backends.base import StatBackend from bedboss.bedstat.backends.gtars_backend import GtarsStatBackend +from bedboss.bedstat.backends.gtars_py_backend import GtarsPyStatBackend from bedboss.bedstat.backends.r_backend import RStatBackend -from bedboss.const import BACKEND_GTARS, BACKEND_R +from bedboss.const import BACKEND_GTARS, BACKEND_GTARS_PY, BACKEND_R -__all__ = ["StatBackend", "RStatBackend", "GtarsStatBackend", "create_backend"] +__all__ = [ + "StatBackend", + "RStatBackend", + "GtarsStatBackend", + "GtarsPyStatBackend", + "create_backend", +] def create_backend(name: str, **kwargs) -> StatBackend: """Create a statistics computation backend by name. - :param name: Backend name (BACKEND_R or BACKEND_GTARS) + :param name: Backend name (BACKEND_R, BACKEND_GTARS, BACKEND_GTARS_PY) :param kwargs: Backend-specific keyword arguments :return: StatBackend instance """ @@ -17,5 +24,9 @@ def create_backend(name: str, **kwargs) -> StatBackend: return RStatBackend(**kwargs) elif name == BACKEND_GTARS: return GtarsStatBackend(**kwargs) + elif name == BACKEND_GTARS_PY: + return GtarsPyStatBackend(**kwargs) else: - raise ValueError(f"Unknown analysis backend: {name!r}. Use 'r' or 'gtars'.") + raise ValueError( + f"Unknown analysis backend: {name!r}. Use 'r', 'gtars', or 'gtars-py'." + ) diff --git a/bedboss/bedstat/backends/gtars_py_backend.py b/bedboss/bedstat/backends/gtars_py_backend.py new file mode 100644 index 00000000..f7b032c9 --- /dev/null +++ b/bedboss/bedstat/backends/gtars_py_backend.py @@ -0,0 +1,465 @@ +"""gtars Python-bindings-direct statistics backend. + +Parallel implementation of GtarsStatBackend that skips the gtars CLI +subprocess and calls gtars Python bindings directly. Caches reference +data (FASTA, GeneModel, PartitionList, SignalMatrix, chrom_sizes) per +backend instance so batch processing amortizes load cost across files. + +Produces the same output dict shape as GtarsStatBackend for downstream +compatibility (compress_distributions + DB insertion unchanged). + +This coexists with GtarsStatBackend during performance evaluation. +After benchmarking, only one backend will remain. +""" + +import logging +import os +import statistics +from dataclasses import dataclass +from pathlib import Path +from typing import Optional, Union + +import pypiper +from gtars.models import ( + GeneModel, + GenomeAssembly, + GenomicDistAnnotation, + PartitionList, + RegionSet, + SignalMatrix, + TssIndex, +) +from gtars.genomic_distributions import ( + calc_expected_partitions, + calc_gc_content, + calc_partitions, + calc_summary_signal, + median_abs_distance, +) + +from bedboss.bedstat.backends.base import StatBackend +from bedboss.bedstat.compress_distributions import ( + compress_distributions, + compress_to_kde, +) +from bedboss.bedstat.gc_content import get_genome_assembly_obj +from bedboss.bedstat.ref_utils import ( + get_chrom_sizes_path, + get_gda_path, + get_osm_path_with_precompile, +) +from bedboss.const import BEDSTAT_OUTPUT, DEFAULT_PRECISION, OUTPUT_FOLDER_NAME +from bedboss.exceptions import BedBossException + +_LOGGER = logging.getLogger("bedboss") + +# Map gtars partition names to legacy DB column name prefixes +PARTITION_NAME_MAP = { + "promoterCore": "promotercore", + "promoterProx": "promoterprox", + "threeUTR": "threeutr", + "fiveUTR": "fiveutr", + "exon": "exon", + "intron": "intron", + "intergenic": "intergenic", +} + + +def round_floats(obj, precision: int = 4): + """Recursively round all floats in a nested dict/list structure.""" + if isinstance(obj, float): + return round(obj, precision) + elif isinstance(obj, dict): + return {k: round_floats(v, precision) for k, v in obj.items()} + elif isinstance(obj, list): + return [round_floats(v, precision) for v in obj] + return obj + + +def _normalize_partitions(raw: Optional[dict]) -> Optional[dict]: + """Convert Python binding's partitions dict to CLI JSON schema. + + Python binding: ``{"partition": [names], "count": [counts], "total": int}`` + CLI schema: ``{"counts": [[name, count], ...], "total": int}`` + """ + if not raw: + return None + names = raw.get("partition") or [] + counts = raw.get("count") or [] + return { + "counts": list(zip(names, counts)), + "total": raw.get("total", 0), + } + + +def _normalize_expected_partitions(raw: Optional[dict]) -> Optional[dict]: + """Pass-through for now — bedboss doesn't consume expected_partitions + scalars beyond storing the blob, so we leave the Python binding's + native shape intact. + """ + return raw or None + + +def _parse_chrom_sizes(path: str) -> dict: + """Parse a chrom.sizes TSV file into a dict {chrom: length}.""" + sizes = {} + with open(path, "r") as f: + for line in f: + parts = line.strip().split("\t") + if len(parts) >= 2: + try: + sizes[parts[0]] = int(parts[1]) + except ValueError: + continue + return sizes + + +@dataclass +class GenomeRefs: + """Cached reference data for a genome. + + Fields are independently optional — any failure to load a given ref + is logged and the field stays None, matching the per-file graceful + degradation of the CLI-based backend. + """ + + assembly: Optional[GenomeAssembly] = None + chrom_sizes: Optional[dict] = None + gene_model: Optional[GeneModel] = None + partition_list: Optional[PartitionList] = None + tss_index: Optional[TssIndex] = None + signal_matrix: Optional[SignalMatrix] = None + + @classmethod + def load( + cls, + genome: str, + rfg_config: Optional[str], + ensdb: Optional[str], + signal_matrix_path: Optional[str], + promoter_upstream: int, + promoter_downstream: int, + ) -> "GenomeRefs": + refs = cls() + + # chrom_sizes (lightweight, try first) + try: + cs_path = get_chrom_sizes_path(genome, rfg_config=rfg_config) + if cs_path: + refs.chrom_sizes = _parse_chrom_sizes(cs_path) + except Exception as e: + _LOGGER.warning(f"chrom_sizes unavailable for {genome}: {e}") + + # GenomeAssembly (FASTA-backed, for GC content) — leverages module-level + # cache in bedboss.bedstat.gc_content + try: + refs.assembly = get_genome_assembly_obj(genome, rfg_config=rfg_config) + except Exception as e: + _LOGGER.warning(f"GenomeAssembly unavailable for {genome}: {e}") + + # Gene model (GTF or GDA .bin) + derived PartitionList + TssIndex + if not ensdb: + try: + ensdb = get_gda_path(genome, rfg_config=rfg_config) + except Exception: + pass + + if ensdb: + try: + if str(ensdb).endswith(".bin"): + # GDA binary: has gene_model + derived tss_index + helper + # partition_list() method + gda = GenomicDistAnnotation.load_bin(str(ensdb)) + refs.gene_model = gda.gene_model() + refs.tss_index = gda.tss_index() + refs.partition_list = gda.partition_list( + promoter_upstream, + promoter_downstream, + refs.chrom_sizes, + ) + else: + # Raw GTF: load GeneModel, derive PartitionList manually, + # TssIndex deferred (would need genes+strands — not in + # Python binding's GeneModel API yet, skip for now) + refs.gene_model = GeneModel.from_gtf(str(ensdb), True, True) + try: + refs.partition_list = PartitionList.from_gene_model( + refs.gene_model, + promoter_upstream, + promoter_downstream, + refs.chrom_sizes, + ) + except Exception as e: + _LOGGER.warning(f"PartitionList from GTF failed: {e}") + except Exception as e: + _LOGGER.warning(f"Gene model load failed: {e}") + + # Signal matrix (optional) + if signal_matrix_path and os.path.exists(signal_matrix_path): + try: + resolved = get_osm_path_with_precompile(signal_matrix_path) + if resolved.endswith(".bin"): + refs.signal_matrix = SignalMatrix.load_bin(resolved) + else: + refs.signal_matrix = SignalMatrix.from_tsv(resolved) + except Exception as e: + _LOGGER.warning(f"SignalMatrix load failed: {e}") + + return refs + + +class GtarsPyStatBackend(StatBackend): + """gtars Python-bindings-direct statistics backend (no CLI subprocess).""" + + def __init__( + self, + region_dist_bins: int = 250, + promoter_upstream: int = 200, + promoter_downstream: int = 2000, + precision: int = DEFAULT_PRECISION, + **kwargs, + ): + self._region_dist_bins = region_dist_bins + self._promoter_upstream = promoter_upstream + self._promoter_downstream = promoter_downstream + self._precision = precision + # Instance-level cache keyed by (genome, ensdb, signal_matrix_path) + self._ref_cache: dict = {} + + def _get_refs( + self, + genome: str, + rfg_config: Optional[str], + ensdb: Optional[str], + signal_matrix_path: Optional[str], + ) -> GenomeRefs: + cache_key = (genome, ensdb, signal_matrix_path) + if cache_key not in self._ref_cache: + self._ref_cache[cache_key] = GenomeRefs.load( + genome=genome, + rfg_config=rfg_config, + ensdb=ensdb, + signal_matrix_path=signal_matrix_path, + promoter_upstream=self._promoter_upstream, + promoter_downstream=self._promoter_downstream, + ) + return self._ref_cache[cache_key] + + def compute( + self, + bedfile: str, + genome: str, + outfolder: str, + bed_digest: str = None, + ensdb: str = None, + open_signal_matrix: str = None, + just_db_commit: bool = False, + rfg_config: Union[str, Path] = None, + pm: pypiper.PipelineManager = None, + ) -> dict: + refs = self._get_refs(genome, rfg_config, ensdb, open_signal_matrix) + + outfolder_stats = os.path.join(outfolder, OUTPUT_FOLDER_NAME, BEDSTAT_OUTPUT) + os.makedirs(outfolder_stats, exist_ok=True) + + stop_pipeline = not pm + + bed_object = RegionSet(bedfile) + if not bed_digest: + bed_digest = bed_object.identifier + + outfolder_stats_results = os.path.abspath( + os.path.join(outfolder_stats, bed_digest) + ) + os.makedirs(outfolder_stats_results, exist_ok=True) + + if not pm: + pm_out_path = os.path.abspath( + os.path.join(outfolder_stats, "pypiper", bed_digest) + ) + os.makedirs(pm_out_path, exist_ok=True) + pm = pypiper.PipelineManager( + name="bedstat-pipeline", + outfolder=pm_out_path, + pipestat_sample_name=bed_digest, + ) + + if not just_db_commit: + try: + gtars_output = self._compute_all(bed_object, refs, pm) + except Exception as e: + _LOGGER.error(f"gtars-py compute failed: {e}") + raise BedBossException(f"gtars-py compute failed: {e}") + else: + gtars_output = {} + + # Extract scalars to flat dict keys + data = {} + scalars = gtars_output.get("scalars", {}) + data["number_of_regions"] = scalars.get("number_of_regions") + data["mean_region_width"] = scalars.get("mean_region_width") + data["median_tss_dist"] = scalars.get("median_tss_dist") + + # Derive median_neighbor_distance from the raw neighbor_distances list + neighbor_distances = gtars_output.get("distributions", {}).get( + "neighbor_distances" + ) + if neighbor_distances: + abs_vals = [abs(d) for d in neighbor_distances if d is not None] + if abs_vals: + data["median_neighbor_distance"] = round(statistics.median(abs_vals), 4) + else: + data["median_neighbor_distance"] = None + else: + data["median_neighbor_distance"] = None + + # Populate legacy partition flat columns + partitions = gtars_output.get("partitions") + if partitions: + total = partitions.get("total", 0) + for name, count in partitions.get("counts", []): + db_name = PARTITION_NAME_MAP.get(name) + if db_name and total > 0: + data[f"{db_name}_frequency"] = count + data[f"{db_name}_percentage"] = round(count / total, 4) + + # GC content: computed inside _compute_all, add mean as a scalar here + gc_contents = gtars_output.pop("_gc_contents", None) + if gc_contents: + gc_mean = round(statistics.mean(gc_contents), 4) + data["gc_content"] = gc_mean + gc_kde = compress_to_kde(gc_contents, n_points=512, log_transform=False) + if gc_kde: + gc_kde["mean"] = gc_mean + if "distributions" not in gtars_output: + gtars_output["distributions"] = {} + gtars_output["distributions"]["gc_content"] = gc_kde + else: + data["gc_content"] = None + + # Compress distributions for DB storage + compress_distributions(gtars_output) + + # Store entire augmented gtars output as distributions blob + data["distributions"] = gtars_output + + if self._precision is not None: + data = round_floats(data, self._precision) + + if stop_pipeline and pm: + pm.stop_pipeline() + + return data + + def _compute_all( + self, + rs: RegionSet, + refs: GenomeRefs, + pm: pypiper.PipelineManager, + ) -> dict: + """Run all gtars statistics via Python bindings, return raw output dict + matching the CLI JSON schema (pre-compression).""" + pm.timestamp("### Computing core statistics") + widths = rs.widths() + neighbor_distances = rs.neighbor_distances() + nearest_neighbors = rs.nearest_neighbors() + chrom_stats_obj = rs.chromosome_statistics() + + # Serialize ChromosomeStatistics objects to dicts + chromosome_stats = {} + for chrom, stats_obj in chrom_stats_obj.items(): + chromosome_stats[chrom] = { + "chromosome": stats_obj.chromosome, + "number_of_regions": stats_obj.number_of_regions, + "start_nucleotide_position": stats_obj.start_nucleotide_position, + "end_nucleotide_position": stats_obj.end_nucleotide_position, + "minimum_region_length": stats_obj.minimum_region_length, + "maximum_region_length": stats_obj.maximum_region_length, + "mean_region_length": stats_obj.mean_region_length, + "median_region_length": stats_obj.median_region_length, + } + + pm.timestamp("### Computing region distribution") + region_distribution = rs.distribution( + n_bins=self._region_dist_bins, + chrom_sizes=refs.chrom_sizes, + ) + + number_of_regions = len(widths) + mean_region_width = ( + sum(widths) / number_of_regions if number_of_regions > 0 else 0.0 + ) + + # TSS distances via cached TssIndex + tss_distances = None + median_tss_dist = None + if refs.tss_index is not None: + pm.timestamp("### Computing TSS distances") + try: + tss_distances = refs.tss_index.feature_distances(rs) + median_tss_dist = median_abs_distance( + [float(d) for d in tss_distances if d is not None] + ) + except Exception as e: + _LOGGER.warning(f"TSS distance computation failed: {e}") + + # Partitions + expected partitions. + # The Python binding returns {partition: [names], count: [counts], total: n} + # but the CLI JSON schema (that downstream code expects) has + # {counts: [[name, count], ...], total: n}. Normalize here. + partitions = None + expected_partitions = None + if refs.partition_list is not None: + pm.timestamp("### Computing partitions") + try: + raw = calc_partitions(rs, refs.partition_list, False) + partitions = _normalize_partitions(raw) + except Exception as e: + _LOGGER.warning(f"Partition classification failed: {e}") + if refs.partition_list is not None and refs.chrom_sizes is not None: + try: + raw_ep = calc_expected_partitions( + rs, refs.partition_list, refs.chrom_sizes, False + ) + expected_partitions = _normalize_expected_partitions(raw_ep) + except Exception as e: + _LOGGER.warning(f"Expected partitions failed: {e}") + + # Signal matrix overlap + open_signal = None + if refs.signal_matrix is not None: + pm.timestamp("### Computing open chromatin signal") + try: + open_signal = calc_summary_signal(rs, refs.signal_matrix) + except Exception as e: + _LOGGER.warning(f"Signal summary failed: {e}") + + # GC content (passed through as _gc_contents for outer fn to handle) + gc_contents = None + if refs.assembly is not None: + pm.timestamp("### Computing GC content") + try: + gc_contents = calc_gc_content(rs, refs.assembly, ignore_unk_chroms=True) + except Exception as e: + _LOGGER.warning(f"GC content failed: {e}") + + # Assemble output matching CLI JSON schema + return { + "scalars": { + "number_of_regions": number_of_regions, + "mean_region_width": mean_region_width, + "median_tss_dist": median_tss_dist, + }, + "partitions": partitions, + "distributions": { + "widths": widths, + "tss_distances": tss_distances, + "neighbor_distances": neighbor_distances, + "nearest_neighbors": nearest_neighbors, + "region_distribution": region_distribution, + "chromosome_stats": chromosome_stats, + }, + "expected_partitions": expected_partitions, + "open_signal": open_signal, + "_gc_contents": gc_contents, # consumed by compute() + } diff --git a/bedboss/cli.py b/bedboss/cli.py index 1816bf76..6f73bd6f 100644 --- a/bedboss/cli.py +++ b/bedboss/cli.py @@ -101,7 +101,7 @@ def run_all( upload_pephub: bool = typer.Option(False, help="Upload to PEPHub"), backend: str = typer.Option( None, - help="Override analysis backend ('r' or 'gtars'). If not set, uses config file value.", + help="Override analysis backend ('r', 'gtars', or 'gtars-py'). If not set, uses config file value.", ), # Universes universe: bool = typer.Option(False, help="Create a universe"), diff --git a/bedboss/const.py b/bedboss/const.py index 9f7c2ecd..64a4d5fb 100644 --- a/bedboss/const.py +++ b/bedboss/const.py @@ -30,6 +30,7 @@ # bedstat BACKEND_R = "r" BACKEND_GTARS = "gtars" +BACKEND_GTARS_PY = "gtars-py" DEFAULT_PRECISION = 3 # bedbuncher From 1615e8e737341800933266d1f28750f5acce5652 Mon Sep 17 00:00:00 2001 From: Sam Park Date: Sun, 5 Apr 2026 17:14:49 -0400 Subject: [PATCH 07/10] Make GtarsStatBackend pure CLI subprocess (no Python binding calls) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Before: GtarsStatBackend was "mixed" — called gtars CLI for most stats but used gtars Python bindings for two things inside the backend: 1. RegionSet(bedfile) to compute bed_digest fallback 2. calc_gc_content(rs, assembly, ...) for GC content This blurred the line between the "CLI" and "Python bindings" backends. After: GtarsStatBackend is a pure CLI subprocess backend. - Pass --fasta to gtars CLI (uses PR #249's --fasta flag that outputs gc_content in the JSON). Read gc_mean + per_region GC from the CLI output instead of calling gtars Python bindings. - Require bed_digest from the caller. Backends no longer parse BED files. Resolution moves up to bedstat() (the orchestration layer), which computes the digest via gtars Python bindings once when absent. - Remove imports of gtars.models.RegionSet and calculate_gc_content. - Add get_fasta_path() helper to ref_utils.py (pure refgenie, no gtars dependency at import time) so the backend can resolve FASTA paths for the --fasta flag. The three backends now have cleanly separated execution domains: - RStatBackend: R subprocess (+ Python bindings helpers for file loading and GC content — kept as-is, R's native GC calc is slow) - GtarsStatBackend: gtars CLI subprocess ONLY - GtarsPyStatBackend: gtars Python bindings ONLY (no subprocess) Parity verified: pure-CLI gtars output matches gtars-py output exactly on all 19 scalar + partition fields for a 126K-region hg38 file. Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedstat/backends/gtars_backend.py | 69 +++++++++++++---------- bedboss/bedstat/bedstat.py | 9 +++ bedboss/bedstat/ref_utils.py | 30 +++++++++- 3 files changed, 78 insertions(+), 30 deletions(-) diff --git a/bedboss/bedstat/backends/gtars_backend.py b/bedboss/bedstat/backends/gtars_backend.py index 0914c790..81e6c24b 100644 --- a/bedboss/bedstat/backends/gtars_backend.py +++ b/bedboss/bedstat/backends/gtars_backend.py @@ -1,3 +1,15 @@ +"""gtars genomicdist CLI backend (pure subprocess, no Python bindings). + +Invokes the `gtars genomicdist` CLI as a subprocess per file. All +computation — scalars, region distribution, GC content, dinucleotide +frequencies, partitions, signal matrix overlap — happens inside the +CLI binary. This backend only orchestrates the subprocess and parses +the JSON output. + +For the Python-bindings counterpart (no subprocess, in-process gtars +calls) see GtarsPyStatBackend. +""" + import json import logging import os @@ -6,16 +18,15 @@ from typing import Union import pypiper -from gtars.models import RegionSet from bedboss.bedstat.backends.base import StatBackend from bedboss.bedstat.compress_distributions import ( compress_distributions, compress_to_kde, ) -from bedboss.bedstat.gc_content import calculate_gc_content from bedboss.bedstat.ref_utils import ( get_chrom_sizes_path, + get_fasta_path, get_gda_path, get_osm_path_with_precompile, ) @@ -48,7 +59,7 @@ def round_floats(obj, precision: int = 4): class GtarsStatBackend(StatBackend): - """gtars genomicdist-based statistics backend.""" + """gtars genomicdist CLI backend — pure subprocess, no Python bindings.""" def __init__( self, @@ -75,6 +86,13 @@ def compute( rfg_config: Union[str, Path] = None, pm: pypiper.PipelineManager = None, ) -> dict: + if bed_digest is None: + raise BedBossException( + "GtarsStatBackend.compute() requires bed_digest; the backend " + "does not parse BED files in Python. Pass bed_digest from the " + "orchestrator (bedstat() resolves it when absent)." + ) + # Auto-fetch GDA/GTF annotation via refgenie if not provided if not ensdb: try: @@ -99,16 +117,15 @@ def compute( "Region distribution will not be normalized." ) + # Resolve FASTA path so the CLI can compute GC content + dinucl freq + fasta_path = get_fasta_path(genome, rfg_config=rfg_config) + outfolder_stats = os.path.join(outfolder, OUTPUT_FOLDER_NAME, BEDSTAT_OUTPUT) os.makedirs(outfolder_stats, exist_ok=True) # Used to stop pipeline if bedstat is used independently stop_pipeline = not pm - bed_object = RegionSet(bedfile) - if not bed_digest: - bed_digest = bed_object.identifier - outfolder_stats_results = os.path.abspath( os.path.join(outfolder_stats, bed_digest) ) @@ -145,6 +162,8 @@ def compute( cmd_parts.extend(["--chrom-sizes", chrom_sizes]) if open_signal_matrix: cmd_parts.extend(["--signal-matrix", open_signal_matrix]) + if fasta_path: + cmd_parts.extend(["--fasta", fasta_path, "--ignore-unk-chroms"]) cmd_parts.extend( [ "--bins", @@ -204,29 +223,21 @@ def compute( data[f"{db_name}_frequency"] = count data[f"{db_name}_percentage"] = round(count / total, 4) - # GC content: compute via Python bindings (requires refgenie FASTA) - try: - gc_contents = calculate_gc_content( - bedfile=bed_object, genome=genome, rfg_config=rfg_config - ) - except BaseException as e: - _LOGGER.warning( - f"GC content calculation skipped for {genome}: {e}. " - "Ensure refgenie is configured with a FASTA asset." - ) - gc_contents = None - - if gc_contents: - gc_mean = round(statistics.mean(gc_contents), 4) + # GC content: read from CLI output (computed by `gtars genomicdist --fasta`) + gc_block = gtars_output.get("gc_content") + if gc_block: + gc_mean = round(gc_block.get("mean", 0.0), 4) data["gc_content"] = gc_mean - - # Compress per-region GC values to 512-pt KDE, inject into distributions - gc_kde = compress_to_kde(gc_contents, n_points=512, log_transform=False) - if gc_kde: - gc_kde["mean"] = gc_mean - if "distributions" not in gtars_output: - gtars_output["distributions"] = {} - gtars_output["distributions"]["gc_content"] = gc_kde + gc_per_region = gc_block.get("per_region") or [] + if gc_per_region: + gc_kde = compress_to_kde( + gc_per_region, n_points=512, log_transform=False + ) + if gc_kde: + gc_kde["mean"] = gc_mean + if "distributions" not in gtars_output: + gtars_output["distributions"] = {} + gtars_output["distributions"]["gc_content"] = gc_kde else: data["gc_content"] = None diff --git a/bedboss/bedstat/bedstat.py b/bedboss/bedstat/bedstat.py index 1565b270..3e8239bc 100755 --- a/bedboss/bedstat/bedstat.py +++ b/bedboss/bedstat/bedstat.py @@ -105,6 +105,15 @@ def bedstat( f"Open Signal Matrix was not found for {genome}. Skipping..." ) + # Resolve bed_digest at the orchestration layer. Backends assume the + # digest is always provided — this keeps backend implementations free + # of file-parsing concerns (especially GtarsStatBackend, which is + # a pure CLI-subprocess backend with no Python-binding dependencies). + if bed_digest is None: + from gtars.models import RegionSet + + bed_digest = RegionSet(bedfile).identifier + return backend.compute( bedfile=bedfile, genome=genome, diff --git a/bedboss/bedstat/ref_utils.py b/bedboss/bedstat/ref_utils.py index 7c47cf71..aca1d576 100644 --- a/bedboss/bedstat/ref_utils.py +++ b/bedboss/bedstat/ref_utils.py @@ -1,7 +1,7 @@ """Reference file resolution for gtars genomicdist. Handles auto-fetching and pre-compilation of GTF annotations, chrom.sizes, -and open signal matrices via refgenie with seqcol API fallback. +open signal matrices, and FASTA files via refgenie with seqcol API fallback. """ import logging @@ -14,6 +14,34 @@ _LOGGER = logging.getLogger("bedboss") +def get_fasta_path(genome: str, rfg_config: str = None) -> Union[str, None]: + """Return path to the FASTA file for the given genome via refgenie. + + Pure refgenie lookup — no gtars dependency. Used by the gtars CLI + backend to get a FASTA path to pass to `gtars genomicdist --fasta`. + + :param genome: genome assembly name + :param rfg_config: path to refgenie config file (optional) + :return: path to the FASTA file, or None if unavailable + """ + from refgenconf import RefgenconfError + from yacman.exceptions import UndefinedAliasError + + from bedboss.bedmaker.utils import get_rgc + + try: + rgc = get_rgc(rfg_config=rfg_config) + return rgc.seek( + genome_name=genome, + asset_name="fasta", + tag_name="default", + seek_key="fasta", + ) + except (UndefinedAliasError, RefgenconfError) as e: + _LOGGER.warning(f"Could not resolve FASTA for {genome}: {e}") + return None + + def _get_chrom_sizes_seqcol(genome: str) -> Union[str, None]: """Fallback: fetch chrom.sizes via seqcol API when refgenie doesn't have it. From 3383efe204db545b7a4c5be04ed9534eb6a1dd2f Mon Sep 17 00:00:00 2001 From: Sam Park Date: Mon, 6 Apr 2026 02:40:16 -0400 Subject: [PATCH 08/10] Remove gtars-py backend, add .fab auto-compilation, fix backend wiring - Remove GtarsPyStatBackend (gtars CLI with .fab is faster and simpler) - Remove BACKEND_GTARS_PY constant and all references - Add get_fab_path() to ref_utils: auto-compiles .fab from FASTA on first use via gtars prep --fasta, cached forever - GtarsStatBackend prefers .fab over plain .fa for GC content - Fix reprocess_bedset: pass backend from config to run_bedbuncher - Update CLI help text and bedbuncher docstrings Two backends remain: "r" (RStatBackend) and "gtars" (GtarsStatBackend). Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedboss.py | 1 + bedboss/bedbuncher/bedbuncher.py | 4 +- bedboss/bedstat/backends/__init__.py | 12 +- bedboss/bedstat/backends/gtars_backend.py | 11 +- bedboss/bedstat/backends/gtars_py_backend.py | 465 ------------------- bedboss/bedstat/ref_utils.py | 40 ++ bedboss/cli.py | 4 +- bedboss/const.py | 1 - 8 files changed, 56 insertions(+), 482 deletions(-) delete mode 100644 bedboss/bedstat/backends/gtars_py_backend.py diff --git a/bedboss/bedboss.py b/bedboss/bedboss.py index cff66330..4ff77164 100644 --- a/bedboss/bedboss.py +++ b/bedboss/bedboss.py @@ -786,4 +786,5 @@ def reprocess_bedset( ) }, lite=False, + backend=bbagent.config.config.analysis.backend, ) diff --git a/bedboss/bedbuncher/bedbuncher.py b/bedboss/bedbuncher/bedbuncher.py index 0a1c700f..6ae814c4 100644 --- a/bedboss/bedbuncher/bedbuncher.py +++ b/bedboss/bedbuncher/bedbuncher.py @@ -122,7 +122,7 @@ def run_bedbuncher( :param upload_s3: whether to upload files to s3 :param force_overwrite: whether to overwrite the record in the database :param lite: whether to run the pipeline in lite mode - :param backend: analysis backend ("r", "gtars", or "gtars-py"). For gtars variants, heavy is ignored. + :param backend: analysis backend ("r" or "gtars"). For gtars, heavy is ignored. :return: """ _LOGGER.info(f"Adding bedset {record_id} to the database") @@ -142,7 +142,7 @@ def run_bedbuncher( "bedsets", ) - gtars_like = backend in ("gtars", "gtars-py") + gtars_like = backend == "gtars" if heavy and not gtars_like: _LOGGER.info("Heavy processing is True. Calculating plots...") plot_value = create_plots( diff --git a/bedboss/bedstat/backends/__init__.py b/bedboss/bedstat/backends/__init__.py index 066b2c38..0d28cdbf 100644 --- a/bedboss/bedstat/backends/__init__.py +++ b/bedboss/bedstat/backends/__init__.py @@ -1,14 +1,12 @@ from bedboss.bedstat.backends.base import StatBackend from bedboss.bedstat.backends.gtars_backend import GtarsStatBackend -from bedboss.bedstat.backends.gtars_py_backend import GtarsPyStatBackend from bedboss.bedstat.backends.r_backend import RStatBackend -from bedboss.const import BACKEND_GTARS, BACKEND_GTARS_PY, BACKEND_R +from bedboss.const import BACKEND_GTARS, BACKEND_R __all__ = [ "StatBackend", "RStatBackend", "GtarsStatBackend", - "GtarsPyStatBackend", "create_backend", "build_backend", ] @@ -21,7 +19,7 @@ def create_backend(name: str, **kwargs) -> StatBackend: should use :func:`build_backend` instead, which handles backend-specific prerequisites (e.g. starting an RServiceManager for the R backend). - :param name: Backend name (BACKEND_R, BACKEND_GTARS, BACKEND_GTARS_PY) + :param name: Backend name ('r' or 'gtars') :param kwargs: Backend-specific keyword arguments :return: StatBackend instance """ @@ -29,11 +27,9 @@ def create_backend(name: str, **kwargs) -> StatBackend: return RStatBackend(**kwargs) elif name == BACKEND_GTARS: return GtarsStatBackend(**kwargs) - elif name == BACKEND_GTARS_PY: - return GtarsPyStatBackend(**kwargs) else: raise ValueError( - f"Unknown analysis backend: {name!r}. Use 'r', 'gtars', or 'gtars-py'." + f"Unknown analysis backend: {name!r}. Use 'r' or 'gtars'." ) @@ -48,7 +44,7 @@ def build_backend(name: str) -> StatBackend: done to release backend-held resources. See `StatBackend` as a context manager for automatic cleanup. - :param name: Backend name (BACKEND_R, BACKEND_GTARS, BACKEND_GTARS_PY) + :param name: Backend name ('r' or 'gtars') :return: StatBackend instance ready for batch processing """ if name == BACKEND_R: diff --git a/bedboss/bedstat/backends/gtars_backend.py b/bedboss/bedstat/backends/gtars_backend.py index 81e6c24b..2a11ee83 100644 --- a/bedboss/bedstat/backends/gtars_backend.py +++ b/bedboss/bedstat/backends/gtars_backend.py @@ -6,8 +6,8 @@ CLI binary. This backend only orchestrates the subprocess and parses the JSON output. -For the Python-bindings counterpart (no subprocess, in-process gtars -calls) see GtarsPyStatBackend. +Uses .fab binary FASTA format when available for optimal GC content +performance (zero-copy mmap). Falls back to plain FASTA otherwise. """ import json @@ -26,6 +26,7 @@ ) from bedboss.bedstat.ref_utils import ( get_chrom_sizes_path, + get_fab_path, get_fasta_path, get_gda_path, get_osm_path_with_precompile, @@ -117,8 +118,10 @@ def compute( "Region distribution will not be normalized." ) - # Resolve FASTA path so the CLI can compute GC content + dinucl freq - fasta_path = get_fasta_path(genome, rfg_config=rfg_config) + # Resolve FASTA path — prefer .fab (zero-copy mmap) over plain .fa + fasta_path = get_fab_path(genome, rfg_config=rfg_config) + if not fasta_path: + fasta_path = get_fasta_path(genome, rfg_config=rfg_config) outfolder_stats = os.path.join(outfolder, OUTPUT_FOLDER_NAME, BEDSTAT_OUTPUT) os.makedirs(outfolder_stats, exist_ok=True) diff --git a/bedboss/bedstat/backends/gtars_py_backend.py b/bedboss/bedstat/backends/gtars_py_backend.py deleted file mode 100644 index f7b032c9..00000000 --- a/bedboss/bedstat/backends/gtars_py_backend.py +++ /dev/null @@ -1,465 +0,0 @@ -"""gtars Python-bindings-direct statistics backend. - -Parallel implementation of GtarsStatBackend that skips the gtars CLI -subprocess and calls gtars Python bindings directly. Caches reference -data (FASTA, GeneModel, PartitionList, SignalMatrix, chrom_sizes) per -backend instance so batch processing amortizes load cost across files. - -Produces the same output dict shape as GtarsStatBackend for downstream -compatibility (compress_distributions + DB insertion unchanged). - -This coexists with GtarsStatBackend during performance evaluation. -After benchmarking, only one backend will remain. -""" - -import logging -import os -import statistics -from dataclasses import dataclass -from pathlib import Path -from typing import Optional, Union - -import pypiper -from gtars.models import ( - GeneModel, - GenomeAssembly, - GenomicDistAnnotation, - PartitionList, - RegionSet, - SignalMatrix, - TssIndex, -) -from gtars.genomic_distributions import ( - calc_expected_partitions, - calc_gc_content, - calc_partitions, - calc_summary_signal, - median_abs_distance, -) - -from bedboss.bedstat.backends.base import StatBackend -from bedboss.bedstat.compress_distributions import ( - compress_distributions, - compress_to_kde, -) -from bedboss.bedstat.gc_content import get_genome_assembly_obj -from bedboss.bedstat.ref_utils import ( - get_chrom_sizes_path, - get_gda_path, - get_osm_path_with_precompile, -) -from bedboss.const import BEDSTAT_OUTPUT, DEFAULT_PRECISION, OUTPUT_FOLDER_NAME -from bedboss.exceptions import BedBossException - -_LOGGER = logging.getLogger("bedboss") - -# Map gtars partition names to legacy DB column name prefixes -PARTITION_NAME_MAP = { - "promoterCore": "promotercore", - "promoterProx": "promoterprox", - "threeUTR": "threeutr", - "fiveUTR": "fiveutr", - "exon": "exon", - "intron": "intron", - "intergenic": "intergenic", -} - - -def round_floats(obj, precision: int = 4): - """Recursively round all floats in a nested dict/list structure.""" - if isinstance(obj, float): - return round(obj, precision) - elif isinstance(obj, dict): - return {k: round_floats(v, precision) for k, v in obj.items()} - elif isinstance(obj, list): - return [round_floats(v, precision) for v in obj] - return obj - - -def _normalize_partitions(raw: Optional[dict]) -> Optional[dict]: - """Convert Python binding's partitions dict to CLI JSON schema. - - Python binding: ``{"partition": [names], "count": [counts], "total": int}`` - CLI schema: ``{"counts": [[name, count], ...], "total": int}`` - """ - if not raw: - return None - names = raw.get("partition") or [] - counts = raw.get("count") or [] - return { - "counts": list(zip(names, counts)), - "total": raw.get("total", 0), - } - - -def _normalize_expected_partitions(raw: Optional[dict]) -> Optional[dict]: - """Pass-through for now — bedboss doesn't consume expected_partitions - scalars beyond storing the blob, so we leave the Python binding's - native shape intact. - """ - return raw or None - - -def _parse_chrom_sizes(path: str) -> dict: - """Parse a chrom.sizes TSV file into a dict {chrom: length}.""" - sizes = {} - with open(path, "r") as f: - for line in f: - parts = line.strip().split("\t") - if len(parts) >= 2: - try: - sizes[parts[0]] = int(parts[1]) - except ValueError: - continue - return sizes - - -@dataclass -class GenomeRefs: - """Cached reference data for a genome. - - Fields are independently optional — any failure to load a given ref - is logged and the field stays None, matching the per-file graceful - degradation of the CLI-based backend. - """ - - assembly: Optional[GenomeAssembly] = None - chrom_sizes: Optional[dict] = None - gene_model: Optional[GeneModel] = None - partition_list: Optional[PartitionList] = None - tss_index: Optional[TssIndex] = None - signal_matrix: Optional[SignalMatrix] = None - - @classmethod - def load( - cls, - genome: str, - rfg_config: Optional[str], - ensdb: Optional[str], - signal_matrix_path: Optional[str], - promoter_upstream: int, - promoter_downstream: int, - ) -> "GenomeRefs": - refs = cls() - - # chrom_sizes (lightweight, try first) - try: - cs_path = get_chrom_sizes_path(genome, rfg_config=rfg_config) - if cs_path: - refs.chrom_sizes = _parse_chrom_sizes(cs_path) - except Exception as e: - _LOGGER.warning(f"chrom_sizes unavailable for {genome}: {e}") - - # GenomeAssembly (FASTA-backed, for GC content) — leverages module-level - # cache in bedboss.bedstat.gc_content - try: - refs.assembly = get_genome_assembly_obj(genome, rfg_config=rfg_config) - except Exception as e: - _LOGGER.warning(f"GenomeAssembly unavailable for {genome}: {e}") - - # Gene model (GTF or GDA .bin) + derived PartitionList + TssIndex - if not ensdb: - try: - ensdb = get_gda_path(genome, rfg_config=rfg_config) - except Exception: - pass - - if ensdb: - try: - if str(ensdb).endswith(".bin"): - # GDA binary: has gene_model + derived tss_index + helper - # partition_list() method - gda = GenomicDistAnnotation.load_bin(str(ensdb)) - refs.gene_model = gda.gene_model() - refs.tss_index = gda.tss_index() - refs.partition_list = gda.partition_list( - promoter_upstream, - promoter_downstream, - refs.chrom_sizes, - ) - else: - # Raw GTF: load GeneModel, derive PartitionList manually, - # TssIndex deferred (would need genes+strands — not in - # Python binding's GeneModel API yet, skip for now) - refs.gene_model = GeneModel.from_gtf(str(ensdb), True, True) - try: - refs.partition_list = PartitionList.from_gene_model( - refs.gene_model, - promoter_upstream, - promoter_downstream, - refs.chrom_sizes, - ) - except Exception as e: - _LOGGER.warning(f"PartitionList from GTF failed: {e}") - except Exception as e: - _LOGGER.warning(f"Gene model load failed: {e}") - - # Signal matrix (optional) - if signal_matrix_path and os.path.exists(signal_matrix_path): - try: - resolved = get_osm_path_with_precompile(signal_matrix_path) - if resolved.endswith(".bin"): - refs.signal_matrix = SignalMatrix.load_bin(resolved) - else: - refs.signal_matrix = SignalMatrix.from_tsv(resolved) - except Exception as e: - _LOGGER.warning(f"SignalMatrix load failed: {e}") - - return refs - - -class GtarsPyStatBackend(StatBackend): - """gtars Python-bindings-direct statistics backend (no CLI subprocess).""" - - def __init__( - self, - region_dist_bins: int = 250, - promoter_upstream: int = 200, - promoter_downstream: int = 2000, - precision: int = DEFAULT_PRECISION, - **kwargs, - ): - self._region_dist_bins = region_dist_bins - self._promoter_upstream = promoter_upstream - self._promoter_downstream = promoter_downstream - self._precision = precision - # Instance-level cache keyed by (genome, ensdb, signal_matrix_path) - self._ref_cache: dict = {} - - def _get_refs( - self, - genome: str, - rfg_config: Optional[str], - ensdb: Optional[str], - signal_matrix_path: Optional[str], - ) -> GenomeRefs: - cache_key = (genome, ensdb, signal_matrix_path) - if cache_key not in self._ref_cache: - self._ref_cache[cache_key] = GenomeRefs.load( - genome=genome, - rfg_config=rfg_config, - ensdb=ensdb, - signal_matrix_path=signal_matrix_path, - promoter_upstream=self._promoter_upstream, - promoter_downstream=self._promoter_downstream, - ) - return self._ref_cache[cache_key] - - def compute( - self, - bedfile: str, - genome: str, - outfolder: str, - bed_digest: str = None, - ensdb: str = None, - open_signal_matrix: str = None, - just_db_commit: bool = False, - rfg_config: Union[str, Path] = None, - pm: pypiper.PipelineManager = None, - ) -> dict: - refs = self._get_refs(genome, rfg_config, ensdb, open_signal_matrix) - - outfolder_stats = os.path.join(outfolder, OUTPUT_FOLDER_NAME, BEDSTAT_OUTPUT) - os.makedirs(outfolder_stats, exist_ok=True) - - stop_pipeline = not pm - - bed_object = RegionSet(bedfile) - if not bed_digest: - bed_digest = bed_object.identifier - - outfolder_stats_results = os.path.abspath( - os.path.join(outfolder_stats, bed_digest) - ) - os.makedirs(outfolder_stats_results, exist_ok=True) - - if not pm: - pm_out_path = os.path.abspath( - os.path.join(outfolder_stats, "pypiper", bed_digest) - ) - os.makedirs(pm_out_path, exist_ok=True) - pm = pypiper.PipelineManager( - name="bedstat-pipeline", - outfolder=pm_out_path, - pipestat_sample_name=bed_digest, - ) - - if not just_db_commit: - try: - gtars_output = self._compute_all(bed_object, refs, pm) - except Exception as e: - _LOGGER.error(f"gtars-py compute failed: {e}") - raise BedBossException(f"gtars-py compute failed: {e}") - else: - gtars_output = {} - - # Extract scalars to flat dict keys - data = {} - scalars = gtars_output.get("scalars", {}) - data["number_of_regions"] = scalars.get("number_of_regions") - data["mean_region_width"] = scalars.get("mean_region_width") - data["median_tss_dist"] = scalars.get("median_tss_dist") - - # Derive median_neighbor_distance from the raw neighbor_distances list - neighbor_distances = gtars_output.get("distributions", {}).get( - "neighbor_distances" - ) - if neighbor_distances: - abs_vals = [abs(d) for d in neighbor_distances if d is not None] - if abs_vals: - data["median_neighbor_distance"] = round(statistics.median(abs_vals), 4) - else: - data["median_neighbor_distance"] = None - else: - data["median_neighbor_distance"] = None - - # Populate legacy partition flat columns - partitions = gtars_output.get("partitions") - if partitions: - total = partitions.get("total", 0) - for name, count in partitions.get("counts", []): - db_name = PARTITION_NAME_MAP.get(name) - if db_name and total > 0: - data[f"{db_name}_frequency"] = count - data[f"{db_name}_percentage"] = round(count / total, 4) - - # GC content: computed inside _compute_all, add mean as a scalar here - gc_contents = gtars_output.pop("_gc_contents", None) - if gc_contents: - gc_mean = round(statistics.mean(gc_contents), 4) - data["gc_content"] = gc_mean - gc_kde = compress_to_kde(gc_contents, n_points=512, log_transform=False) - if gc_kde: - gc_kde["mean"] = gc_mean - if "distributions" not in gtars_output: - gtars_output["distributions"] = {} - gtars_output["distributions"]["gc_content"] = gc_kde - else: - data["gc_content"] = None - - # Compress distributions for DB storage - compress_distributions(gtars_output) - - # Store entire augmented gtars output as distributions blob - data["distributions"] = gtars_output - - if self._precision is not None: - data = round_floats(data, self._precision) - - if stop_pipeline and pm: - pm.stop_pipeline() - - return data - - def _compute_all( - self, - rs: RegionSet, - refs: GenomeRefs, - pm: pypiper.PipelineManager, - ) -> dict: - """Run all gtars statistics via Python bindings, return raw output dict - matching the CLI JSON schema (pre-compression).""" - pm.timestamp("### Computing core statistics") - widths = rs.widths() - neighbor_distances = rs.neighbor_distances() - nearest_neighbors = rs.nearest_neighbors() - chrom_stats_obj = rs.chromosome_statistics() - - # Serialize ChromosomeStatistics objects to dicts - chromosome_stats = {} - for chrom, stats_obj in chrom_stats_obj.items(): - chromosome_stats[chrom] = { - "chromosome": stats_obj.chromosome, - "number_of_regions": stats_obj.number_of_regions, - "start_nucleotide_position": stats_obj.start_nucleotide_position, - "end_nucleotide_position": stats_obj.end_nucleotide_position, - "minimum_region_length": stats_obj.minimum_region_length, - "maximum_region_length": stats_obj.maximum_region_length, - "mean_region_length": stats_obj.mean_region_length, - "median_region_length": stats_obj.median_region_length, - } - - pm.timestamp("### Computing region distribution") - region_distribution = rs.distribution( - n_bins=self._region_dist_bins, - chrom_sizes=refs.chrom_sizes, - ) - - number_of_regions = len(widths) - mean_region_width = ( - sum(widths) / number_of_regions if number_of_regions > 0 else 0.0 - ) - - # TSS distances via cached TssIndex - tss_distances = None - median_tss_dist = None - if refs.tss_index is not None: - pm.timestamp("### Computing TSS distances") - try: - tss_distances = refs.tss_index.feature_distances(rs) - median_tss_dist = median_abs_distance( - [float(d) for d in tss_distances if d is not None] - ) - except Exception as e: - _LOGGER.warning(f"TSS distance computation failed: {e}") - - # Partitions + expected partitions. - # The Python binding returns {partition: [names], count: [counts], total: n} - # but the CLI JSON schema (that downstream code expects) has - # {counts: [[name, count], ...], total: n}. Normalize here. - partitions = None - expected_partitions = None - if refs.partition_list is not None: - pm.timestamp("### Computing partitions") - try: - raw = calc_partitions(rs, refs.partition_list, False) - partitions = _normalize_partitions(raw) - except Exception as e: - _LOGGER.warning(f"Partition classification failed: {e}") - if refs.partition_list is not None and refs.chrom_sizes is not None: - try: - raw_ep = calc_expected_partitions( - rs, refs.partition_list, refs.chrom_sizes, False - ) - expected_partitions = _normalize_expected_partitions(raw_ep) - except Exception as e: - _LOGGER.warning(f"Expected partitions failed: {e}") - - # Signal matrix overlap - open_signal = None - if refs.signal_matrix is not None: - pm.timestamp("### Computing open chromatin signal") - try: - open_signal = calc_summary_signal(rs, refs.signal_matrix) - except Exception as e: - _LOGGER.warning(f"Signal summary failed: {e}") - - # GC content (passed through as _gc_contents for outer fn to handle) - gc_contents = None - if refs.assembly is not None: - pm.timestamp("### Computing GC content") - try: - gc_contents = calc_gc_content(rs, refs.assembly, ignore_unk_chroms=True) - except Exception as e: - _LOGGER.warning(f"GC content failed: {e}") - - # Assemble output matching CLI JSON schema - return { - "scalars": { - "number_of_regions": number_of_regions, - "mean_region_width": mean_region_width, - "median_tss_dist": median_tss_dist, - }, - "partitions": partitions, - "distributions": { - "widths": widths, - "tss_distances": tss_distances, - "neighbor_distances": neighbor_distances, - "nearest_neighbors": nearest_neighbors, - "region_distribution": region_distribution, - "chromosome_stats": chromosome_stats, - }, - "expected_partitions": expected_partitions, - "open_signal": open_signal, - "_gc_contents": gc_contents, # consumed by compute() - } diff --git a/bedboss/bedstat/ref_utils.py b/bedboss/bedstat/ref_utils.py index aca1d576..28ac1cdf 100644 --- a/bedboss/bedstat/ref_utils.py +++ b/bedboss/bedstat/ref_utils.py @@ -42,6 +42,46 @@ def get_fasta_path(genome: str, rfg_config: str = None) -> Union[str, None]: return None +def get_fab_path(genome: str, rfg_config: str = None) -> Union[str, None]: + """Return path to a .fab binary FASTA for the given genome. + + Checks for a .fab file alongside the refgenie FASTA. If not found, + auto-compiles it via `gtars prep --fasta`. Returns the .fab path + for zero-copy mmap access in the CLI. + + :param genome: genome assembly name + :param rfg_config: path to refgenie config file (optional) + :return: path to the .fab file, or None if FASTA unavailable + """ + fasta_path = get_fasta_path(genome, rfg_config=rfg_config) + if not fasta_path: + return None + + fab_path = fasta_path + ".fab" + if os.path.exists(fab_path): + return fab_path + + # Auto-compile + _LOGGER.info(f"Compiling .fab for {genome}: {fasta_path} -> {fab_path}") + try: + result = subprocess.run( + ["gtars", "prep", "--fasta", fasta_path, "-o", fab_path], + capture_output=True, + text=True, + ) + if result.returncode == 0 and os.path.exists(fab_path): + _LOGGER.info(f"Created .fab: {fab_path}") + return fab_path + else: + _LOGGER.warning( + f"Failed to compile .fab: {result.stderr.strip()}" + ) + return None + except FileNotFoundError: + _LOGGER.warning("gtars CLI not found — cannot compile .fab") + return None + + def _get_chrom_sizes_seqcol(genome: str) -> Union[str, None]: """Fallback: fetch chrom.sizes via seqcol API when refgenie doesn't have it. diff --git a/bedboss/cli.py b/bedboss/cli.py index e3126bfb..d943e062 100644 --- a/bedboss/cli.py +++ b/bedboss/cli.py @@ -101,7 +101,7 @@ def run_all( upload_pephub: bool = typer.Option(False, help="Upload to PEPHub"), backend: str = typer.Option( None, - help="Override analysis backend ('r', 'gtars', or 'gtars-py'). If not set, uses config file value.", + help="Override analysis backend ('r' or 'gtars'). If not set, uses config file value.", ), # Universes universe: bool = typer.Option(False, help="Create a universe"), @@ -403,7 +403,7 @@ def run_stats( just_db_commit: bool = typer.Option(False, help="Just commit to the database?"), backend: str = typer.Option( "r", - help="Analysis backend ('r', 'gtars', or 'gtars-py'). Default: 'r'.", + help="Analysis backend ('r' or 'gtars'). Default: 'r'.", ), # PipelineManager multi: bool = typer.Option(False, help="Run multiple samples"), diff --git a/bedboss/const.py b/bedboss/const.py index 64a4d5fb..9f7c2ecd 100644 --- a/bedboss/const.py +++ b/bedboss/const.py @@ -30,7 +30,6 @@ # bedstat BACKEND_R = "r" BACKEND_GTARS = "gtars" -BACKEND_GTARS_PY = "gtars-py" DEFAULT_PRECISION = 3 # bedbuncher From 79ecfbab075a09c60c5d1a28899ea0c8b7677421 Mon Sep 17 00:00:00 2001 From: Sam Park Date: Mon, 6 Apr 2026 03:10:48 -0400 Subject: [PATCH 09/10] Apply black formatting to backends/__init__.py and ref_utils.py Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedstat/backends/__init__.py | 4 +--- bedboss/bedstat/ref_utils.py | 4 +--- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/bedboss/bedstat/backends/__init__.py b/bedboss/bedstat/backends/__init__.py index 0d28cdbf..f0c724a5 100644 --- a/bedboss/bedstat/backends/__init__.py +++ b/bedboss/bedstat/backends/__init__.py @@ -28,9 +28,7 @@ def create_backend(name: str, **kwargs) -> StatBackend: elif name == BACKEND_GTARS: return GtarsStatBackend(**kwargs) else: - raise ValueError( - f"Unknown analysis backend: {name!r}. Use 'r' or 'gtars'." - ) + raise ValueError(f"Unknown analysis backend: {name!r}. Use 'r' or 'gtars'.") def build_backend(name: str) -> StatBackend: diff --git a/bedboss/bedstat/ref_utils.py b/bedboss/bedstat/ref_utils.py index 28ac1cdf..5a3f6f85 100644 --- a/bedboss/bedstat/ref_utils.py +++ b/bedboss/bedstat/ref_utils.py @@ -73,9 +73,7 @@ def get_fab_path(genome: str, rfg_config: str = None) -> Union[str, None]: _LOGGER.info(f"Created .fab: {fab_path}") return fab_path else: - _LOGGER.warning( - f"Failed to compile .fab: {result.stderr.strip()}" - ) + _LOGGER.warning(f"Failed to compile .fab: {result.stderr.strip()}") return None except FileNotFoundError: _LOGGER.warning("gtars CLI not found — cannot compile .fab") From 49f21c6234652e2d1d66f40d089647b7390dc230 Mon Sep 17 00:00:00 2001 From: Sam Park Date: Mon, 6 Apr 2026 17:44:52 -0400 Subject: [PATCH 10/10] Fix proportional region distribution bins using chrom_sizes compress_region_distribution now takes chrom_sizes to compute correct per-chromosome array lengths (longest chr = n_bins, shorter chrs proportionally fewer). Without this, arrays were sized by max observed rid which varied by region coverage. Co-Authored-By: Claude Opus 4.6 (1M context) --- bedboss/bedstat/backends/gtars_backend.py | 12 ++++++- bedboss/bedstat/compress_distributions.py | 44 ++++++++++++++++++----- 2 files changed, 47 insertions(+), 9 deletions(-) diff --git a/bedboss/bedstat/backends/gtars_backend.py b/bedboss/bedstat/backends/gtars_backend.py index 2a11ee83..20a2691d 100644 --- a/bedboss/bedstat/backends/gtars_backend.py +++ b/bedboss/bedstat/backends/gtars_backend.py @@ -244,8 +244,18 @@ def compute( else: data["gc_content"] = None + # Read chrom_sizes into a dict for proportional region distribution bins + chrom_sizes_dict = None + if chrom_sizes and os.path.exists(chrom_sizes): + chrom_sizes_dict = {} + with open(chrom_sizes) as f: + for line in f: + parts = line.strip().split("\t") + if len(parts) >= 2: + chrom_sizes_dict[parts[0]] = int(parts[1]) + # Compress distributions for DB storage - compress_distributions(gtars_output) + compress_distributions(gtars_output, chrom_sizes=chrom_sizes_dict) # Store entire augmented gtars JSON as distributions blob data["distributions"] = gtars_output diff --git a/bedboss/bedstat/compress_distributions.py b/bedboss/bedstat/compress_distributions.py index b014155c..ecfb5e64 100644 --- a/bedboss/bedstat/compress_distributions.py +++ b/bedboss/bedstat/compress_distributions.py @@ -165,21 +165,40 @@ def compress_tss_histogram( } -def compress_region_distribution(raw: dict) -> Optional[dict]: +def compress_region_distribution( + raw: dict, + chrom_sizes: Optional[dict] = None, + n_bins: int = 250, +) -> Optional[dict]: """Compress per-chromosome region distribution to dense count arrays. Input: gtars format {"chr1": [{"start": ..., "end": ..., "rid": ...}, ...], ...} Output: {"chr1": [count_at_rid_0, count_at_rid_1, ...], ...} - The array index is the rid (bin index) used by the UI's faceted chart. + When chrom_sizes is provided, each chromosome's array length is + proportional to its size (matching the WASM regionDistribution behavior). + The longest chromosome gets n_bins bins; shorter ones get fewer. + + :param raw: gtars region_distribution grouped by chromosome + :param chrom_sizes: {chr: length} for proportional bin counts + :param n_bins: max bins for the longest chromosome (default 250) """ if not raw: return None + # Compute proportional bin count per chromosome + bins_per_chr = {} + if chrom_sizes: + max_len = max(chrom_sizes.values()) + if max_len > 0: + for chrom, length in chrom_sizes.items(): + bins_per_chr[chrom] = max(1, round(length / max_len * n_bins)) + result = {} for chrom, regions in raw.items(): if not regions: - result[chrom] = [] + n = bins_per_chr.get(chrom, 0) + result[chrom] = [0] * n continue rids = np.array( [r.get("rid", 0) if isinstance(r, dict) else 0 for r in regions], @@ -189,20 +208,27 @@ def compress_region_distribution(raw: dict) -> Optional[dict]: [r.get("n", 1) if isinstance(r, dict) else 1 for r in regions], dtype=np.int32, ) - bins = np.zeros(rids.max() + 1, dtype=np.int64) - np.add.at(bins, rids, counts_arr) + # Array length: proportional to chrom size, or fallback to max_rid + 1 + n = bins_per_chr.get(chrom, int(rids.max()) + 1) + bins = np.zeros(n, dtype=np.int64) + valid = rids < n + np.add.at(bins, rids[valid], counts_arr[valid]) result[chrom] = bins.tolist() return result -def compress_distributions(gtars_output: dict) -> dict: +def compress_distributions( + gtars_output: dict, chrom_sizes: Optional[dict] = None +) -> dict: """Compress all distributions in a gtars genomicdist output. Modifies gtars_output["distributions"] in place, replacing raw arrays with compressed formats. - Returns the modified gtars_output. + :param gtars_output: raw gtars genomicdist JSON output + :param chrom_sizes: {chr: length} for proportional region distribution bins + :return: the modified gtars_output """ dists = gtars_output.get("distributions", {}) @@ -236,7 +262,9 @@ def compress_distributions(gtars_output: dict) -> dict: grouped.setdefault(chrom, []).append(entry) rd = grouped if isinstance(rd, dict): - dists["region_distribution"] = compress_region_distribution(rd) + dists["region_distribution"] = compress_region_distribution( + rd, chrom_sizes=chrom_sizes + ) # chromosome_stats: unchanged (already compact)