bio-gene-regulatory-networks-perturbation-simulation

What this skill does

## Version Compatibility

Reference examples tested with: CellOracle 0.18+, scanpy 1.10+, anndata 0.10+; Dynamo (dynamo-release 1.4+), GEARS, CPA (scvi-tools ecosystem) where used.

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"** -> Propagate a forced expression change through a learned GRN (or a learned velocity vector field), then project the shift onto the cell-state manifold to predict the DIRECTION cells move.
- Python: `celloracle.Oracle` for GRN-based in silico KO/overexpression
- Python: `dynamo` for vector-field (Jacobian) perturbation; `GEARS`/`CPA` for response prediction

## The Single Most Important Modern Insight -- These Methods Predict a Direction, Locally, and Inherit Every Error of the Underlying Model

In silico perturbation outputs a **shift vector whose orientation is the deliverable**; its absolute magnitude is uncalibrated. CellOracle propagates an expression *shift* (not absolute counts) and is explicit that it is for **hypothesis generation**, predicting direction not magnitude. The machinery is a **local linear approximation**: CellOracle's per-cluster GRN is a regularized linear model iterated a small fixed number of steps (`n_propagation` default 3), and Dynamo's perturbation is a first-order Taylor expansion of the learned vector field via its Jacobian. So predictions are valid only for **small perturbations near the observed manifold, for one or a few steps** -- they cannot capture long-range feedback, trans effects, or transitions to attractors not already in the data (CellOracle moves probability mass among observed states on the KNN graph; it cannot invent a new cell type). And every prediction **inherits the errors of its substrate**: CellOracle is only as good as its GRN, Dynamo only as good as its RNA-velocity estimate (which has documented reliability problems, Bergen 2021; Gorin 2022).

The 2025 baseline discipline is non-negotiable: for perturbation *response prediction*, deep and foundation models do **not** consistently beat trivial baselines -- Ahlmann-Eltze, Huber & Anders 2025 (*Nat Methods* 22:1657) show that for unseen single perturbations no DL model beats predicting the training **mean** response, and for combinations none beats the **additive** baseline (summing single-gene effects). Always report the mean baseline (unseen singles) and additive baseline (combinations); a method that does not beat them adds no value. And ask the validation question of any result: how many predictions were tested against a real knockout, and were the failures reported?

## Method Taxonomy

| Method | Citation | Mechanism | Predicts | Inherits errors of |
|--------|----------|-----------|----------|--------------------|
| CellOracle | Kamimoto 2023 *Nature* | shift propagated through per-cluster linear GRN, projected on KNN graph | direction of cell-state shift | the GRN (base + regression) |
| Dynamo | Qiu 2022 *Cell* | Jacobian of a learned analytical vector field (Df = J*Dx) | redirected fate / least-action path | RNA-velocity estimate |
| GEARS | Roohani 2024 *Nat Biotechnol* | GNN over gene + GO graphs | unseen single & combinatorial responses, GI type | training Perturb-seq + graphs |
| scGen / CPA | Lotfollahi 2019/2023 | latent-space arithmetic / disentangled composable embeddings | response to known/composed perturbations | training distribution |
| Boolean / ODE (BoolNet, BoolODE) | Müssel 2010 | mechanistic logic/kinetic simulation; attractors | extrapolative fate, multistability | the hand-curated wiring |

Data-driven methods (CellOracle/Dynamo/DL) are genome-scale, local, linear, direction-only; mechanistic Boolean/ODE models are hand-curated and small but nonlinear and extrapolative. State which regime a result belongs to.

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Predict TF KO/OE effect on cell fate from scRNA + accessibility prior | CellOracle | GRN-based shift onto the developmental flow; direction + perturbation score |
| Have high-quality (labeling-based) velocity, want fate redirection | Dynamo | learned vector field + Jacobian; no GRN prior needed |
| Predict response to unseen perturbations / combinations | GEARS or CPA | generalize via gene/GO graphs or composable embeddings -- but check baselines |
| Need mechanistic attractor/multistability reasoning | Boolean/ODE (BoolNet) | only family that extrapolates beyond observed states |
| Build the base GRN that CellOracle perturbs | -> multiomics-grn | base GRN construction lives there |
| Validate predictions experimentally | -> single-cell/perturb-seq | Perturb-seq is the interventional ground truth |

## CellOracle: Build the GRN

**Goal:** Learn context-specific (per-cluster) regulatory weights on top of a base-GRN prior.

**Approach:** Import scRNA-seq and the base GRN, impute, then fit a regularized linear model per cluster and filter links by significance.

```python
import scanpy as sc
import celloracle as co
import pandas as pd

adata = sc.read_h5ad('clustered.h5ad')          # normalized, log, PCA, clustered
oracle = co.Oracle()
oracle.import_anndata_as_raw_count(adata=adata, cluster_column_name='cell_type',
                                   embedding_name='X_umap')
oracle.import_TF_data(TF_info_matrix=pd.read_parquet('base_grn.parquet'))  # prior; see multiomics-grn
oracle.perform_PCA()
oracle.knn_imputation(n_pca_dims=50, k=None, balanced=True, b_sight=3000, b_maxl=1500)

links = oracle.get_links(cluster_name_for_GRN_unit='cell_type', alpha=10)
links.filter_links(p=0.001, weight='coef_abs', threshold_number=2000)
oracle.get_cluster_specific_TFdict_from_Links(links_object=links)
oracle.fit_GRN_for_simulation(alpha=10, use_cluster_specific_TFdict=True)
```

## CellOracle: Simulate and Score

**Goal:** Predict the direction cells move under a TF knockout or overexpression.

**Approach:** Set the gene of interest to 0 (KO) or above-max (OE), propagate the shift a few steps, estimate KNN transition probabilities, compute the per-cell embedding shift, and score it against the developmental vector field by inner product.

```python
import numpy as np

# perturb_condition sets the clamped value: 0.0 = knockout; above observed max = overexpression.
# n_propagation is small BY DESIGN (3); increasing it amplifies linear-approximation error.
oracle.simulate_shift(perturb_condition={'GATA1': 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)

# Perturbation score = inner product of the simulated shift with the developmental flow.
# It is meaningless without a defined development vector field (Gradient_calculator below).
from celloracle.applications import Gradient_calculator
grad = Gradient_calculator(oracle_object=oracle, pseudotime_key='pseudotime')
grad.calculate_p_mass(smooth=0.8, n_grid=40, n_neighbors=200)
grad.calculate_mass_filter(min_mass=0.01, plot=False)   # required before calculate_gradient
grad.calculate_gradient()
shift = np.sqrt((oracle.adata.obsm['delta_embedding'] ** 2).sum(axis=1))   # magnitude is uncalibrated
```

## Dynamo: Vector-Field Perturbation

**Goal:** Predict fate redirection from a learned dynamical system rather than a GRN prior.

**Approach:** Reconstruct the analytical vector field, compute its Jacobian, then apply a genetic perturbation as a local linear response.

```python
import dynamo as dyn

dyn.vf.VectorField(adata, basis='umap')          # learn the analytical vector field first
dyn.vf.jacobian(adata, regulators=['GATA1'