PyPeakRankR: Reproducible Peak Feature Extraction

• PyPeakRankR is an open-source Python package for reproducible extraction of quantitative features from chromatin accessibility peaks.
• It generates a standardized peak-feature matrix including BigWig signal summaries, GC content, conservation scores, and higher-order moments.
• The tool separates deterministic feature extraction from downstream ranking, offering both CLI and Python API support and cross-assembly liftOver functionality.

Searching arXiv for the PyPeakRankR paper and directly related work to ground the article in current preprints. PyPeakRankR is an open-source Python package for reproducible peak-level feature extraction in regulatory element ranking. It was introduced to address the lack of a standardized tool for assembling the diverse quantitative features needed to prioritize chromatin accessibility peaks for functional validation, particularly in workflows based on ATAC-seq. The package extracts BigWig signal summaries, GC content, PhyloP conservation scores, distribution moments, and cell-type specificity rankings into a single reproducible peak-by-feature matrix stored as a tab-separated values (TSV) file; it separates deterministic feature extraction from downstream ranking, provides both a command-line interface and a matching Python API, supports cross-assembly scoring via liftOver, and runs in minutes on thousands of peaks. In the reported validation contexts, its predecessor PeakRankR ranked among the top 3 of 16 methods in the Brain Initiative Cell Census Network (BICCN) community challenge, and PyPeakRankR was used within the Cross-species Enhancer Ranking Pipeline (CERP) in a basal ganglia study to identify enhancer-AAV tools achieving greater than 70% on-target specificity across cell types (Somasundaram et al., 16 Jun 2026).

1. Scope and problem setting

High-throughput chromatin accessibility assays such as bulk or single-cell ATAC-seq generate tens of thousands of candidate regulatory elements (“peaks”). Downstream validation, including enhancer selection for reporter assays and design of cell-type-specific AAV tools, requires consistent and reproducible measures of signal intensity, sequence content, conservation, shape, and specificity. The motivating observation behind PyPeakRankR is that many laboratories assemble these quantities through ad hoc scripts, producing inconsistent definitions, poor portability, and difficulty in benchmarking (Somasundaram et al., 16 Jun 2026).

PyPeakRankR was designed as a “table-first” pipeline. It starts from a BED or TSV of peaks and appends columns for each feature without rewriting coordinates. This design is paired with deterministic feature extraction: given the same peak list, BigWig files, FASTA, and conservation tracks, the package returns the same TSV. A central conceptual distinction in the system is the separation between “feature extraction” and “ranking.” Features can be computed once, after which different prioritization strategies or machine-learning procedures can be evaluated on the same upstream matrix. This separation is significant because it constrains upstream variability and makes benchmarking more transparent (Somasundaram et al., 16 Jun 2026).

In typical workflows, PyPeakRankR sits immediately downstream of peak calling, for example after MACS2, and upstream of the ranking or machine-learning model selected for a particular study. This suggests that the package is intended as a standardizing intermediate layer rather than a complete end-to-end inference system (Somasundaram et al., 16 Jun 2026).

2. Feature model and extracted quantities

For each peak, PyPeakRankR can append signal-derived, sequence-derived, conservation-derived, and specificity-derived columns. The extracted feature classes reported for the package are summarized below.

Feature class	Reported quantities
BigWig signal summaries	mean, maximum, sum, user-selectable percentiles
GC content	GC fraction
PhyloP conservation	mean, max, sum, etc.
Signal-distribution moments	skewness, kurtosis, bimodality
Cell-type specificity ranking	ATAC specificity score per peak per group, min-max normalized to $[0,1]$

The BigWig signal summaries include mean signal over the peak interval, maximum signal, sum (total) signal, and user-selectable percentiles such as the 25th and 75th percentiles. Under the hood, PyPeakRankR uses pyBigWig to read per-base coverage from one or more BigWig tracks and computes summary statistics in NumPy. GC content is obtained by extracting the underlying sequence via pyfaidx and computing

$\text{GC fraction} = \frac{G + C}{A + T + G + C}.$

For conservation, the package extracts per-base PhyloP scores from a BigWig track and computes summaries such as mean, max, and sum exactly as for signal (Somasundaram et al., 16 Jun 2026).

A distinctive aspect of PyPeakRankR is the inclusion of higher-order moments of the per-base signal distribution. Let $x_1 \ldots x_n$ denote coverage values across $n$ bases in a peak, with mean $\mu$ and standard deviation $\sigma$. The package computes skewness,

$$s \;=\;\frac{1}{n}\sum_{i=1}^n\left(\frac{x_i-\mu}{\sigma}\right)^3,$$

kurtosis,

$$k \;=\;\frac{1}{n}\sum_{i=1}^n\left(\frac{x_i-\mu}{\sigma}\right)^4,$$

and the bimodality coefficient (Sarle’s coefficient),

$$b \;=\;\frac{s^2 + 1}{k + 3\,\frac{(n-1)^2}{(n-2)(n-3)}}.$$

SciPy provides the implementations of these metrics that PyPeakRankR calls directly. The inclusion of skewness, kurtosis, and bimodality extends peak representation beyond magnitude-based summaries and formalizes “shape” as a quantitative component of regulatory element ranking (Somasundaram et al., 16 Jun 2026).

For multi-group experiments, PyPeakRankR can compute an “ATAC specificity” score per peak per group according to

$$\text{specificity}_{g} = \frac{\text{signal}_{g}}{\operatorname{mean}_{h\neq g}\text{ signal}_{h}},$$

followed by min-max normalization across all peaks to the range $\text{GC fraction} = \frac{G + C}{A + T + G + C}.$0. The paper states that this matches the ranking formula used by CERP and that it was one of the winning metrics in the BICCN community challenge (Somasundaram et al., 16 Jun 2026).

3. Software architecture and interfaces

PyPeakRankR is implemented in pure Python (version $\text{GC fraction} = \frac{G + C}{A + T + G + C}.$1) with the dependencies pandas, NumPy, SciPy, pyBigWig, and pyfaidx. It is distributed with both a command-line interface and a Python API. The dual-interface design is not incidental: the same functionality can be embedded into shell-based pipelines or notebook-centric analysis environments without rewriting the core logic (Somasundaram et al., 16 Jun 2026).

The command-line interface follows a subcommand pattern. init-table creates a TSV from a peak list; add-signal, add-gc, add-phylop, and add-moments append the corresponding feature columns; and rank-by-specificity computes group-specific specificity ranks. Each subcommand reads an input TSV, appends new columns, and writes a new TSV. Because the tool never mutates coordinates, feature extraction steps can be run in any order, and a single step can be rerun if parameters change. This coordinate invariance is a key part of the package’s reproducibility model (Somasundaram et al., 16 Jun 2026).

At the API level, the main functions are exposed in pypeakranker.core, with a minimal workflow based on init_table, add_signal, add_gc, add_phylop, add_moments, and rank_by_specificity. Each function returns a pandas.DataFrame, making it possible to interleave custom feature columns or filtering steps before writing the final TSV. This suggests that PyPeakRankR is designed not only as a fixed pipeline but also as a composable library for analytical workflows that need a stable feature-extraction substrate (Somasundaram et al., 16 Jun 2026).

A recurrent misconception in peak-ranking workflows is that ranking logic and feature extraction must be tightly coupled. PyPeakRankR explicitly rejects that coupling: it computes the feature matrix once and leaves ranking or machine-learning strategy selection to downstream analysis. In the terms used by the paper, the package provides deterministic feature extraction while enabling transparent benchmarking of prioritization strategies on the same upstream data (Somasundaram et al., 16 Jun 2026).

4. Cross-assembly scoring and reproducibility model

When assay data are aligned to one assembly but annotations are available on another, PyPeakRankR can call UCSC liftOver under the hood. The reported usage adds --liftover-chain and --min-match to any of the add-* commands. Internally, the package invokes the UCSC liftOver binary on the input intervals, filters out intervals that fail mapping or do not satisfy the minimum match threshold, and then proceeds with feature extraction on the mapped intervals (Somasundaram et al., 16 Jun 2026).

The package documentation identifies several caveats for this cross-assembly mode. A compatible chain file must be supplied; some peaks may be lost or split during conversion; the software logs how many intervals failed to map; and coordinate accuracy depends on chain completeness, so the liftOver report should be inspected. These caveats delimit the scope of cross-assembly comparability and indicate that liftOver support is operationally integrated but not treated as a lossless transformation (Somasundaram et al., 16 Jun 2026).

Reproducibility is framed as both a software property and an experimental design principle. PyPeakRankR was designed for speed and determinism: extraction of all standard features on approximately 10,000 peaks takes on the order of 2–5 minutes on a modern 8-core workstation, and because each subcommand is pure Python plus vectorized NumPy/SciPy, runs can be parallelized by processing separate BigWig files concurrently. At the same time, the package reports deterministic outputs with no random seeds and no multithreading nondeterminism, such that rerunning the same pipeline produces byte-identical TSVs (Somasundaram et al., 16 Jun 2026).

The significance of this design is explicit in the paper. When different ranking strategies are benchmarked on exactly the same feature matrix, performance differences cannot be attributed to upstream variability. In this sense, PyPeakRankR treats reproducibility not merely as software hygiene but as a condition for valid comparative evaluation (Somasundaram et al., 16 Jun 2026).

5. Validation contexts and empirical use

The paper situates PyPeakRankR within two validation contexts. The first is the BICCN community challenge, described as a benchmark for predicting cell-type-specific enhancers in the mammalian cortex. In that challenge, PeakRankR, the original R implementation with three features—specificity, magnitude, and coverage—placed in the top 3 of 16 submitted methods. PyPeakRankR is presented as a Python re-implementation and extension of that strategy, adding GC content, PhyloP, and distribution moments for contemporary scATAC-seq studies (Somasundaram et al., 16 Jun 2026).

The second context is a cross-species basal ganglia study by Wirthlin et al. in 2026. There, PyPeakRankR was used within the Cross-species Enhancer Ranking Pipeline (CERP) to extract composite features across mouse and macaque basal ganglia cell types. The paper states that composite ranking of these features outperformed simple fold-change approaches and that the resulting enhancer-AAV constructs achieved greater than 70% on-target specificity, with top candidates exceeding 90% specificity in vivo (Somasundaram et al., 16 Jun 2026).

These examples delimit the type of evidence currently attached to PyPeakRankR. The package is not presented as a ranking algorithm with a single canonical scoring rule; rather, it is validated through its role in a reproducible feature-extraction layer that supports successful downstream prioritization. A plausible implication is that its primary scientific contribution lies in standardizing and stabilizing upstream feature generation so that ranking methods can be compared or composed under controlled conditions (Somasundaram et al., 16 Jun 2026).

6. Installation, licensing, and operational constraints

PyPeakRankR installs via pip from GitHub and can also be installed from source by cloning the repository and running pip install . in the python-package directory. The reported prerequisites are Python $\text{GC fraction} = \frac{G + C}{A + T + G + C}.$2 and a UNIX-compatible environment, specifically macOS or Linux, for liftOver support. The package repository is hosted at the Allen Institute GitHub location specified in the paper, and the software is freely available under the MIT license (Somasundaram et al., 16 Jun 2026).

The MIT licensing terms are described as permissive: the code and derived feature matrices may be used, modified, and redistributed provided the original notice is retained. The paper also notes a Zenodo archive, DOI 10.5281/zenodo.15238527, that permanently snapshots the described version. This archival detail is relevant for reproducible computational scholarship because it provides a versioned software target aligned with the reported behavior (Somasundaram et al., 16 Jun 2026).

Operationally, the package depends on pandas, NumPy, SciPy, pyBigWig, and pyfaidx, all of which are installed automatically. Since the software is pure Python and the CLI subcommands wrap pure-Python functions, the implementation emphasizes portability and a direct correspondence between library and command-line usage. This suggests a design intended to minimize divergence between exploratory analysis and production execution environments (Somasundaram et al., 16 Jun 2026).