Scenic regulons

Version Compatibility

Reference examples tested with: pySCENIC 0.12+, ctxcore 0.2+, arboreto 0.1.6+, scanpy 1.10+, loompy 3.0+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

The motif-DB machinery lives in ctxcore; a ctxcore/feather-format version mismatch is the most common silent failure. pySCENIC is most reliable on a dedicated Python 3.10 environment.

SCENIC Regulons

"Identify transcription factor regulons and score TF activity from my scRNA-seq data" -> Run the pySCENIC three-step pipeline: infer TF-target co-expression with GRNBoost2, prune to direct targets by cis-regulatory motif enrichment with cisTarget, then score per-cell regulon activity with AUCell.

• CLI: pyscenic grn -> pyscenic ctx -> pyscenic aucell
• Python: arboreto_with_multiprocessing.py for the GRN step (avoids the dask breakage)

The Single Most Important Modern Insight -- Motif Pruning Is What Converts Co-expression into Directed Regulation

Step 1 (GRNBoost2) produces undirected co-expression only -- it is no better than WGCNA and inherits all of co-expression's confounding (indirect edges, batch, cell-cycle). The entire conceptual payload of SCENIC is Step 2 (cisTarget): for each module it asks whether the candidate TF's binding motif is significantly enriched (NES >= 3.0) in the cis-regulatory space of the module's targets, and keeps only the targets in the motif's leading edge. This (a) imposes a mechanistic prior -- the TF can physically bind near its retained targets, (b) breaks the symmetry of co-expression into a TF -> target direction, and (c) discards indirect targets. A "regulon" is by definition only the post-cisTarget TF plus its direct targets. Modules that were never pruned are co-expression modules, and calling them regulons misuses the word.

The second non-obvious consequence is AUCell: regulon activity is not TF expression. AUCell ranks genes within each cell and computes the area under the recovery curve for the regulon's gene set, so activity can be high even when the TF's own mRNA is dropout-zero (TF transcripts are sparse). Showing TF expression in place of regulon AUC -- or "validating" activity by its correlation with TF expression -- misses the method's point and is circular. SCENIC regulons remain motif-supported co-expression: a strong, directed hypothesis worth a knockdown, not proof of causal regulation.

Pipeline Taxonomy

Step	Tool	Produces	Key parameter	Watch out for
1. GRN	GRNBoost2 (or GENIE3)	TF-target co-expression adjacencies	--seed, --num_workers	stochastic; not reproducible without a fixed seed
2. Prune	cisTarget (ctxcore)	regulons (direct targets)	--nes_threshold 3.0, --rank_threshold 5000	feather DB + motif2TF version must match
3. Score	AUCell	per-cell regulon activity (AUC)	--auc_threshold 0.05	this is the top-fraction, NOT the binarization cut

Decision Tree by Scenario

Scenario	Recommended	Why
scRNA-seq, want TF regulons + per-cell activity	pySCENIC grn/ctx/aucell	the canonical workflow
GRN step hangs / KilledWorker	arboreto_with_multiprocessing.py	arboreto's dask backend breaks on newer dask
Need reproducible regulons	run GRN 10-100x, keep links recurring >80%	GRNBoost2/GENIE3 are stochastic
Which regulons mark a cell type	Regulon Specificity Score (RSS)	JSD-based specificity, not just magnitude
Paired scRNA + scATAC available	-> multiomics-grn (SCENIC+)	accessibility defines enhancers; eRegulons add the region layer
Bulk RNA-seq / want protein activity	-> grn-inference (ARACNe + VIPER)	SCENIC is single-cell; VIPER reads TF activity from bulk
Compare activity across conditions	run SCENIC once on the integrated object	raw AUC is population-relative; batch survives into regulons

Required Databases

cisTarget needs three matched resources: ranking database(s), motif-to-TF annotations, and the TF list -- all the same species/assembly/symbol namespace. Download from resources.aertslab.org/cistarget/.

# Human hg38 gene-based rankings (~1.5 GB each). Run ctx with BOTH search-space DBs
# (500bp+100bp around TSS, and TSS +/-10kb) so the leading-edge logic pools them.
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc9nr/gene_based/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl
# The ranking-DB version (mc9nr / v10) and the motif2tf annotation version MUST match.

Step 1: GRN Inference (use the multiprocessing wrapper)

Goal: Infer TF-target co-expression adjacencies as candidate regulatory modules.

Approach: Run GRNBoost2 via the bundled multiprocessing script (single-node, stable) rather than the dask backend, and fix the seed so the stochastic boosting is reproducible.

# arboreto's dask backend breaks on dask>=2.x (silent hangs, KilledWorker).
# The bundled multiprocessing wrapper is the supported workaround.
python arboreto_with_multiprocessing.py \
    filtered.loom allTFs_hg38.txt \
    --method grnboost2 --output adj.tsv \
    --num_workers 8 --seed 42

Step 2: Prune to Regulons by Motif Enrichment

Goal: Keep only TF-target links whose target genes are enriched for the TF's binding motif -- the step that confers directness and direction.

Approach: Load the ranking databases and motif2TF annotations, build candidate modules from the adjacencies, and run cisTarget pruning; targets surviving motif enrichment (NES >= 3.0) form the regulon.

import glob, pickle, pandas as pd
from pyscenic.utils import modules_from_adjacencies
from pyscenic.prune import prune2df, df2regulons
from ctxcore.rnkdb import FeatherRankingDatabase

adjacencies = pd.read_csv('adj.tsv', sep='\t')
expr = pd.read_csv('expr.csv', index_col=0)            # cells x genes
modules = list(modules_from_adjacencies(adjacencies, expr))

dbs = [FeatherRankingDatabase(f, name=f) for f in glob.glob('*.genes_vs_motifs.rankings.feather')]
# rank_threshold=5000 matches the CLI default (the prune2df Python default is 1500).
df = prune2df(dbs, modules, 'motifs-v9-nr.hgnc-m0.001-o0.0.tbl', rank_threshold=5000)
regulons = df2regulons(df)                              # TF + direct targets only

with open('regulons.pkl', 'wb') as fh:
    pickle.dump(regulons, fh)

CLI equivalent for steps 1-2 (pyscenic grn, then pyscenic ctx adj.tsv DB.feather --annotations_fname motifs.tbl --expression_mtx_fname filtered.loom -o reg.csv). ctx verified defaults: --rank_threshold 5000, --auc_threshold 0.05, --nes_threshold 3.0, --min_genes 20. --mask_dropouts now defaults to False (matching R SCENIC); it changes the TF-target correlation sign that splits activating (+) from repressing (-) regulons, so report the setting used.

Step 3: AUCell Per-Cell Activity

Goal: Score each regulon's activity in every cell, robustly to dropout.

Approach: Rank genes within each cell, integrate the recovery curve over the top fraction (auc_threshold, default 0.05 = top 5%), and emit a cell-by-regulon AUC matrix.

from pyscenic.aucell import aucell

# auc_threshold = top 5% of the ranking integrated for the AUC -- NOT a binarization cut.
auc_mtx = aucell(expr, regulons, auc_threshold=0.05, num_workers=8)
auc_mtx.to_csv('auc_matrix.csv')

Interpretation: Specificity and Binarization

Goal: Surface the regulons that define each cell type and convert activity to on/off states for clustering.

Approach: Use the Regulon Specificity Score (Jensen-Shannon divergence vs an idealized cell-type-specific distribution) for identity regulators, and binarize the AUC distribution (bimodal -> density threshold) for state heatmaps.

from pyscenic.rss import regulon_specificity_scores
from pyscenic.binarization import binarize

cell_types = pd.read_csv('cell_types.csv', index_col=0)['cell_type']
rss = regulon_specificity_scores(auc_mtx, cell_types)     # high RSS = identity regulator
binary_mtx, thresholds = binarize(auc_mtx)                # per-regulon on/off

RSS (rewards specificity) and a per-cluster AUC z-score (rewards magnitude) can disagree; prefer RSS for "which regulon marks this cluster."

Per-Method Failure Modes

Calling unpruned modules "regulons"

Trigger: skipping ctx, or dropping the NES threshold to admit everything. Mechanism: without motif enrichment the output is co-expression, not direct regulation. Symptom: no motif DB/version reported; implausibly large "regulons." Fix: always run cisTarget; report DB + motif2TF versions and the search-space windows.

Dask hang in the GRN step

Trigger: native arboreto on dask>=2.x. Mechanism: scheduler incompatibility. Symptom: silent hang or KilledWorker. Fix: use arboreto_with_multiprocessing.py (single-node, stable).

Species / assembly mismatch

Trigger: mouse genes against an hg38 ranking DB, or HGNC vs MGI symbol mismatch. Mechanism: gene IDs do not map into the database. Symptom: near-empty regulon set. Fix: match expression IDs, ranking DB, and motif2TF to one species/assembly/namespace.

Cross-condition AUC comparison without batch control

Trigger: comparing raw AUC across separately-run SCENIC analyses or strong batches. Mechanism: AUC is relative to the population it was ranked within; batch-driven co-expression can pass motif enrichment by chance. Symptom: a "condition-specific regulator" that tracks the batch. Fix: run SCENIC once on the integrated object; sanity-check condition regulons against batch.

Over-reading _extended or _- regulons

Trigger: using _extended regulons for direct-binding claims, or building a story on (-) repressor activity. Mechanism: _extended adds orthology/similarity-inferred (low-confidence) motif annotations; negative regulons are sparse and weakly enriched. Symptom: direct-regulation claims from low-confidence edges. Fix: default to high-confidence positive regulons; treat _extended/(-) as hypotheses.

Quantitative Thresholds

Threshold	Source	Rationale
NES >= 3.0 (motif enrichment)	Aibar 2017 / iRegulon (Janky 2014)	recovery-curve enrichment cutoff defining a supported motif
auc_threshold = 0.05 (top 5%)	pySCENIC default	fraction of the ranking integrated for the AUC
GRN reruns: keep links recurring >80% of runs	Van de Sande 2020	GRNBoost2/GENIE3 are stochastic; recurrence = high confidence
min_genes = 20 per regulon	pySCENIC default	smaller target sets give unstable AUC
>= a few hundred cells per cell type	practical	rare clusters and doublets inflate spurious regulons

Common Errors

Error / symptom	Cause	Solution
"not a cisTarget Feather database in v1 or v2 format"	ctxcore/DB version mismatch	download current DB; align ctxcore version
empty regulon set	species/assembly or symbol mismatch	match gene IDs to the DB namespace
different regulons each run	unset seed in GRN step	fix --seed; run multiple seeds and intersect
activity != TF expression confuses the reader	conflating regulon AUC with TF mRNA	report AUCell activity; that independence is the point
ctx returns nothing	missing/mismatched --annotations_fname	supply matching motif2TF; check DB is gene-based (not region-based)

References

• Aibar S, et al. 2017. SCENIC: single-cell regulatory network inference and clustering. Nat Methods 14(11):1083-1086.
• Van de Sande B, et al. 2020. A scalable SCENIC workflow for single-cell gene regulatory network analysis. Nat Protoc 15(7):2247-2276.
• Moerman T, et al. 2019. GRNBoost2 and Arboreto. Bioinformatics 35(12):2159-2161.
• Janky R, et al. 2014. iRegulon: cisTarget ranking-and-recovery framework. PLoS Comput Biol 10(7):e1003731.
• Suo S, et al. 2018. Revealing critical regulators of cell identity (Regulon Specificity Score). Cell Rep 25(6):1436-1445.e3.
• Huynh-Thu VA, et al. 2010. GENIE3. PLoS ONE 5(9):e12776.

Related Skills

• multiomics-grn - enhancer-driven eRegulons from paired scRNA+scATAC (SCENIC+)
• grn-inference - bulk GRN inference and VIPER TF protein-activity (the Califano lineage)
• coexpression-networks - undirected co-expression modules (what step 1 produces alone)
• single-cell/clustering - cluster cells before regulon and RSS analysis
• single-cell/preprocessing - QC, doublet removal, and normalization of scRNA-seq inputs
• single-cell/doublet-detection - remove doublets that inflate spurious regulons