rustscenic 0.4.6

Rust and PyO3 implementation of SCENIC-style regulatory-network analysis. Includes GRN, AUCell, topics, cistarget, peak calling, cell QC, enhancer-gene links, and eRegulon assembly. Installs without dask, Java, or CUDA.

rustscenic

A Rust + PyO3 implementation of the practical SCENIC and SCENIC+ compute path: one install, modern Python, low-memory CPU execution, and measured validation across GRN, AUCell, cisTarget, topics, and multiome pipeline stages without Java, dask, or CUDA.

pip install rustscenic

Run the full pipeline:

rustscenic pipeline --rna data.h5ad --tfs tfs.txt --output out/

Five runtime dependencies (numpy, pandas, pyarrow, scipy, anndata). Python 3.10 to 3.13, Linux + macOS (x86_64 + aarch64); Windows x64 is covered by CI and v0.4.6 release wheels. No dask, no Java, no CUDA.

The practical SCENIC+ compute path in one package:

flowchart LR
    rna["RNA data"] --> network["find gene regulation"]
    network --> regulons["regulons"]
    regulons --> activity["score each cell"]

    motifs["motif data"] --> supported["filter by motif support"]
    regulons --> supported

    atac["ATAC data"] --> topics["discover topics"]
    rna --> enhancers["link enhancers to genes"]
    atac --> enhancers

    supported --> programs["build enhancer-linked regulons"]
    enhancers --> programs

Status

Current release: v0.4.6 on PyPI. This micro patch fixes the monolithic region-cistarget path so file-backed region motif rankings are projected to the current run's linked peaks before reading. See CHANGELOG and validation/ for evidence and caveats.

Active limitations are listed under Scope and alternatives, with full detail in site_docs/limitations.md.

Goal

rustscenic is being built as the single-install replacement for the practical SCENIC / SCENIC+ workflow: RNA GRN inference, AUCell regulon activity, motif enrichment, ATAC fragment preprocessing, topic modelling, enhancer-gene linking, and eRegulon assembly in one package.

The project is intentionally not a thin wrapper around the old stack. The target is a simpler architecture that makes regulatory-network analysis easier to install, cheaper to run on CPU, deterministic under a fixed seed, and robust to real atlas conventions such as ENSEMBL var_names, duplicate gene symbols, backed AnnData, and UCSC/Ensembl chromosome mismatches.

What it does

Rust-native replacements for the compute stages plus the glue that scenicplus builds eRegulons from:

Stage	rustscenic	Replaces
Gene-regulatory network inference	rustscenic.grn.infer	arboreto.grnboost2
Per-cell regulon activity scoring	rustscenic.aucell.score	pyscenic.aucell.aucell
Topic modelling on scATAC peaks (Online VB)	rustscenic.topics.fit	pycisTopic (gensim VB)
Topic modelling K ≥ 30 (Mallet-class collapsed Gibbs)	rustscenic.topics.fit_gibbs	pycisTopic (Mallet, Java)
Motif-regulon enrichment	rustscenic.cistarget.enrich	pycistarget AUC kernel
ATAC fragments → cells × peaks matrix	rustscenic.preproc.fragments_to_matrix	pycisTopic fragment loader
Cell QC (TSS enrichment, FRiP, insert size)	rustscenic.preproc.qc	pycisTopic.qc
Enhancer → gene correlation	rustscenic.enhancer.link_peaks_to_genes	scenicplus p2g linking
eRegulon assembly (TF × enhancers × target genes)	rustscenic.eregulon.build_eregulons	scenicplus eRegulon builder
End-to-end pipeline orchestrator	rustscenic.pipeline.run	scenicplus snakemake

Bundled with the wheel: HGNC (1,839 human) and MGI (1,721 mouse) TF lists via rustscenic.data.tfs(species). Motif rankings can be fetched and cached via rustscenic.data.download_motif_rankings. Cellxgene-curated h5ads (ENSEMBL IDs in var_names, gene symbols in var["feature_name"]) are auto-detected so atlas data works without manual patching.

Quick example (PBMC-3k, RNA GRN + AUCell)

import anndata as ad
import rustscenic.grn, rustscenic.aucell

import rustscenic.data

adata = ad.read_h5ad("rna.h5ad")
tfs = rustscenic.data.tfs("hs")  # bundled HGNC list (1,839 TFs)

# 1. GRN inference
grn = rustscenic.grn.infer(adata, tf_names=tfs, n_estimators=5000, seed=777)

# 2. Build top-50-target regulons and score per-cell activity
regulons = [
    (f"{tf}_regulon", grn[grn["TF"] == tf].nlargest(50, "importance")["target"].tolist())
    for tf in grn["TF"].unique()
]
auc = rustscenic.aucell.score(adata, regulons, top_frac=0.05)

Full RNA example script: examples/pbmc3k_end_to_end.py. Runs in about 3 minutes on an 8-core laptop with n_estimators=500. docs/tester-quickstart.md is the collaborator smoke-test path.

Focused external-reader docs are in site_docs/ and are built by the docs workflow with MkDocs.

Measured against the pyscenic / arboreto reference

Same input on both sides. Every row has a log file under validation/. For the public benchmark matrix with dataset, command, hardware, baseline, runtime, memory, parity metric and biological sanity check, see site_docs/benchmarks.md.

Axis	pyscenic / arboreto	rustscenic
Installs on fresh Python 3.10 to 3.13 venv	arboreto: TypeError: Must supply at least one delayed object (dask_expr); pyscenic: ModuleNotFoundError: pkg_resources in current stacks	PyPI wheels and sdist install; core APIs import
AUCell wall-time, Ziegler 2021 atlas (31,602 × 59; measured 2026-04 pre-v0.4.x; refresh deferred to v0.5)	6.81 s (pyscenic)	0.25 s
AUCell wall-time, 10x Multiome (10,290 × 1,457; measured 2026-04 pre-v0.4.x; refresh deferred to v0.5)	18.6 s (pyscenic)	0.21 s
Peak RSS, 4 stages on 100,000 cells × 20,292 genes	> 40 GB (reported)	6.3 GB
Cistarget kernel vs ctxcore.recovery.aucs	reference	Pearson 1.0000, mean abs diff 2.4 × 10⁻⁵
AUCell per-cell Pearson vs pyscenic (Ziegler, 31,602 cells; measured 2026-04 pre-v0.4.x; refresh deferred to v0.5)	reference	0.984 mean, 91.7 % of cells > 0.95
Canonical airway TFs matching literature (Ziegler, n=14)	8 / 14 (pyscenic, unit weights)	8 / 14 - same hits, same 5/14 misses
Bit-identical output under same seed across threaded runs	no (dask non-determinism)	yes
Runtime dependencies	40 +	5

Tool-to-tool variation (same hits, same misses on the same 14 canonical TFs) is smaller than the dataset-inherent noise, consistent with rustscenic being numerically equivalent to pyscenic at the per-cell level.

Community validation reports

External collaborator reports complement the maintainer benchmark set. Each row links the public issue or PR plus committed JSON evidence.

Reporter	Dataset	Stages	Result	Status
@Skycr	Kamath et al. 2022 midbrain dopaminergic neurons, 15,684 cells	GRN + cisTarget	266,805 GRN edges, 9 regulons, 9/9 expected DA-neuron TFs recovered	issue #68, PR #71, validation/community/kamath_da_grn.json
@lmVl12	10x Multiome GEM-X 10k human brain, immune-subsetted 8,215 cells	GRN + AUCell + topics	4,293,902 GRN edges, 1,748 regulons, AUCell/topic outputs non-empty; neural signal improved after immune subsetting	issues #69, #70, PR #74, validation/community/human_brain_10k_v0.4.1.json

Per-stage detail

Numbers are rustscenic's values. The measurement context (dataset, n_cells, version) is in each row. The parity refresh against current upstream stacks (six-dataset sweep) is now planned for v0.5+; see docs/v0.4.x-benchmark-plan.md for the dataset list and success criteria.

GRN - arboreto.grnboost2 replacement

Measurement	Value
Per-edge Spearman vs arboreto (PBMC-3k scanpy, n_estimators=5000, 480,680 shared edges, v0.3.10)	0.611
Within-TF Spearman, mean across 1,274 TFs (same fixture)	0.632 (median 0.649)
Per-edge Spearman vs arboreto (multiome3k, n_estimators=5000, 816 k common edges, 2026-04)	0.58
Per-target TF-ranking Spearman mean	0.57
TRRUST known TF→target edges recovered (PBMC-3k)	17 / 18 (94 %)
Lineage TFs correctly enriched in expected cell types (PBMC-10k)	8 / 8 (SPI1, PAX5, EBF1, TCF7, LEF1, TBX21, CEBPD, IRF8)
Cortex marker TFs present in regulon set (E18 multiome, 4,770 cells, v0.3.10; name-presence, not cell-type enrichment)	9 / 9 (Pax6, Neurod2, Sox2, Ascl1, Tbr1, Neurog2, Fezf2, Eomes, Foxg1)
MITF regulon activity, Tirosh 2016 melanoma - malignant vs TME	3.48×
Wall vs pyscenic on PBMC-3k (n_estimators=5000, seed 777, Apple M5, v0.3.10; pyscenic in sync mode - not apples-to-apples against dask-parallel)	214 s vs 381 s (1.78×)
100k-cell bootstrap, n_estimators=100	17 min / 5.0 GB peak RSS

At high cell counts, GRN target blocking is adaptive by default. Users can force a specific response-block width with rustscenic.grn.infer(..., target_block_size=32) or rustscenic grn --target-block-size 32 when benchmarking cache/RSS behaviour.

Edge rankings disagree with arboreto at fine grain (per-edge Spearman 0.611 on PBMC-3k v0.3.10 / 0.58 on multiome3k 2026-04, top-10k Jaccard 0.20) - expected consequence of independent histogram-GBM quantisation. Coarse biology converges (per-TF Spearman ≈ 0.65, all canonical lineage TFs recovered on both human PBMC and mouse cortex). Downstream AUCell is 0.99 per-cell with pyscenic, so edge-ranking differences do not propagate.

AUCell - pyscenic.aucell replacement

Measurement	Value
Per-cell Pearson vs pyscenic (10x Multiome, 2,588 × 1,457)	0.988 mean, 99.5 % of cells > 0.95
Per-cell Pearson vs pyscenic (Ziegler atlas, 31,602 × 59)	0.984 mean, 91.7 % of cells > 0.95
Per-regulon Pearson (10x Multiome)	0.87 mean, 90.5 % > 0.80
Exact top-regulon-per-cell match (Multiome)	88.4 %
Wall-time, 10k cells × 1,457 regulons	0.21 s (vs 18.6 s pyscenic)
100 k cells × 500 regulons	10 s, 5.6 GB peak RSS

Topics - pycisTopic LDA replacement (Online VB + collapsed Gibbs)

Two algorithms ship side-by-side:

• rustscenic.topics.fit - Online VB LDA, fastest at K ≤ 10.
• rustscenic.topics.fit_gibbs - collapsed Gibbs (Mallet's algorithm class). Add n_threads=N for parallel AD-LDA.

Real PBMC 3k Multiome ATAC, 1,500 cells × 98,319 peaks, K = 30, intrinsic top-10 NPMI on the training corpus:

Tool	Wall	Unique topics (of 30)	Top-10 NPMI mean
rustscenic.topics.fit (Online VB)	104 s	2 / 30 (collapsed)	+0.012
rustscenic.topics.fit_gibbs (serial)	191 s	22 / 30	+0.031
rustscenic.topics.fit_gibbs (8-thread)	84 s	25 / 30	+0.019
Mallet (pycisTopic reference)	n/a	24 / 30	0.196 (extrinsic)

Collapsed Gibbs gives ~11× more distinct topics than Online VB on sparse scATAC at K = 30 and ~2.7× higher intrinsic NPMI; the parallel AD-LDA path adds a 2.56× wall-clock speedup at 8 threads while preserving topic diversity. Mallet's published 0.196 is an extrinsic NPMI (different protocol, not directly comparable in absolute scale). See docs/topic-collapse.md and docs/bench-vs-references.md. Reproduce with python validation/scaling/bench_npmi_head_to_head.py and python validation/scaling/bench_gibbs_parallel.py.

Cistarget - pycistarget AUC kernel replacement

Validated on the aertslab hg38 v10 feather database (5,876 motifs × 27,015 genes):

Measurement	Value
Per-regulon Pearson vs ctxcore.recovery.aucs (58 TRRUST regulons)	1.0000 (all > 0.9999, abs diff 2.4 × 10⁻⁵)
Self-consistency (motif's own top-500 genes → rank #1)	10 / 10
TRRUST at scale (166 TFs ≥ 10 targets): TF-annotated motif ranks #1	19 %
Same benchmark: any TF-motif in top-100	68 to 100 % (rises with regulon size)
Mouse mm10 cross-species (5 TRRUST TFs)	2 / 5 rank #1, 4 / 5 in top-5
100 k-cell workload × 100 regulons	2.6 s, 6.3 GB peak RSS

Bit-identical to ctxcore.recovery.aucs at float32 precision. The 19 % rank-#1 rate is the scaled-out TRRUST-vs-motif-binding benchmark, a property of the gold-standard mismatch, not the implementation.

End-to-end + determinism

Pipeline	Wall	Peak RSS	Stages
Reference (arboreto + pyscenic + tomotopy), 10x Multiome 3k	11.8 min	n/a	4
rustscenic, 10x Multiome 3k	9.1 min	n/a	4
rustscenic, 10x PBMC 3k multiome real-data (v0.3.9, measured 2026-05-02)	7.5 min	3.67 GB	7 (all)
rustscenic, 10x brain E18 5k multiome real-data (v0.3.10, measured 2026-05-04)	13.8 min	4.01 GB	7 (all)
rustscenic, 10x PBMC granulocyte 10k multiome real-data (v0.4.3, measured 2026-05-11)	38.1 min	5.39 GB	7 (all)
rustscenic, 100k synthetic multiome E2E (measured v0.3.10, 2026-04-27)	12.7 min	7.09 GB	7 (all)
rustscenic, 200k synthetic multiome E2E (measured v0.3.10, 2026-04-27)	16.8 min	7.44 GB	7 (all)

Real 10x multiome scaling from 2,767 to 11,620 cells:

• cell count: 4.2x
• wall time: 5.1x, slope about 1.21 over the full span
• peak RSS: 1.47x
• 10k PBMC granulocyte run recovered 10 of 10 canonical TFs by name
• brain E18 5k run recovered 9 of 9 cortex TFs by name

Name-presence checks are not cell-type enrichment tests. Synthetic 100k and 200k runs are scale gates, not biological validation. Full commands, hardware, baseline status, and caveats are in site_docs/benchmarks.md.

Scope and alternatives

rustscenic covers the practical SCENIC / SCENIC+ compute path on CPU. Adjacent tools with different scope:

• GPU, CUDA - flashSCENIC (uses RegDiffusion, a different algorithm from GENIE3 / GRNBoost2, so outputs are not pyscenic-numerical).
• Multiomic enhancer-aware GRN - scenicplus (joint scRNA + scATAC enhancer inference; superset of this scope).
• TF-activity scoring from prebuilt regulons, no GRN inference - decoupler-py with CollecTRI.
• R Bioconductor ecosystem - the original R-SCENIC or Epiregulon.

rustscenic does not bundle the aertslab motif ranking feather databases (300 MB to 35 GB). Users fetch them from resources.aertslab.org and pass the resulting DataFrame to cistarget.enrich.

Current limitations before treating rustscenic as a full SCENIC+ replacement:

• refreshed AUCell timings against current upstream stacks
• region-cisTarget parity checks on real region-ranking databases
• six-dataset benchmark sweep planned for v0.5+
• cell-type enrichment checks for biology claims, not only TF-name recovery
• smoother raw 10x pipeline.run input without caller-side ATAC subsetting

Per-stage CLI

rustscenic grn       --expression data.h5ad --tfs tfs.txt --output grn.parquet
rustscenic aucell    --expression data.h5ad --regulons grn.parquet --output auc.parquet
rustscenic topics    --expression atac.h5ad --output topics --n-topics 30
rustscenic cistarget --rankings motifs.feather --regulons grn.parquet --output enrichment.tsv

Repo layout

• crates/ - Rust workspace: rustscenic-{grn, aucell, topics, preproc, py}
• python/rustscenic/ - Python package, CLI entry point, type stubs
• examples/pbmc3k_end_to_end.py - RNA GRN + AUCell script on real PBMC-3k
• validation/ - reproducible benchmark scripts + measurement reports for every number above, plus VALIDATION_SUMMARY.md
• tests/ - pytest suite (169 Python tests, 1 skipped) + Rust crate tests (57)
• manuscript/ - preprint source
• docs/topic-collapse.md - known algorithmic caveat

License

MIT. Algorithm implementations follow the aertslab Python references - original method credit to Aibar et al. 2017 (SCENIC), Bravo González-Blas et al. 2023 (SCENIC+), Hoffman-Blei-Bach 2010 (Online VB LDA).

Citation and attribution

If you use rustscenic in a paper, report, benchmark, derivative package, or lab workflow, cite the exact release used. GitHub citation metadata is in CITATION.cff. Zenodo concept DOI: 10.5281/zenodo.20246040. Zenodo mints release-specific DOIs from the tagged GitHub releases.

rustscenic was created and is maintained by Ekin Kahraman. See AUTHORS.md and docs/collaboration-and-authorship.md for contribution and authorship expectations.

Contact

File issues at github.com/Ekin-Kahraman/rustscenic/issues. Bug, correctness, and validation-report templates pre-fill the fields we need. If you ran the pipeline on real data and want the result folded into the v0.4.x sweep, see docs/tester-reporting.md. If reporting ARI or related clustering metrics, include the comparator; see docs/evaluation-metrics.md. Coordinated vulnerability disclosure: see SECURITY.md.