Gene Regulatory Network Inference (pySCENIC)

Gene Regulatory Network Inference (pySCENIC)

Infer gene regulatory networks (GRNs) de novo from single-cell RNA-seq data using pySCENIC. This workflow discovers transcription factor (TF) regulons directly from expression patterns and calculates cell-level TF activity scores.

When to Use This Skill

Use pySCENIC GRN inference when you need to:

• ✅ Discover TF-target relationships de novo from single-cell RNA-seq data
• ✅ Calculate cell-level TF activity scores for regulatory programs
• ✅ Identify cell-type-specific transcriptional programs
• ✅ Find key regulators driving cell state transitions or differentiation
• ✅ Compare regulatory programs across conditions, tissues, or species

Don't use this skill for:

• ❌ Bulk RNA-seq with few samples → Use functional-enrichment-from-degs or tf-activity workflows
• ❌ Quick TF activity from DE results → Use curated network approaches (faster, less computational)
• ❌ <500 cells → Insufficient for robust GRN inference
• ❌ Limited computational resources → Requires 16GB+ RAM, several hours runtime

Key Concept: Unlike curated network approaches, pySCENIC infers TF-target relationships directly from your data using co-expression analysis, then validates them using motif enrichment in cis-regulatory regions.

The SCENIC Pipeline:

1. GRN Inference (GRNBoost2): Identify co-expression modules linking TFs to potential target genes
2. Regulon Prediction (cisTarget): Prune targets to those with TF binding motifs in promoters
3. Cell Activity Scoring (AUCell): Score each cell for regulon activity

Quick Start

Fastest way to test the workflow (~10-15 minutes):

# Step 1: Load example PBMC data (500 cells)
from load_example_data import load_pbmc3k_example
adata, ex_matrix = load_pbmc3k_example(preprocess=True, subsample=500)

# Step 2: Run complete GRN workflow
from run_grn_workflow import run_complete_grn_workflow
results = run_complete_grn_workflow(
    ex_matrix=ex_matrix,
    tf_list_file="pyscenic_databases/allTFs_hg38.txt",
    database_glob="pyscenic_databases/*.feather",
    motif_annotations_file="pyscenic_databases/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl",
    output_dir="scenic_test_results",
    n_workers=4
)

# Step 3: Visualize and export
from integrate_with_adata import integrate_with_adata
from plot_regulon_visualizations import generate_all_visualizations
from export_all import export_all

adata = integrate_with_adata(adata, results['auc_matrix'], results['regulons'])
generate_all_visualizations(results['auc_matrix'], results['regulons'], adata,
                            output_dir="scenic_test_results/plots")
export_all(results['regulons'], results['auc_matrix'], results['auc_summary'],
           results['adjacencies'], output_dir="scenic_test_results")

Expected output: 20-40 regulons, 13 output files (CSVs, plots, integrated H5AD)

Note: Requires databases to be downloaded first (see Installation → Reference Databases section).

Installation

Required Software

Minimum Python version: Python ≥3.8

Quick install:

pip install pyscenic scanpy pandas numpy loompy networkx seaborn matplotlib reportlab

Note: pySCENIC automatically installs arboreto (for GRNBoost2) and ctxcore (for cisTarget) as dependencies.

Reference Databases (Required)

Download species-specific databases from SCENIC resources:

For human (hg38):

wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl
wget https://resources.aertslab.org/cistarget/tf_lists/allTFs_hg38.txt

For other species: See references/database_downloads.md

Database requirements:

• Motif ranking database (.feather format, ~1-2GB)
• Motif annotation file (.tbl format)
• TF list for your species

Inputs

Required Input

1. Single-cell expression matrix (one of):

   • AnnData object (.h5ad) with raw or normalized counts
   • Loom file (.loom)
   • CSV/TSV matrix (genes × cells)

2. Reference databases (downloaded from SCENIC resources):

   • Motif ranking database (e.g., hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather)
   • Motif annotation file (e.g., motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl)
   • TF list (e.g., allTFs_hg38.txt)

Data Requirements

• Minimum cells: 500 (1,000+ recommended for robust GRN inference)
• Minimum genes: 2,000+ expressed genes
• Memory: 16GB+ RAM (more for large datasets)
• Runtime: 1-4 hours depending on dataset size
• QC: Apply basic quality filtering before SCENIC (remove low-quality cells/genes)

Outputs

Files Generated

Analysis objects (Pickle):

• regulons.pkl - Regulon objects for downstream analysis
  • Load with: regulons = pickle.load(open('regulons.pkl', 'rb'))
  • Required for: Downstream TF activity analysis, network visualization
• auc_matrix.pkl - AUCell activity matrix (cells × regulons)
  • Load with: auc_matrix = pickle.load(open('auc_matrix.pkl', 'rb'))

Results (CSV):

• adjacencies.csv - Raw TF-target co-expression (GRNBoost2 output)
• regulons.csv - Final TF-target relationships after motif pruning
• aucell_matrix.csv - Cell × Regulon activity scores (values 0-1)
• aucell_summary.csv - Per-regulon statistics
• scenic_regulon_summary.csv - Comprehensive regulon summary

Integrated data:

• adata_with_scenic.h5ad - AnnData with integrated regulon activities

Visualizations (PNG + SVG):

• regulon_heatmap.png/.svg - Top regulons by variance
• regulon_network.png/.svg - TF-target network visualization

Reports:

• scenic_report.md - Analysis summary with top regulons
• scenic_analysis_report.pdf - Publication-quality PDF with Introduction, Methods, Results (embedded figures), Conclusions
  • Requires: pip install reportlab (optional — markdown report generated regardless)

Clarification Questions

Before running, confirm:

1. Input Files (ASK THIS FIRST):

   • Do you have single-cell RNA-seq data to analyze?
   • If uploaded: Is this the .h5ad/.loom/matrix file you'd like to use?
   • Expected formats: AnnData (.h5ad), Loom (.loom), CSV/TSV matrix
   • Or use example/demo data? Use load_example_data.py for PBMC 3k (~2,700 cells, 30-45 min test)

2. Species?

   • Human (hg38) - most common, databases readily available
   • Mouse (mm10) - databases available
   • Other species - check database availability first

3. Dataset size?

   • 500-2,000 cells → Minimum viable, ~1-2 hours
   • 2,000-10,000 cells → Good, ~2-3 hours
   • 10,000+ cells → Excellent, may need subsampling, 3-4 hours

4. Do you have the cisTarget databases downloaded?

   • Yes → Provide paths to .feather, .tbl, and TF list files
   • No → Will download during setup (~2-3GB, 10-15 min)

5. What outputs do you need?

   • Regulons only (TF-target relationships)
   • AUCell scores (cell-level TF activity)
   • Both with visualizations (recommended)

Standard Workflow

Note: Run from the OmicsClaw root directory and add the workflow scripts to sys.path:

import sys; import os; sys.path.insert(0, os.path.abspath('knowledge_base/scripts/grn-pyscenic'))

🚨 MANDATORY: USE SCRIPTS EXACTLY AS SHOWN - DO NOT WRITE INLINE CODE 🚨

Step 1 - Load data and run GRN inference:

# Load expression data
# Option A: Load your own data
from load_expression_data import load_expression_data
adata, ex_matrix = load_expression_data("scrnaseq_data.h5ad")

# Option B: Load example PBMC 3k data for testing
# from load_example_data import load_pbmc3k_example
# adata, ex_matrix = load_pbmc3k_example()

# Run GRN inference with GRNBoost2
from run_grn_workflow import run_complete_grn_workflow
results = run_complete_grn_workflow(
    ex_matrix=ex_matrix,
    tf_list_file="allTFs_hg38.txt",
    database_glob="pyscenic_databases/*.feather",
    motif_annotations_file="pyscenic_databases/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl",
    output_dir="scenic_results"
)

DO NOT write inline GRNBoost2, cisTarget, or AUCell code. Use the script.

✅ VERIFICATION: You should see:

• "✓ Data loaded successfully: X cells, Y genes"
• "✓ GRN inference completed: X TF-target pairs"
• "✓ cisTarget pruning completed: X regulons"
• "✓ AUCell scoring completed"

Step 2 - Integrate with AnnData:

from integrate_with_adata import integrate_with_adata
adata = integrate_with_adata(
    adata,
    results['auc_matrix'],
    results['regulons'],
    output_file="scenic_results/adata_with_scenic.h5ad"
)

✅ VERIFICATION: "✓ Integration completed: regulon activities added to adata.obsm['X_aucell']"

Step 3 - Generate visualizations:

from plot_regulon_visualizations import generate_all_visualizations
generate_all_visualizations(
    results['auc_matrix'],
    results['regulons'],
    adata,
    top_n=20,
    output_dir="scenic_results/plots"
)

🚨 DO NOT write inline plotting code (matplotlib, seaborn, etc.). Just use the script. 🚨

The script handles PNG + SVG export with graceful fallback.

✅ VERIFICATION: "✓ All visualizations generated successfully!"

Step 4 - Export results:

from export_all import export_all
export_all(
    regulons=results['regulons'],
    auc_matrix=results['auc_matrix'],
    auc_summary=results['auc_summary'],
    adjacencies=results['adjacencies'],
    output_dir="scenic_results"
)

DO NOT write custom export code. Use export_all().

✅ VERIFICATION: "=== Export Complete ===" Note: If reportlab is installed, a PDF report (scenic_analysis_report.pdf) is also generated.

⚠️ CRITICAL - DO NOT:

• ❌ Write inline GRNBoost2/cisTarget code → STOP: Use run_complete_grn_workflow()
• ❌ Write inline plotting code → STOP: Use generate_all_visualizations()
• ❌ Write custom export code → STOP: Use export_all()
• ❌ Try to install system dependencies → Script checks availability

⚠️ IF SCRIPTS FAIL - Script Failure Hierarchy:

1. Fix and Retry (90%) - Install missing package, re-run script
2. Modify Script (5%) - Edit the script file itself, document changes
3. Use as Reference (4%) - Read script, adapt approach, cite source
4. Write from Scratch (1%) - Only if genuinely impossible, explain why

NEVER skip directly to writing inline code without trying the script first.

Common Issues

Suggested Next Steps

After completing pySCENIC analysis:

1. Identify key regulators: Focus on high-variance regulons with cell-type-specific activity
2. Validate regulons: Compare with literature, ChIP-seq data, or perturbation experiments
3. Downstream analysis:
   • Differential regulon activity between conditions (Scanpy/Seurat)
   • Trajectory analysis with regulon dynamics (RNA velocity + SCENIC)
   • Integration with other modalities (ATAC-seq, ChIP-seq)
4. Functional enrichment: Analyze target genes of top regulons
5. Network analysis: Identify TF-TF interactions and regulatory hierarchies

Related Skills

• scrnaseq-scanpy-core-analysis - Upstream: Single-cell preprocessing and clustering
• scrnaseq-seurat-core-analysis - Upstream: Alternative single-cell preprocessing (R)
• functional-enrichment-from-degs - Related: Pathway analysis of regulon targets
• de-results-to-plots - Related: Visualizing differential activity results

References

• Aibar et al. (2017). SCENIC: single-cell regulatory network inference and clustering. Nature Methods. doi:10.1038/nmeth.4463
• Van de Sande et al. (2020). A scalable SCENIC workflow for single-cell gene regulatory network analysis. Nature Protocols. doi:10.1038/s41596-020-0336-2
• Huynh-Thu et al. (2010). Inferring regulatory networks from expression data using tree-based methods. PLoS ONE. doi:10.1371/journal.pone.0012776
• pySCENIC Documentation: https://pyscenic.readthedocs.io/
• SCENIC Resources: https://resources.aertslab.org/cistarget/