SCENIC tutorial

· 3 min read · Tutorials

SCENIC tutorial

TL;DR: SCENIC infers Transcription Factor (TF) activity from single-cell data. Use pySCENIC for Python/CLI workflows. Key steps: GRN (co-expression), CTX (motif enrichment), and AUCell (activity scoring).

What is SCENIC?

Single-Cell rEgulatory Network Inference and Clustering.

Most scRNA-seq analysis looks at gene expression (mRNA levels). SCENIC looks at regulon activity (is the TF and its targets active?).

Why use it?

  • Identify master regulators driving cell states.
  • Robust to dropouts (aggregates signal across target genes).
  • Works across species (Human, Mouse, Fly).

SCENIC Workflow

Part 1: Installation (The Right Way)

pySCENIC dependency management can be tricky. To avoid common errors (like np.object deprecation), use this exact conda setup:

# Create environment
conda create -y -n pyscenic python=3.10

# Activate
conda activate pyscenic

# Install main packages
pip install pyscenic scanpy loompy matplotlib seaborn

# Critical dependency fixes
pip install "numpy<1.24" "pandas<2.0"

Verify installation:

pyscenic --help

Part 2: Prepare Input Data

1. Expression Matrix (.loom)

Export your Seurat/Scanpy object to .loom.

From Seurat (R):

library(Seurat)
library(SeuratDisk)

# Load object
obj <- readRDS("seurat_obj.rds")

# Convert to loom
as.loom(obj, filename = "input.loom", verbose = FALSE)

From Scanpy (Python):

import scanpy as sc

adata = sc.read("adata.h5ad")
adata.write_loom("input.loom")

2. Download Databases

You need three types of files from cisTarget resources:

A. TF List:

wget https://resources.aertslab.org/cistarget/tf_lists/allTFs_hg38.txt

B. Motif Rankings (Gene-based):
Warning: Large files (1GB+)

# Human 10kb up/down
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

C. Motif Annotation:

wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl

Part 3: Running the Pipeline (CLI)

Run these steps in a screen or tmux session (takes hours).

Step 1: GRN (Gene Regulatory Network)

Infers co-expression modules.

pyscenic grn \
    --num_workers 12 \
    --output adj.tsv \
    --method grnboost2 \
    input.loom \
    allTFs_hg38.txt

Step 2: CTX (Cis-Target enrichment)

Prunes modules for direct targets (motif enrichment).

pyscenic ctx \
    adj.tsv \
    hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
    --annotations_fname motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl \
    --expression_mtx_fname input.loom \
    --output regulons.csv \
    --mask_dropouts \
    --num_workers 12

Step 3: AUCell (Activity Scoring)

Scores regulon activity in each cell.

pyscenic aucell \
    input.loom \
    regulons.csv \
    --output output_scenic.loom \
    --num_workers 12

Part 4: Visualization in Python

Load the results into Scanpy for visualization.

import scanpy as sc
import loompy as lp
import pandas as pd

# Load original data
adata = sc.read_loom("input.loom")

# Load SCENIC results (AUCell matrix)
lf = lp.connect("output_scenic.loom", mode='r', validate=False)
auc_mtx = pd.DataFrame(lf.ca.RegulonsAUC, index=lf.ca.CellID)
lf.close()

# Add to AnnData
adata.obsm['X_aucell'] = auc_mtx

1. UMAP on Regulons

Cluster cells based on TF activity rather than gene expression.

# Compute neighbors on AUCell matrix
sc.pp.neighbors(adata, use_rep='X_aucell')
sc.tl.umap(adata)
sc.pl.umap(adata, color='celltype', title='Regulon-based UMAP')

2. Visualize Specific TFs

# Copy specific regulon scores to .obs for plotting
regulon = 'PAX5(+)'
adata.obs[regulon] = auc_mtx[regulon]

sc.pl.umap(adata, color=regulon, cmap='viridis')

3. Regulon Specificity

# Matrix plot of top regulons per cluster
sc.pl.matrixplot(adata, 
                 var_names=['PAX5(+)', 'EBF1(+)', 'TCF3(+)'], 
                 groupby='celltype', 
                 use_raw=False)

Troubleshooting

"AttributeError: module 'numpy' has no attribute 'object'"

Cause: NumPy 1.24+ deprecated np.object.
Fix: Downgrade NumPy:

pip install "numpy<1.24"

"TypeError: Must supply at least one delayed object"

Cause: dask version incompatibility.
Fix:

pip install dask==2023.5.0 distributed==2023.5.0

"loompy.connect validate=False error"

Cause: Newer loompy versions validate heavily.
Fix: Regenerate the loom file using scanpy instead of SeuratDisk if possible, or use the validate=False flag in python scripts.

Resources

See also:

Originally posted on WeChat Official Account: bioinfo_studyroom

Gene Regulatory Network TF Activity