Bio Gene Regulatory Networks Perturbation Simulation

Version Compatibility

Reference examples tested with: anndata 0.10+, matplotlib 3.8+, numpy 1.26+, pandas 2.2+, scanpy 1.10+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• CLI: <tool> --version then <tool> --help to confirm flags

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Perturbation Simulation

"Predict what happens if I knock out this transcription factor" → Simulate TF perturbation effects on cell identity by combining a base GRN from accessible chromatin with learned regulatory weights from scRNA-seq, then propagating the perturbation signal to predict cell state shifts.

• Python: celloracle.Oracle() for GRN construction and perturbation simulation

Simulate transcription factor perturbation effects on cell state using CellOracle. Integrates GRN inference from scRNA-seq with base GRN from chromatin accessibility to predict cell identity shifts from TF knockouts or overexpression.

CellOracle Overview

CellOracle constructs a GRN by combining:

1. A base GRN from accessible chromatin regions + motif scanning (defines possible TF-target links)
2. scRNA-seq expression data (learns active regulatory weights)

The base GRN can come from scATAC-seq, bulk ATAC-seq, or published chromatin data. CellOracle does NOT require paired multiome data -- any source of accessible regions works.

Installation

pip install celloracle

Step 1: Base GRN from Accessible Regions

From scATAC-seq Peaks

import celloracle as co
import pandas as pd
import numpy as np

# Load peak data (BED format: chr, start, end)
peaks = pd.read_csv('atac_peaks.bed', sep='\t', header=None, names=['chr', 'start', 'end'])

# Scan peaks for TF binding motifs using CellOracle's built-in scanner
# Uses gimmemotifs internally
tfi = co.motif_analysis.TFinfo(peak_data_frame=peaks, ref_genome='hg38')

# Scan for motifs
tfi.scan(fpr=0.02)  # false positive rate for motif matching

# Filter and format as base GRN
tfi.filter_motifs_by_score(threshold=10)
tfi.make_TFinfo_dataframe_and_target_gene_dataframe()

base_grn = tfi.to_dataframe()
base_grn.to_parquet('base_grn.parquet')
print(f'Base GRN: {len(base_grn)} TF-target links')

From Published Chromatin Data

# CellOracle provides pre-built base GRNs for common cell types
# Download from: https://github.com/morris-lab/CellOracle/wiki/
base_grn = co.data.load_mouse_scATAC_atlas_base_GRN(
    organism='Mouse',
    tissue='whole_brain'
)

Step 2: GRN Construction from scRNA-seq

import scanpy as sc
import celloracle as co

adata = sc.read_h5ad('clustered.h5ad')

# Ensure data is preprocessed: normalized, log-transformed, with PCA and clustering
oracle = co.Oracle()
oracle.import_anndata_as_raw_count(
    adata=adata,
    cluster_column_name='cell_type',
    embedding_name='X_umap'
)

# Load base GRN
base_grn = pd.read_parquet('base_grn.parquet')
oracle.import_TF_data(TF_info_matrix=base_grn)

# Fit GRN models per cluster
# Uses regularized linear regression (Bayesian Ridge) to learn TF-target weights
oracle.perform_PCA()
oracle.knn_imputation(n_pnn=30, balanced=True, b_sight=3000, b_maxl=1500)

# Fit GRN for all clusters
links = oracle.get_links(cluster_name_for_GRN_unit='cell_type', alpha=10, verbose_level=0)

# Filter links by statistical significance
# p-value threshold for keeping TF-target connections
links.filter_links(p=0.001, weight='coef_abs', threshold_number=2000)

# Inspect top regulatory connections
links.links_dict['T_cell'].sort_values('coef_abs', ascending=False).head(20)

Step 3: Perturbation Simulation

Knockout Simulation

Goal: Predict how cells change state when a transcription factor is knocked out by simulating the perturbation through the learned GRN.

Approach: Set the target TF expression to zero, propagate the effect through the regulatory network for n steps, estimate transition probabilities to neighboring cell states, and compute embedding shifts that quantify predicted cell identity changes.

# Simulate TF knockout (set expression to 0)
oracle.simulate_shift(perturb_condition={'GATA1': 0.0}, n_propagation=3)

# Get perturbation scores
oracle.estimate_transition_prob(n_neighbors=200, knn_random=True, sampled_fraction=1)
oracle.calculate_embedding_shift(sigma_corr=0.05)

# Perturbation score: magnitude of predicted cell state shift
# Higher score = more affected by the perturbation
perturbation_scores = oracle.adata.obsm['delta_embedding']
shift_magnitude = np.sqrt((perturbation_scores ** 2).sum(axis=1))
oracle.adata.obs['GATA1_KO_shift'] = shift_magnitude

Overexpression Simulation

# Simulate TF overexpression (set to high value)
# Value is relative to max observed expression
oracle.simulate_shift(perturb_condition={'PAX5': 3.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)

Multi-TF Perturbation

# Simulate multiple TF perturbations simultaneously
oracle.simulate_shift(
    perturb_condition={'GATA1': 0.0, 'SPI1': 0.0},  # double knockout
    n_propagation=3
)
oracle.estimate_transition_prob(n_neighbors=200, knn_random=True, sampled_fraction=1)
oracle.calculate_embedding_shift(sigma_corr=0.05)

Visualization

Quiver Plot (Vector Field)

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8, 8))

# Quiver plot shows predicted direction of cell state change
oracle.plot_quiver(
    ax=ax, scale=30,
    color=oracle.adata.obs['cell_type'],
    plot_whole_cells=True
)
ax.set_title('GATA1 KO - predicted cell state shifts')
plt.savefig('gata1_ko_quiver.pdf', bbox_inches='tight')

Gradient Plot

fig, ax = plt.subplots(1, 2, figsize=(16, 8))

# Perturbation score on embedding
sc.pl.embedding(oracle.adata, basis='umap', color='GATA1_KO_shift',
                cmap='Reds', ax=ax[0], show=False, title='Shift magnitude')

# Cell type reference
sc.pl.embedding(oracle.adata, basis='umap', color='cell_type',
                ax=ax[1], show=False, title='Cell types')

plt.savefig('gata1_ko_gradient.pdf', bbox_inches='tight')

Systematic TF Screen

Goal: Rank candidate transcription factors by their predicted impact on cell fate to prioritize perturbation experiments.

Approach: Loop knockout simulations over a list of TFs, compute the mean and max embedding shift magnitude for each, and rank by overall cell state disruption.

# Screen multiple TFs to find drivers of cell fate
tfs_to_screen = ['GATA1', 'SPI1', 'CEBPA', 'PAX5', 'TCF7', 'RUNX1']

results = {}
for tf in tfs_to_screen:
    oracle.simulate_shift(perturb_condition={tf: 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.adata.obsm['delta_embedding'] ** 2).sum(axis=1))
    results[tf] = {
        'mean_shift': shift.mean(),
        'max_shift': shift.max(),
        'affected_cells': (shift > shift.quantile(0.9)).sum()
    }

screen_df = pd.DataFrame(results).T.sort_values('mean_shift', ascending=False)
print(screen_df)

Parameter Reference

Related Skills

• scenic-regulons - TF regulon inference from scRNA-seq with pySCENIC
• multiomics-grn - Enhancer-driven GRNs with SCENIC+
• single-cell/trajectory-inference - Trajectory analysis for cell fate context
• single-cell/perturb-seq - Experimental perturbation data analysis