Tutorial

Identifying Gene Regulatory Networks from scRNA-seq using CellOracle

Difficulty: Advanced Time: 25 min read

Introduction

While single-cell RNA sequencing (scRNA-seq) provides a snapshot of gene expression, understanding the underlying regulatory mechanics—how transcription factors (TFs) control target genes—requires advanced network modeling. CellOracle is a machine learning tool that infers Gene Regulatory Networks (GRNs) and simulates how cells behave when specific genes are knocked out or overexpressed.

Architecture Overview


graph LR
    scRNA(scRNA-seq Data) --> Net(CellOracle Net Object)
    ATAC(scATAC-seq / Motif Data) --> BaseGRN(Base GRN)
    BaseGRN --> Net
    Net --> Fit(Ridge Regression Fit)
    Fit --> InferredGRN(Inferred Active GRN)
    InferredGRN --> Oracle(CellOracle Simulation Object)
    scRNA --> Oracle
    Oracle --> Shift(Simulate KO / Overexpression)
    Shift --> UMAP(Project Shift on UMAP)

Prerequisites

  • Python 3.9+
  • CellOracle, Scanpy
  • A preprocessed scRNA-seq dataset (e.g., in AnnData format)
  • Base GRN data (e.g., TF binding motif data from Cicero/ATAC-seq)

Step 1: Install Dependencies

pip install celloracle scanpy

Step 2: Load and Prepare Data

CellOracle requires both a scRNA-seq gene expression matrix and a Base GRN (which dictates which TFs can potentially bind to which target gene promoters).

import celloracle as co
import scanpy as sc

# Load preprocessed scRNA-seq data
adata = sc.read_h5ad("processed_pbmc.h5ad")

# Instantiate a Net object
net = co.Net(
    tensor_dictionary=None, # Used for scATAC-seq integration
    adata=adata,
    cluster_column_name="cluster_labels" # Column defining cell types
)

Step 3: Inferring the Gene Regulatory Network

We load a base GRN and then fit a predictive model (regularized Bayesian Ridge Regression) for each target gene. The model infers the active regulatory connections in each specific cell cluster.

# Load a Base GRN (e.g., from mouse or human promoter motif scans)
base_grn = co.data.load_human_promoter_base_GRN()

# Add the base GRN to the Net object
net.add_tfinfo_dictionary(base_grn)

# Fit the GRN models for all clusters
net.fit_All_genes(
    p=0.001, # p-value threshold for connections
    alpha=10, # Regularization strength
    verbose=True
)

# Extract the inferred GRN for a specific cluster (e.g., 'B cells')
b_cell_grn = net.cluster_specific_GRNs['B cells']
print(b_cell_grn.head())

Step 4: Simulating Transcription Factor Knockout & Overexpression

Once the GRN is established, CellOracle can simulate the transcriptomic shift resulting from altering a specific transcription factor. We project this shift back onto our UMAP embedding to visualize the simulated cell state trajectory.

# Instantiate an Oracle object for simulation
oracle = co.Oracle()
oracle.import_anndata_as_raw_count(adata, cluster_column_name="cluster_labels")
oracle.import_TF_data(TF_info_matrix=b_cell_grn)

# Perform PCA and KNN
oracle.perform_PCA()
oracle.knn_imputation(n_pca_dims=50, k=50)

# 1. Simulate Knockout (KO) of the 'PAX5' transcription factor (a known B cell regulator)
oracle.simulate_shift(TF_to_simulate="PAX5", val=0.0)

# Project the simulated shift onto the UMAP
oracle.estimate_transition_prob(n_neighbors=200, knn_random=True, sampled_fraction=1)
oracle.calculate_embedding_shift(sigma_corr=0.05)

# Plot the KO vector field
co.plot.plot_quiver(oracle, scale=5, title="PAX5 Knockout Simulation")

# 2. Simulate Overexpression of PAX5
# We can set the expression value artificially high (e.g., 5.0)
oracle.simulate_shift(TF_to_simulate="PAX5", val=5.0)
oracle.estimate_transition_prob(n_neighbors=200, knn_random=True, sampled_fraction=1)
oracle.calculate_embedding_shift(sigma_corr=0.05)
co.plot.plot_quiver(oracle, scale=5, title="PAX5 Overexpression Simulation")

Conclusion

CellOracle transforms descriptive single-cell data into a predictive in silico laboratory. By identifying key regulatory nodes and simulating their disruption or over-activation, researchers can uncover novel therapeutic targets for complex diseases directly from computational models before committing to expensive wet-lab experiments.