Version Compatibility
Reference examples tested with: pySCENIC 0.12+, arboreto 0.1.6+, ctxcore 0.2+, pycisTopic 2.0+, pycistarget 1.0+, SCENIC+ 1.0a1 (Snakemake CLI), CellOracle 0.18+, anndata 0.10+, pandas 2.2+, scanpy 1.10+, scipy 1.12+
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 |
Pipeline Overview
Processed AnnData (QC'd, normalized, clustered)
|
+----- RNA only? -------> Path A: pySCENIC (3-step)
| |
| v
| [1. GRNBoost2] ----> TF-target adjacencies
| |
| v
| [2. RcisTarget] ---> Regulon pruning (motif enrichment)
| |
| v
| [3. AUCell] -------> Regulon activity scoring
|
+----- Multiome? -------> Path B: SCENIC+
| |
| v
| [1. cisTopic] -----> Topic modeling on ATAC
| |
| v
| [2. pycistarget] --> Enhancer-TF mapping
| |
| v
| [3. SCENIC+] ------> eGRN construction
|
+---> [CellOracle Perturbation Simulation] (either path)
|
v
Perturbation scores + predicted cell state shifts
Path A: pySCENIC (RNA-Only)
Step 1: GRN Inference with GRNBoost2
import scanpy as sc
import pandas as pd
from arboreto.algo import grnboost2
adata = sc.read_h5ad('processed.h5ad')
expr_matrix = pd.DataFrame(
adata.raw.X.toarray() if hasattr(adata.raw.X, 'toarray') else adata.raw.X,
index=adata.obs_names, columns=adata.raw.var_names
)
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
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'
modules = list(modules_from_adjacencies(adjacencies, expr_matrix))
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
def validate_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]
if 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)')
known_tfs = ['PAX6', 'SOX2', 'GATA1', 'SPI1', 'FOXP3', 'TBX21', 'EBF1']
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}')
scipy.stats stats
cell_types = adata.obs[cell_type_key].unique()
(cell_types) >= :
ct1_idx = adata.obs[cell_type_key] == cell_types[]
ct2_idx = adata.obs[cell_type_key] == cell_types[]
n_differential =
i, rname (regulon_names[:(, (regulon_names))]):
stat, pval = stats.mannwhitneyu(
auc_matrix.values[ct1_idx, i], auc_matrix.values[ct2_idx, i]
)
pval < :
n_differential +=
()
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
cistopic_obj = create_cistopic_object(
fragment_matrix=adata_atac.X.T,
cell_names=adata_atac.obs_names.tolist(),
region_names=adata_atac.var_names.tolist()
)
models = run_cgs_models(
cistopic_obj,
n_topics=[10, 20, 30, 40, 50],
n_cpu=8, n_iter=300, random_state=42
)
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')
region_sets = {topic: pr.PyRanges(region_names_to_coordinates(region_bin[topic].index.tolist()))
for topic in region_bin}
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.
scenicplus init_snakemake --out_dir scenicplus_run
cd scenicplus_run/Snakemake && snakemake --cores 16
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 = co.data.load_human_promoter_base_GRN()
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)
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)
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
def validate_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 = np.sqrt((oracle.delta_embedding ** 2).sum(axis=1))
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
def run_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() if hasattr(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=)
regulons = df2regulons(df_motifs)
()
()
auc_matrix = aucell(expr_matrix, regulons, num_workers=)
adata.obsm[] = auc_matrix.loc[adata.obs_names].values
adata.uns[] = [r.name r regulons]
adata.write()
auc_matrix.to_csv()
()
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 |
| ctx step returns empty | Missing/mismatched motif2TF annotation (most common) | Confirm the .tbl matches the DB vintage + species |
| SCENIC+ peaks miss rare types | Called peaks before/without cell-type labels | 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.)
Related Skills
- gene-regulatory-networks/scenic-regulons - pySCENIC implementation details
- gene-regulatory-networks/multiomics-grn - SCENIC+ enhancer-driven GRNs
- gene-regulatory-networks/perturbation-simulation - CellOracle details
- single-cell/clustering - Upstream cell type annotation
- single-cell/preprocessing - QC and normalization before GRN inference
- atac-seq/single-cell-atac - scATAC preprocessing for SCENIC+ Multiome input
- atac-seq/co-accessibility - Cicero / SCENIC+ cis-regulatory connections
- atac-seq/enhancer-gene-linking - ABC / ENCODE-rE2G enhancer-gene mapping
- atac-seq/motif-deviation - chromVAR for TF motif accessibility
- workflows/scrnaseq-pipeline - Upstream: provides the cleaned, annotated RNA object for Path A (pySCENIC)
- workflows/multiome-pipeline - Upstream: provides the paired RNA+ATAC object for Path B (SCENIC+)
