Metadata-Version: 2.5
Name: pydreg
Version: 0.3.1.post2
Summary: End-to-end dREG peak calling from PRO-seq/GRO-seq bigWig data
Project-URL: Homepage, https://github.com/adamyhe/pydreg
Project-URL: Repository, https://github.com/adamyhe/pydreg
Project-URL: Issues, https://github.com/adamyhe/pydreg/issues
Project-URL: Paper, https://doi.org/10.64898/2026.09.06.745329
License-Expression: GPL-3.0-only
License-File: LICENSE
Keywords: bigwig,bioinformatics,chro-seq,genomics,gro-seq,peak calling,pro-seq,svm,transcriptional regulation
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: GNU General Public License v3 (GPLv3)
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.11
Requires-Dist: huggingface-hub>=0.24
Requires-Dist: numba>=0.64.0
Requires-Dist: numpy>=1.26
Requires-Dist: pandas>=2.0
Requires-Dist: pybigtools>=0.3
Requires-Dist: pysam>=0.22
Requires-Dist: safetensors>=0.4
Requires-Dist: scikit-learn>=1.4
Requires-Dist: scipy>=1.16
Requires-Dist: statsmodels>=0.14
Requires-Dist: threadpoolctl>=3.6.0
Requires-Dist: tqdm>=4.68.4
Requires-Dist: zstandard>=0.22
Provides-Extra: gpu
Requires-Dist: cupy-cuda12x[ctk]>=13.0; (sys_platform == 'linux' and platform_machine in 'x86_64 aarch64') and extra == 'gpu'
Provides-Extra: mlx
Requires-Dist: mlx>=0.20; (sys_platform == 'darwin' and platform_machine == 'arm64') and extra == 'mlx'
Description-Content-Type: text/markdown

# pydreg

[![PyPI](https://img.shields.io/pypi/v/pydreg)](https://pypi.org/project/pydreg/)
[![Tests](https://github.com/adamyhe/pydreg/actions/workflows/ci.yml/badge.svg)](https://github.com/adamyhe/pydreg/actions/workflows/ci.yml)
[![Weights](https://img.shields.io/badge/%F0%9F%A4%97-Weights-yellow)](https://huggingface.co/adamyhe/pydreg)
[![PyPI Downloads](https://static.pepy.tech/personalized-badge/pydreg?period=total&units=INTERNATIONAL_SYSTEM&left_color=BLACK&right_color=GREEN&left_text=downloads)](https://pepy.tech/projects/pydreg)

An accelerated Python port of [dREG](https://github.com/Danko-Lab/dREG) (Danko Lab) — detects active transcriptional regulatory elements (promoters and enhancers) from PRO-seq/GRO-seq nascent-transcription data.

Given a pair of strand-specific bigWig files, pydreg scores every informative genomic position with a pretrained SVR model, then calls significant peaks with FDR control, mirroring the original R package's recommended `run_dREG.R` pipeline end to end.

## Performance

pydreg is substantially faster and less memory demanding than dREG. On an NVIDIA Titan Xp with 16 cores, across 12 real PRO-seq/GRO-seq/ChRO-seq libraries we see a median 4.51x walltime speedup (range 4.40x-4.71x) and a median 5.38x reduction in peak RSS (range 3.78x-6.75x):

<p>
  <img src="https://raw.githubusercontent.com/adamyhe/pydreg/main/figures/plots/walltime.svg" alt="dREG versus pydreg walltime" width="45%">
  <img src="https://raw.githubusercontent.com/adamyhe/pydreg/main/figures/plots/memory.svg" alt="dREG versus pydreg peak RSS" width="45%">
</p>

Peak calls from pydreg have ~0.999 Jaccard similarity with official dREG peak calls.

## Installation

```bash
pip install pydreg[gpu]
pydreg --help
```

Or with [uv](https://docs.astral.sh/uv/) — `uv tool install` if you only need the `pydreg` CLI (isolated environment, nothing else to manage):

```bash
uv tool install pydreg[gpu]
pydreg --help
```

If you want the Python API (`from pydreg import pipeline`, see below) available in your own project instead, use `uv add pydreg[gpu]` there, or `uv pip install pydreg[gpu]` into an already-active environment.

`[gpu]` installs [CuPy](https://cupy.dev/) and enables GPU-accelerated scoring on Linux with an NVIDIA GPU (auto-selected whenever one's detected). On macOS with Apple Silicon, install `[mlx]` instead (`pip install pydreg[mlx]`) for GPU-accelerated scoring via [MLX](https://github.com/ml-explore/mlx) on the machine's Metal GPU (also auto-selected whenever usable). Neither is required, but CPU-only scoring is **much** slower and so is not generally recommended.

**GPU requirement (NVIDIA/`cupy`):** compute capability ≥3.0 (essentially any CUDA-capable NVIDIA GPU from the last decade-plus, including older Pascal-class cards). Scoring runs pydreg's own RBF kernel implementation directly on a CuPy device array (`pydreg.backend._build_cupy_predict_fn`) and enables broad CUDA compatibility and fp32 support. v0.1.x used cuML, which came with significant GPU restrictions (no Pascal support) and no fp32 support — see `docs/OPTIMIZATION.md` for the full writeup, including why that library was dropped in favor of a handrolled `cupy` SVM implementation.

**GPU requirement (Apple Silicon/`mlx`):** any Apple Silicon Mac (M1 or later) with a working Metal GPU. Scoring runs the same RBF kernel formula directly on an MLX device array (`pydreg.backend._build_mlx_predict_fn`); see `docs/OPTIMIZATION.md` for the full writeup and real-hardware benchmark numbers.

Pretrained model weights (an RBF-kernel SVR scorer and a small random-forest peak-splitter) are downloaded automatically from [`adamyhe/pydreg`](https://huggingface.co/adamyhe/pydreg) on Hugging Face the first time they're needed, and cached locally by `huggingface_hub`.

If you already have compatible exported model files, pass them with `--svr-model /path/to/svm.model.safetensors.zst` and `--rf-model /path/to/rf.model.safetensors.zst`, or use `pipeline.run(..., svr_model_path=..., rf_model_path=...)` from Python.

## Usage

### CLI

```bash
pydreg plus.bw minus.bw out_prefix --verbose --cores 8
```

- `plus.bw`/`minus.bw`: strand-specific bigWig files (3′-mapped, point-mode, unnormalized **read counts** — the same input format the original dREG expects). See [proseq2.0](https://github.com/Danko-Lab/proseq2.0/) for the Danko lab's pipeline. `minus.bw` may be positive- or negative-signed — `pydreg` takes the absolute value of both strands during feature extraction (matching the original C implementation), so sign convention doesn't affect scoring.
- `out_prefix`: prefix for all output files (see below).

Options:

**Primary options.** These are the main options that most people will want to set based on their system.

| flag                                  | default          | meaning                                                                                                                                                                                                                                                                                                                                       |
| ------------------------------------- | ---------------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `--backend {auto,cupy,mlx,sklearn,numpy}` | `auto`       | Scoring backend. `auto` uses `cupy` when a usable CUDA device is detected, `mlx` when a usable Apple Silicon (Metal) GPU is detected, otherwise pure NumPy. scikit-learn is selectable explicitly but is not auto-selected (see `docs/OPTIMIZATION.md` for why). An explicit choice raises if that backend isn't actually usable, rather than silently falling back. |
| `-p`, `--cores N`                     | `1`              | Cores pydreg is allowed to use, applied consistently across the pipeline: worker processes for the final peak-calling stage (embarrassingly parallel across broad candidate peaks) *and* numba's thread count for the parallelized feature-extraction/informative-position-scanning kernels — one number, not two independently-tunable ones. **Set this to the max number of cores you can spare on your machine** — the default of `1` is just a safe, non-surprising starting point, not a recommendation.                                                                     |
| `--svr-model PATH`                    | pretrained       | Local SVR scorer model file (`.safetensors` or `.safetensors.zst`) instead of downloading/loading the default Hugging Face weight.                                                                                                                                                                                                           |
| `--rf-model PATH`                     | pretrained       | Local random-forest peak-splitter model file (`.safetensors` or `.safetensors.zst`) instead of downloading/loading the default Hugging Face weight.                                                                                                                                                                                           |
| `--no-progress`                       | off              | Disable tqdm progress bars (auto-hidden anyway when stdout isn't a terminal).                                                                                                                                                                                                                                                                 |
| `-v`, `--verbose`                     | off              | Log progress at INFO level.                                                                                                                                                                                                                                                                                                                   |

**Additional tuning arguments**. These change how pydreg handles batching. These are safe to tune, but the default settings used will generally work well across most systems. In particular, there's no real performance gain from increasing chunk size, at least in my testing.

| flag                                  | default          | meaning      |
| ------------------------------------- | ---------------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `--peak-calling-block-width N`        | `100`            | Candidate broad peaks handed to each peak-calling worker per task; smaller blocks improve load balancing on uneven peak sizes. Tune alongside `--cores`.                                                                                                                                                                         |
| `--query-chunk N`                     | backend-specific | Positions scored per batch; defaults to a size tuned per backend (`pydreg.backend.DEFAULT_QUERY_CHUNK`). Pure batching — does not change scores.                                                                                                                                                                                              |
| `--cupy-sv-chunk N`                   | `16384`          | Support vectors (of 605,187) evaluated per GPU kernel/GEMM call for the `cupy` backend specifically. The main lever for trading GPU memory for fewer, larger, better-amortized kernel launches — measured memory-bandwidth-bound (not launch-overhead-bound) on both a TITAN Xp and an A100, so growing this past the default buys no throughput, only more VRAM (see `docs/OPTIMIZATION.md`). Pure batching — does not change scores. |
| `--mlx-sv-chunk N`                    | `16384`          | The same lever as `--cupy-sv-chunk`, but for the `mlx` backend specifically. Unlike `cupy`, throughput does measurably degrade well past the default (confirmed on a real Apple M4), and an unchunked call can crash outright past Metal's hard per-buffer allocation limit — see `docs/OPTIMIZATION.md`. Pure batching at/below the default — does not change scores. |

**Arguments that will change results.** `--pv-adjust`/`--pv-threshold` are genuine statistical choices. The rest reproduce specific constants/tolerances from legacy dREG's own hardcoded behavior; pydreg's >0.999 Jaccard agreement with real dREG was measured at their defaults, so moving off those defaults is unexplored territory — nothing downstream validates the result against R at other values.

| flag                         | default | meaning                                                                                                                                                                                                         |
| ---------------------------- | ------- | --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |
| `--pv-adjust METHOD`         | `fdr`   | Multiple-testing correction method (any `statsmodels.stats.multitest.multipletests` method name).                                                                                                               |
| `--pv-threshold P`           | `0.05`  | Significance threshold applied after correction.                                                                                                                                                                |
| `--smoothwidth N`            | `4`     | Smoothing window used during peak-splitting; matches legacy dREG's own hardcoded `smoothwidth=4` in `find_rf_peaks`.                                                                                            |
| `--pmv-laplace-cdf-maxpts N` | `25000` | Max integration points for the per-summit p-value's quasi-Monte-Carlo integral; matches R's `mvtnorm::pmvnorm()`/`GenzBretz()` default. Only lower this if you want to trade fidelity with R for further speed. |
| `--pmv-laplace-cdf-eps EPS`  | `0.001` | Absolute/relative tolerance for the same integral; also matches R's default. Only lower this (i.e. tighten precision) if you specifically want to exceed R's own reference precision, at real speed cost.       |
| `--pmv-laplace-tail-tol TOL` | `1e-6`  | Stops the per-summit p-value's z-grid integration early once the remaining tail's provable worst-case error falls below `TOL`, trading a bounded amount of fidelity for speed — unlike the two rows above, this changes *how much* of the integral is evaluated, not just how precisely. The `1e-6` default was validated on a full real production run against real dREG with no measurable fidelity cost (same peaks called, indistinguishable Jaccard agreement — see `docs/OPTIMIZATION.md`). Pass `0.0` for exact, bit-for-bit-with-R behavior. |
| `--pmv-laplace-exact`        | off     | Shorthand for `--pmv-laplace-tail-tol 0.0` (exact, slower). Ignored if `--pmv-laplace-tail-tol` is also given.                                                                                                    |

### Python API

```python
from pydreg import pipeline

result = pipeline.run(
    "plus.bw",
    "minus.bw",
    "out_prefix",
    backend_name=None,
    svr_model_path=None,
    rf_model_path=None,
    cores=8
)
# result: {"dense_infp": ..., "raw_peak": ..., "peak_bed": ..., "min_score": ...}
```

Pass `write_outputs=False` to get the result dict back without writing files, if you just want to work with the DataFrames directly.

### Output files

Given `out_prefix`, pydreg writes:

| file                                                  | contents                                                      |
| ----------------------------------------------------- | ------------------------------------------------------------- |
| `{out_prefix}.dREG.infp.bw`                           | Every informative position and its raw dREG score.            |
| `{out_prefix}.dREG.raw.peak.bed.gz` (+`.tbi`)         | All candidate peaks before FDR filtering.                     |
| `{out_prefix}.dREG.peak.full.bed.gz` (+`.tbi`)        | Significant peaks: chrom, start, end, score, p-value, center. |
| `{out_prefix}.dREG.peak.score.bed.gz`/`.bw` (+`.tbi`) | Significant peaks' scores only.                               |
| `{out_prefix}.dREG.peak.prob.bed.gz`/`.bw` (+`.tbi`)  | Significant peaks' `1 - p-value`.                             |

`.bed.gz` files are bgzipped and tabix-indexed; `.bw` files are standard bigWig tracks.

## How it works

1. **Informative-position scan** — tiles the genome looking for positions with any transcriptional signal on either strand, to avoid scoring silent regions.
2. **Feature extraction** — for each informative position, bins nearby read counts into multiple nested window sizes ("zoom levels") per strand, producing a fixed-length feature vector.
3. **Scoring** — an RBF-kernel SVR (605,187 support vectors, trained on the original dREG data) maps each feature vector to a dREG score in ~[0, 1].
4. **Peak calling** — merges scored positions into broad candidate regions, refines local maxima with a small random-forest model to decide where to split adjacent peaks, computes a per-peak p-value, and applies FDR control to select significant peaks.

See `docs/METHODS.md` for a plain-language walkthrough of each stage, and `docs/OPTIMIZATION.md` for the performance design choices layered on top without changing any of the above. `docs/PLANNING.md`/`docs/PERF_LOG.md` are the underlying comprehensive design/research records (full algorithmic spec, every upstream R quirk and why it's kept, every benchmark) for anyone going deeper.

This is validated directly against the original: on real test data, pydreg's called peaks agree with real dREG's at a >0.999 Jaccard index.

## Caveats

- `minus.bw`'s sign doesn't matter: both informative-position detection and feature extraction take the absolute value of the minus-strand signal (the latter matching the original C's `bigwig_readi(..., abs=1, ...)` read call, which strips sign from both strands before any binning) — see `docs/PLANNING.md` for the sourced trace.
- Chromosome/contig selection follows the original plus-strand-driven scan: contigs present only in `plus.bw` can still be scored, with the missing `minus.bw` contig treated as all-zero signal during feature extraction; contigs present only in `minus.bw` are not discovered by the initial informative-position scan.
- A handful of upstream R bugs/quirks are faithfully replicated rather than fixed, because the pretrained model's expected behavior was produced by that exact code (e.g. a `mean()`-argument-binding bug in the p-value calculation, an off-by-one in broad-peak merging that drops the last group per chromosome, and others) — see `docs/PLANNING.md` for the full list and reasoning.
- Peak-calling p-values have small inherent run-to-run noise (the per-summit p-value's underlying quasi-Monte-Carlo integral is unseeded, matching the original R implementation's `mvtnorm::pmvnorm`, which is also unseeded) — this doesn't affect which peaks are called significant in practice, and is reflected in the >0.999 (not exactly 1.0) Jaccard index above.

## Contributing

See [`CONTRIBUTING.md`](CONTRIBUTING.md) for development setup, running tests, and what to read before making algorithmic or performance changes.

## License

GPL-3.0 (matching the original dREG R package, which is GPL-3-licensed).

## Citation

If you use this package, please cite the pydreg preprint:

> He, A. Y., & Danko, C. G. (2026). pydreg: a fast Python package for identifying active cis-regulatory elements from nascent transcription. *bioRxiv*. https://doi.org/10.64898/2026.09.06.745329

Since pydreg is a port that reuses dREG's pretrained models and peak-calling procedure, please also cite the original dREG papers:

> Danko, C. G., Hyland, S. L., Core, L. J., Martins, A. L., Waters, C. T., Lee, H. W., Cheung, V. G., Kraus, W. L., Lis, J. T., & Siepel, A. (2015). Identification of active transcriptional regulatory elements from GRO-seq data. *Nature Methods*, 12(5), 433-438. https://doi.org/10.1038/nmeth.3329

> Wang, Z., Chu, T., Choate, L. A., & Danko, C. G. (2018). Identification of regulatory elements from nascent transcription using dREG. *Genome Research*, 29, 293–303. https://doi.org/10.1101/gr.238279.118

Please also cite the version number of this port to improve reproducibility.
