diff --git a/bedboss/bedboss.py b/bedboss/bedboss.py index 4f1b1140..4ff77164 100644 --- a/bedboss/bedboss.py +++ b/bedboss/bedboss.py @@ -434,9 +434,8 @@ def insert_pep( # Build the stats backend once for the whole batch. The backend holds # per-backend resources (persistent R service, gtars reference caches, # etc.) that should be reused across files. - stat_backend = ( - build_backend(bbagent.config.config.analysis.backend) if not lite else None - ) + backend_name = bbagent.config.config.analysis.backend + stat_backend = build_backend(backend_name) if not lite else None for i, pep_sample in enumerate(pep.samples): is_processed = skipper.is_processed(pep_sample.sample_name) @@ -520,6 +519,7 @@ def insert_pep( force_overwrite=force_overwrite, annotation=bedset_annotation, lite=lite, + backend=backend_name, ) else: _LOGGER.info( @@ -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 0a9c9ecd..6ae814c4 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"). For 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,8 @@ def run_bedbuncher( "bedsets", ) - if heavy: + gtars_like = backend == "gtars" + if heavy and not gtars_like: _LOGGER.info("Heavy processing is True. Calculating plots...") plot_value = create_plots( bedset=bed_set, @@ -148,7 +151,12 @@ def run_bedbuncher( ) plots = BedSetPlots(region_commonality=FileModel(**plot_value)) else: - _LOGGER.info("Heavy processing is False. Plots won't be calculated") + if gtars_like and heavy: + _LOGGER.info( + f"Heavy processing ignored for {backend} backend (no R plots available)" + ) + else: + _LOGGER.info("Heavy processing is False. Plots won't be calculated") plots = None bbagent.bedset.create( @@ -178,6 +186,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 +233,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 96e7810d..f0c724a5 100644 --- a/bedboss/bedstat/backends/__init__.py +++ b/bedboss/bedstat/backends/__init__.py @@ -1,8 +1,15 @@ 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", "build_backend"] +__all__ = [ + "StatBackend", + "RStatBackend", + "GtarsStatBackend", + "create_backend", + "build_backend", +] def create_backend(name: str, **kwargs) -> StatBackend: @@ -12,14 +19,14 @@ 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 or BACKEND_GTARS) + :param name: Backend name ('r' or 'gtars') :param kwargs: Backend-specific keyword arguments :return: StatBackend instance """ 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'.") @@ -35,7 +42,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 or BACKEND_GTARS) + :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 new file mode 100644 index 00000000..20a2691d --- /dev/null +++ b/bedboss/bedstat/backends/gtars_backend.py @@ -0,0 +1,269 @@ +"""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. + +Uses .fab binary FASTA format when available for optimal GC content +performance (zero-copy mmap). Falls back to plain FASTA otherwise. +""" + +import json +import logging +import os +import statistics +from pathlib import Path +from typing import Union + +import pypiper + +from bedboss.bedstat.backends.base import StatBackend +from bedboss.bedstat.compress_distributions import ( + compress_distributions, + compress_to_kde, +) +from bedboss.bedstat.ref_utils import ( + get_chrom_sizes_path, + get_fab_path, + get_fasta_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 CLI backend — pure subprocess, no Python bindings.""" + + 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: + 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: + 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." + ) + + # 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) + + # Used to stop pipeline if bedstat is used independently + stop_pipeline = not pm + + 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]) + if fasta_path: + cmd_parts.extend(["--fasta", fasta_path, "--ignore-unk-chroms"]) + 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") + + # 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: + 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: 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 + 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 + + # 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, chrom_sizes=chrom_sizes_dict) + + # 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/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/compress_distributions.py b/bedboss/bedstat/compress_distributions.py new file mode 100644 index 00000000..ecfb5e64 --- /dev/null +++ b/bedboss/bedstat/compress_distributions.py @@ -0,0 +1,272 @@ +"""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, + 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, ...], ...} + + 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: + 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], + dtype=np.int32, + ) + counts_arr = np.array( + [r.get("n", 1) if isinstance(r, dict) else 1 for r in regions], + dtype=np.int32, + ) + # 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, 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. + + :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", {}) + + # 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, chrom_sizes=chrom_sizes + ) + + # 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..5a3f6f85 --- /dev/null +++ b/bedboss/bedstat/ref_utils.py @@ -0,0 +1,270 @@ +"""Reference file resolution for gtars genomicdist. + +Handles auto-fetching and pre-compilation of GTF annotations, chrom.sizes, +open signal matrices, and FASTA files 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_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_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. + + 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 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( + genome_name=genome, + asset_name="fasta", + tag_name="default", + seek_key="chrom_sizes", + ) + except (UndefinedAliasError, RefgenconfError): + 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 None + + +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 e9ee7c25..d943e062 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, @@ -395,7 +402,8 @@ 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' or 'gtars'). Default: 'r'." + "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 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"