Orchestrates gene regulatory network inference from processed single-cell data to regulons and in-silico perturbation, via pySCENIC (RNA-only GRNBoost2 -> cisTarget -> AUCell), SCENIC+ (multiome cisTopic -> pycistarget -> eGRN), and CellOracle perturbation. Use when recognizing that an inferred GRN is UNDIRECTED by default and reporting only the evidence tier delivered (co-expression vs motif-pruned vs enhancer-resolved vs perturbation), matching species/assembly/namespace across the TF-list + cisTarget DB + motif2TF annotation, feeding RAW counts of the cleaned/doublet-free/batch-controlled cells (never imputed/batch-corrected values), running the cisTarget pruning that buys directionality (modules are not regulons without it), or choosing the RNA-only vs multiome path. Hands mechanism to the gene-regulatory-networks component skills; not a re-teach of any single step.
Orchestrates gene regulatory network inference from processed single-cell data to regulons and in-silico perturbation, via pySCENIC (RNA-only GRNBoost2 -> cisTarget -> AUCell), SCENIC+ (multiome cisTopic -> pycistarget -> eGRN), and CellOracle perturbation. Use when recognizing that an inferred GRN is UNDIRECTED by default and reporting only the evidence tier delivered (co-expression vs motif-pruned vs enhancer-resolved vs perturbation), matching species/assembly/namespace across the TF-list + cisTarget DB + motif2TF annotation, feeding RAW counts of the cleaned/doublet-free/batch-controlled cells (never imputed/batch-corrected values), running the cisTarget pruning that buys directionality (modules are not regulons without it), or choosing the RNA-only vs multiome path. Hands mechanism to the gene-regulatory-networks component skills; not a re-teach of any single step.
[{"after_grn_inference":"50-500 regulons detected, known TFs present"},{"after_activity_scoring":"AUCell scores separate known cell types"},{"after_perturbation":"Predicted shifts match known biology"}]
Before using code patterns, verify installed versions match. If versions differ:
Python: pip show <package> then help(module.function) to check signatures
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Note: SCENIC+ is now a Snakemake pipeline (scenicplus init_snakemake); the pre-2024 manual create_SCENICPLUS_object/build_grn object API is deprecated. GRNBoost2's arboreto dask backend is the #1 operational landmine — use the bundled multiprocessing if the dask cluster hangs. The TF-list, cisTarget ranking DB, and motif2TF .tbl must all be the SAME species + assembly + collection vintage. Confirm in-tool before quoting.
Gene Regulatory Network Pipeline
"Infer gene regulatory networks from my single-cell data" -> Orchestrate pySCENIC regulon inference (GRNBoost2, cisTarget, AUCell), CellOracle perturbation simulation, and regulon-based cell type characterization.
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step.
The governing principle
An inferred GRN is an UNDIRECTED association graph by default; directionality is IMPORTED, and only the evidence tier actually delivered should be reported. GRNBoost2 co-expression is tier-1 (undirected); the cisTarget motif-pruning step (Path A step 2) is what buys directionality and discards indirect edges — modules before ctx are NOT regulons, calling them so is a category error. SCENIC+ adds enhancer resolution; perturbation adds causal direction. Do not write tier-6 "master regulator drives X" prose over a tier-1 co-expression result.
Species/assembly/namespace must match across the TF-list, the cisTarget ranking DB, and the motif2TF .tbl — all three. A mismatch (mouse genes in an hg38 DB; feather v1 DB with a v10 motif annotation) yields near-empty regulons. Also commit and report the search-space window (500bp/100bp proximal vs TSS±10kb) — results are not comparable across windows.
Feed RAW counts of the CLEANED, doublet-free, batch-controlled cells — never imputed or batch-corrected values. GRNBoost2 on imputed counts inflates correlations (imputation smooths neighbors into agreement); on batch-corrected values a batch module can pass motif enrichment by chance; doublets create a fake "hybrid regulator." Run SCENIC ONCE on the integrated object.
Validate a regulon with an ORTHOGONAL modality, not the TF's own mRNA. AUCell activity != TF expression; correlating a regulon's activity with its TF's expression is circular (and activity is dropout-robust while the TF mRNA may read zero). GRNBoost2 is stochastic — run multiple seeds and keep recurrent links (Van de Sande 2020).
Made-once commitments
Commitment
Consequence inherited downstream
Species + assembly + gene namespace (HGNC vs MGI; hg38 vs mm10)
TF list, cisTarget DB, motif2TF .tbl must all match, or regulons are near-empty
cisTarget DB vintage + search window (proximal vs TSS±10kb; gene- vs region-based)
Which edges survive ctx pruning; results not comparable across windows/vintages
Input matrix identity (RAW counts, from the cleaned/doublet-free/integrated object)
Every adjacency and regulon; imputed/batch values fabricate edges
RNA-only vs multiome availability
Which path is possible: SCENIC+ REQUIRES paired multiome; RNA-only -> pySCENIC or CellOracle-with-prebuilt-base-GRN
import scanpy as sc
import pandas as pd
from arboreto.algo import grnboost2
adata = sc.read_h5ad('processed.h5ad')
# Extract expression matrix (raw counts recommended for GRNBoost2)
expr_matrix = pd.DataFrame(
adata.raw.X.toarray() ifhasattr(adata.raw.X, 'toarray') else adata.raw.X,
index=adata.obs_names, columns=adata.raw.var_names
)
# TF list from cisTarget resources# Human: https://resources.aertslab.org/cistarget/tf_lists/
tf_names = pd.read_csv('allTFs_hg38.txt', header=None)[0].tolist()
tf_names = [tf for tf in tf_names if tf in expr_matrix.columns]
adjacencies = grnboost2(expr_matrix, tf_names=tf_names, seed=42, verbose=True)
adjacencies.to_csv('adjacencies.tsv', sep='\t', index=False, header=False)
Step 2: Regulon Pruning with RcisTarget
from pyscenic.prune import prune2df, df2regulons
from pyscenic.utils import modules_from_adjacencies
from ctxcore.rnkdb import FeatherRankingDatabase
# cisTarget databases (~10 GB each, download once)# Human: hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather# Mouse: mm10_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather# FeatherRankingDatabase(fname, name) -- `name` is REQUIRED (no default) in ctxcore
dbs = [FeatherRankingDatabase(db, name=os.path.splitext(os.path.basename(db))[0]) for db in [
'hg38_500bp_up_100bp_down.genes_vs_motifs.rankings.feather',
'hg38_10kbp_up_10kbp_down.genes_vs_motifs.rankings.feather'
]]
motif_annotations = 'motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl'# Build co-expression modules as Regulon objects. prune2df reads module.transcription_factor,# which a bare GeneSignature lacks (AttributeError); modules_from_adjacencies applies pySCENIC's# standard top-target/importance thresholds and returns the Regulon objects prune2df expects.
modules = list(modules_from_adjacencies(adjacencies, expr_matrix))
# Prune modules using motif enrichment# NES threshold 3.0 (default); rank_threshold=5000 matches the CLI (prune2df default is 1500).
df_motifs = prune2df(dbs, modules, motif_annotations, rank_threshold=5000, num_workers=8)
regulons = df2regulons(df_motifs)
print(f'Discovered {len(regulons)} regulons')
Step 3: AUCell Activity Scoring
from pyscenic.aucell import aucell
auc_matrix = aucell(expr_matrix, regulons, num_workers=8)
adata.obsm['X_aucell'] = auc_matrix.loc[adata.obs_names].values
adata.uns['regulon_names'] = [r.name for r in regulons]
QC Checkpoint: GRN Inference
defvalidate_grn(regulons, auc_matrix, adata, cell_type_key='cell_type'):
'''
QC gates after GRN inference.
- 50-500 regulons is typical range
- Known lineage TFs should appear (e.g., PAX6 in neurons, GATA1 in erythroid)
- AUCell scores should separate known cell types
'''
n_regulons = len(regulons)
regulon_names = [r.name for r in regulons]
# Gate 1: Regulon countif n_regulons < 50:
print(f'WARNING: Only {n_regulons} regulons. Check TF list or lower NES threshold.')
elif n_regulons > 500:
print(f'WARNING: {n_regulons} regulons found. Consider stricter pruning.')
else:
print(f'OK: {n_regulons} regulons in expected range (50-500)')
# Gate 2: Known TFs present
known_tfs = ['PAX6', 'SOX2', 'GATA1', 'SPI1', 'FOXP3', 'TBX21', 'EBF1']
# df2regulons names regulons 'PAX6(+)' / 'PAX6(-)'; strip the suffix or this gate never fires
regulon_bases = {name.split('(')[0] for name in regulon_names}
found = [tf for tf in known_tfs if tf in regulon_bases]
print(f'Known lineage TFs found: {found}')
# Gate 3: AUCell separates cell typesimport scipy.stats as stats
cell_types = adata.obs[cell_type_key].unique()
iflen(cell_types) >= 2:
ct1_idx = adata.obs[cell_type_key] == cell_types[0]
ct2_idx = adata.obs[cell_type_key] == cell_types[1]
n_differential = 0for i, rname inenumerate(regulon_names[:min(50, len(regulon_names))]):
stat, pval = stats.mannwhitneyu(
auc_matrix.values[ct1_idx, i], auc_matrix.values[ct2_idx, i]
)
if pval < 0.01:
n_differential += 1print(f'Differentially active regulons between top 2 types: {n_differential}/50')
return n_regulons
Path B: SCENIC+ (Multiome)
Step 1: ATAC Topic Modeling with cisTopic
import pycisTopic
from pycisTopic.cistopic_class import create_cistopic_object
from pycisTopic.lda_models import run_cgs_models
# Create cisTopic object from fragments
cistopic_obj = create_cistopic_object(
fragment_matrix=adata_atac.X.T, # cisTopic wants regions x cells; AnnData .X is cells x regions
cell_names=adata_atac.obs_names.tolist(),
region_names=adata_atac.var_names.tolist()
)
# Run LDA topic modeling# n_topics: test range around expected cell types (e.g., 2x number of clusters)
models = run_cgs_models(
cistopic_obj,
n_topics=[10, 20, 30, 40, 50],
n_cpu=8, n_iter=300, random_state=42
)
# evaluate_models plots the model-selection metrics; read the elbow and pass that topic COUNT# as select_model (an int, not True -- True==1 would select a non-existent 1-topic model).from pycisTopic.lda_models import evaluate_models
model = evaluate_models(models, select_model=40, return_model=True)
cistopic_obj.add_LDA_model(model)
Step 2: Enhancer-TF Mapping
import pyranges as pr
from pycistarget.utils import region_names_to_coordinates
from pycistarget.motif_enrichment_cistarget import run_cistarget
from pycisTopic.topic_binarization import binarize_topics
region_bin = binarize_topics(cistopic_obj, method='otsu') # dict of DataFrames keyed by topic (region names in the index)# run_cistarget needs a dict of pyranges.PyRanges, not the raw binarized DataFrames.
region_sets = {topic: pr.PyRanges(region_names_to_coordinates(region_bin[topic].index.tolist()))
for topic in region_bin}
# Run motif enrichment on accessible regions. The first arg is the cisTarget ranking DB: pass the# feather path and run_cistarget instantiates cisTargetDatabase itself. Prebuilt DBs at# https://resources.aertslab.org/cistarget/ . The parameter is spelled `specie`, not `species`.
CTX_DB = '/path/to/hg38_screen_v10_clust.regions_vs_motifs.rankings.feather'
cistarget_results = run_cistarget(
CTX_DB,
region_sets=region_sets,
specie='homo_sapiens',
auc_threshold=0.005,
nes_threshold=3.0,
rank_threshold=0.05,
n_cpu=8
)
Step 3: eGRN Construction
Goal: Assemble eRegulons (TF -> enhancer -> gene triplets) from the multiome data.
Approach: Current SCENIC+ runs topic modeling, motif enrichment, and eGRN construction through one Snakemake pipeline; the deprecated manual create_SCENICPLUS_object/build_grn API (and pre-2024 tutorials) should not be used. See gene-regulatory-networks/multiomics-grn for the full pipeline and the peak-to-gene caveats.
# Scaffold, edit the config (point at fragments, scRNA AnnData, cell-type labels, databases),# then run from inside the Snakemake directory.
scenicplus init_snakemake --out_dir scenicplus_run
# edit scenicplus_run/Snakemake/config/config.yamlcd scenicplus_run/Snakemake && snakemake --cores 16
# Read the resulting direct (high-confidence) eRegulon table (filename is config-/version-# dependent, so resolve it by glob).import glob, pandas as pd
eregulons = pd.read_csv(glob.glob('scenicplus_run/**/eRegulon*direct*.tsv', recursive=True)[0], sep='\t')
print(f'eRegulons: {eregulons["TF"].nunique()} enhancer-driven regulators')
CellOracle Perturbation Simulation
Goal: Predict the direction cells move under a TF knockout, as a hypothesis (direction, not calibrated magnitude).
Approach: CellOracle needs a base GRN (a TF-target scaffold from motif scanning of accessible regions, not the pySCENIC adjacencies), then learns per-cluster weights, propagates a forced expression shift, and projects it onto the cell-state graph. See gene-regulatory-networks/perturbation-simulation for the base-GRN construction and the local-linear / direction-only caveats.
import celloracle as co
import numpy as np
oracle = co.Oracle()
oracle.import_anndata_as_raw_count(adata=adata, cluster_column_name='cell_type',
embedding_name='X_umap')
# Base GRN = motif-scanned accessible regions (preferred) or a prebuilt CellOracle base GRN;# this is NOT the pySCENIC adjacencies. See multiomics-grn / perturbation-simulation.
base_grn = co.data.load_human_promoter_base_GRN() # `version` must match the genome build
oracle.import_TF_data(TF_info_matrix=base_grn)
oracle.perform_PCA()
k = int(0.025 * oracle.adata.n_obs)
oracle.knn_imputation(n_pca_dims=50, k=k, balanced=True, b_sight=k * 8, b_maxl=k * 4)
# Learn context-specific weights, then fit the simulation GRN.
links = oracle.get_links(cluster_name_for_GRN_unit='cell_type', alpha=10)
links.filter_links(p=0.001, weight='coef_abs', threshold_number=2000)
oracle.get_cluster_specific_TFdict_from_Links(links_object=links)
oracle.fit_GRN_for_simulation(alpha=10, use_cluster_specific_TFdict=True)
# Simulate TF knockout (0.0) and project the shift onto the embedding.
oracle.simulate_shift(perturb_condition={'MYC': 0.0}, n_propagation=3)
oracle.estimate_transition_prob(n_neighbors=200, knn_random=True, sampled_fraction=1)
oracle.calculate_embedding_shift(sigma_corr=0.05)
shift = np.sqrt((oracle.delta_embedding ** 2).sum(axis=1))
QC Checkpoint: Perturbation
defvalidate_perturbation(oracle, perturbed_tf, expected_affected_cluster=None):
'''
QC gate: perturbation shifts should match known biology.
- Transition probabilities should show directional shift
- If expected_affected_cluster known, check it shows largest change
'''import numpy as np, pandas as pd
# Shift magnitude per cell from the simulated embedding shift (delta_embedding).
shift = np.sqrt((oracle.delta_embedding ** 2).sum(axis=1))
# observed=True: cell_type is categorical, and the default retains filtered-out categories as NaN rows.# Sort here, not at print time: the gate below reads index[:3], which is category order until sorted.
mean_shift = pd.Series(shift, index=oracle.adata.obs_names).groupby(
oracle.adata.obs['cell_type'].values, observed=True).mean().sort_values(ascending=False)
print(f'Mean shift magnitude by cell type after {perturbed_tf} KO:')
print(mean_shift)
if expected_affected_cluster:
if expected_affected_cluster in mean_shift.index[:3]:
print(f'OK: {expected_affected_cluster} among top affected clusters')
else:
print(f'WARNING: {expected_affected_cluster} not among top affected')
return mean_shift
Complete Pipeline Script
import scanpy as sc
import pandas as pd
from arboreto.algo import grnboost2
from pyscenic.prune import prune2df, df2regulons
from pyscenic.aucell import aucell
from pyscenic.utils import modules_from_adjacencies
from ctxcore.rnkdb import FeatherRankingDatabase
defrun_scenic_pipeline(adata_path, tf_list_path, db_paths, motif_annotations_path, output_prefix):
'''Run complete pySCENIC pipeline.'''
adata = sc.read_h5ad(adata_path)
expr_matrix = pd.DataFrame(
adata.raw.X.toarray() ifhasattr(adata.raw.X, 'toarray') else adata.raw.X,
index=adata.obs_names, columns=adata.raw.var_names
)
tf_names = pd.read_csv(tf_list_path, header=None)[0].tolist()
tf_names = [tf for tf in tf_names if tf in expr_matrix.columns]
print(f'Step 1: GRN inference with {len(tf_names)} TFs')
adjacencies = grnboost2(expr_matrix, tf_names=tf_names, seed=42, verbose=True)
print('Step 2: Regulon pruning')
dbs = [FeatherRankingDatabase(db, name=os.path.splitext(os.path.basename(db))[0]) for db in db_paths]
modules = list(modules_from_adjacencies(adjacencies, expr_matrix))
df_motifs = prune2df(dbs, modules, motif_annotations_path, rank_threshold=5000, num_workers=8)
regulons = df2regulons(df_motifs)
print(f'Discovered {len(regulons)} regulons')
print('Step 3: AUCell scoring')
auc_matrix = aucell(expr_matrix, regulons, num_workers=8)
adata.obsm['X_aucell'] = auc_matrix.loc[adata.obs_names].values
adata.uns['regulon_names'] = [r.name for r in regulons]
adata.write(f'{output_prefix}_scenic.h5ad')
auc_matrix.to_csv(f'{output_prefix}_aucell.csv')
print(f'Pipeline complete: {len(regulons)} regulons, AUCell matrix saved')
return adata, regulons, auc_matrix
Parameter Recommendations
Step
Parameter
Recommendation
GRNBoost2
min_targets
10 (minimum targets per TF module)
RcisTarget
NES threshold
3.0 (standard), 2.5 (permissive)
RcisTarget
databases
Use both 500bp and 10kbp upstream databases
AUCell
auc_threshold
0.05 (fraction of ranked genes)
cisTopic
n_topics
Test 2x expected cell types
CellOracle
n_propagation
3 (default signal propagation steps)
CellOracle
k (imputation)
int(0.025 * n_cells) (CellOracle tutorial rule; ~1250 at 50k cells)
Common Errors
Symptom
Cause
Fix
Near-empty regulons
Species/namespace/DB-vintage mismatch across TF-list, ranking DB, motif2TF
Pin all three to the SAME species + assembly + collection vintage
"Hybrid-state regulator" artifact
Ran GRN on a doublet-contaminated or un-integrated object
Infer on cleaned, doublet-free, batch-controlled cells; run SCENIC once on the integrated object
Inflated adjacencies / everything correlates
Inferred on imputed/smoothed counts
Use RAW counts; imputation only inside CellOracle's simulation scope
Regulon "validated" by TF-expression correlation
AUCell activity <-> TF mRNA circularity
Validate with an orthogonal modality (perturbation/ChIP), not the TF's own mRNA
Label cells first; pycisTopic calls per-celltype pseudobulk peaks
Perturbation magnitudes reported as quantitative
Over-read the direction-only local model
Report direction + a baseline; never a quantitative KO magnitude
Modules called "regulons" without directionality
Skipped the cisTarget ctx pruning step
Run ctx; co-expression modules become regulons only after motif pruning
< 50 regulons / > 500 regulons
Strict pruning-wrong TF list / permissive thresholds
Lower NES to 2.5 (verify species) / raise NES to 3.5
GRNBoost2 hangs or memory error
arboreto dask backend / large dataset
Use bundled multiprocessing; subsample to ~50k cells for GRNBoost2
References
Van de Sande B, Flerin C, Davie K, et al (2020) A scalable SCENIC workflow for single-cell gene regulatory network analysis. Nature Protocols 15:2247-2276. DOI 10.1038/s41596-020-0336-2. (pySCENIC 3-step; multi-run stability.)
Bravo González-Blas C, De Winter S, Hulselmans G, et al (2023) SCENIC+: single-cell multiomic inference of enhancers and gene regulatory networks. Nature Methods 20:1355-1367. DOI 10.1038/s41592-023-01938-4. (eRegulons; needs paired multiome + cell-type labels before peak calling.)
Kamimoto K, Stringa B, Hoffmann CM, et al (2023) Dissecting cell identity via network inference and in silico gene perturbation. Nature 614:742-751. DOI 10.1038/s41586-022-05688-9. (CellOracle; direction-only in-silico perturbation.)
