Bio Gene Regulatory Networks Perturbation Simulation

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

# 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

# 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