An accelerated Python port of 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.
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):
Peak calls from pydreg have ~0.999 Jaccard similarity with official dREG peak calls.
pip install pydreg[gpu]
pydreg --helpOr with uv — uv tool install if you only need the pydreg CLI (isolated environment, nothing else to manage):
uv tool install pydreg[gpu]
pydreg --helpIf 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 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 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 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.
pydreg plus.bw minus.bw out_prefix --verboseplus.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 for the Danko lab's pipeline.minus.bwmay be positive- or negative-signed —pydregtakes 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. |
from pydreg import pipeline
result = pipeline.run(
"plus.bw",
"minus.bw",
"out_prefix",
backend_name=None,
svr_model_path=None,
rf_model_path=None,
)
# 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.
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.
- Informative-position scan — tiles the genome looking for positions with any transcriptional signal on either strand, to avoid scoring silent regions.
- 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.
- 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].
- 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.
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'sbigwig_readi(..., abs=1, ...)read call, which strips sign from both strands before any binning) — seedocs/PLANNING.mdfor the sourced trace.- Chromosome/contig selection follows the original plus-strand-driven scan: contigs present only in
plus.bwcan still be scored, with the missingminus.bwcontig treated as all-zero signal during feature extraction; contigs present only inminus.bware 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) — seedocs/PLANNING.mdfor 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.
See CONTRIBUTING.md for development setup, running tests, and what to read before making algorithmic or performance changes.
GPL-3.0 (matching the original dREG R package, which is GPL-3-licensed).
If you use this package, please 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.