A comprehensive guide to state-of-the-art approaches combining single-cell transcriptomics, CRISPR screens, and machine learning to predict how cells respond to perturbations. Covering methods from 2016-2025, including recent critical evaluations of foundation models.
Perturbation modeling is like a tennis rally. The model serves a prediction about how a cell will respond to a CRISPR edit or drug. The wet lab returns the ball with the observed cellular response. The difference between the predicted and measured return tells us how the model should improve.
Understanding how cells respond to perturbations (genetic knockouts, drug treatments, etc.) is fundamental to biology and medicine. However:
Computational prediction methods aim to:
Pioneering VAE-based approach using vector arithmetic in latent space for perturbation prediction. Assumes homogeneous responses and additive perturbation effects.
Combines compositionality, disentanglement, and interpretability through sparse global perturbation variables. Each perturbation targets sparse subset of latent dimensions.
Conditional VAE with additive composition of perturbation and covariate embeddings. Predicts counterfactual distributions but requires training data for each perturbation.
Integrates dual knowledge graphs (gene coexpression + Gene Ontology perturbation network) with deep learning.
Uses Input Convex Neural Networks (ICNNs) to learn transport maps between unpaired control and perturbed cell distributions. Predicts single-cell-level responses.
Foundation model pretrained on 33+ million cells. Uses masked gene modeling with binned expression values. Fine-tuned for perturbation prediction with adapter layers.
Large-scale foundation model with read-depth-aware pretraining on 50M cells. Designed for multiple downstream tasks including perturbation prediction.
Set-based transformer trained on 100M+ perturbation cells across 70 cell lines and 167M observational cells. Arc Institute's first virtual cell model using bidirectional attention over cell populations.
Two-stage pipeline: (1) infer base GRN from scATAC-seq, (2) refine with scRNA-seq using Bayesian/Bagging Ridge regression. Linear propagation through GRN.
Physics-inspired approach using linear response theory: Δx = Σu where Σ is covariance matrix from control cell fluctuations. Bayesian inference with Horseshoe priors.
Industrial-scale perturbation atlases provide massive training data for foundation models. These represent significant experimental advances in scale and throughput.
Largest public perturbation atlas: 100M cells, 50 cancer cell lines, 1,100 drugs. Vevo Therapeutics' Mosaic platform with Parse GigaLab sequencing. Open-sourced on Arc Virtual Cell Atlas and HuggingFace.
Harmonized repository of 44 datasets from 25 publications. 32 CRISPR + 9 drug datasets, average 160K cells per dataset. Introduces Energy statistics for perturbation quantification.
Landmark genome-wide study: K562 (>2.5M cells, 9,867 genes targeted with CRISPRi), RPE1 Essential. Key finding: only 41% of perturbations show measurable transcriptome-wide effects.
Fix-Cryopreserve-ScRNAseq platform achieving 8M cells targeting all human protein-coding genes. sgRNA abundance as dose-dependent proxy (R=0.91 with KD efficiency).
Continuous quantification of perturbation strength replacing binary classification. Systematic pathway profiling across 6 cell lines × 5 signaling contexts (2.6M cells).
Arc Institute's recurring benchmark competition for perturbation prediction. Aims to establish rigorous standards like CASP did for protein structure prediction. First challenge focuses on context generalization (H1 hESC).
Population-Level (Fit-Based):
Distribution-Based:
Rank-Based (Critical for Mode Collapse Detection):
Perturb-seq (Dixit 2016, Adamson 2016) establishes experimental paradigm for pooled perturbation screens with single-cell readouts. CROP-seq (Datlinger 2017) provides alternative guide capture method. ~200K cells per study.
scGen (2019) pioneers VAE approach. Norman et al. (2019) provide key combinatorial perturbation dataset. Replogle et al. (2022) scale to genome-wide with >2.5M cells. CellOracle (2023) leverages GRNs.
Neural OT (CellOT), VAEs (SAMS-VAE, CPA), and GNNs (GEARS) provide alternative frameworks. Each method claims superiority on different metrics and datasets.
scFoundation, scGPT adapted for perturbation tasks.
100M+ cell datasets emerging. Physics-inspired methods (CIPHER) and biological priors (scLAMBDA with LLM embeddings) show up. Field recognizes need for proper baselines and evaluation.
| Method | Year | Architecture | Scale | Key Innovation |
|---|---|---|---|---|
| Perturb-Seq | 2016 | Experimental | ~10K-200K cells | CRISPR + scRNA-seq pooled screens |
| scGen | 2019 | VAE | ~50K cells | Latent space vector arithmetic |
| CellOT | 2023 | Neural OT | ~50K cells | ICNN transport maps, single-cell predictions |
| SAMS-VAE | 2023 | VAE | ~120K cells | Sparse additive mechanism shifts |
| CPA | 2023 | VAE | ~100K cells | Compositional perturbation + covariates |
| GEARS | 2023/2024 | GNN | ~600K cells | Dual knowledge graphs (GO + coexpression) |
| CellOracle | 2023 | GRN | ~100K cells | scATAC-seq informed GRN inference |
| scGPT | 2024 | Transformer (100M) | 33M cells pretrain | Generative pretraining (BUT often underperforms) |
| scFoundation | 2024 | Transformer (100M) | 50M cells | Read-depth-aware (BUT mode collapse issues) |
| CIPHER | 2025 | Statistical Physics | ~1.4M cells | Linear response theory (Δx = Σu) |
Note: Publication years reflect official journal publication dates. Some papers (e.g., GEARS) appeared online earlier than print publication.
Step 1: Always implement simple baselines first:
For Drug Response Prediction:
For Genetic Perturbations:
For Interpretability:
The walkthrough below runs a complete, purely computational perturbation-prediction experiment end to end: load a Perturb-seq dataset (Norman 2019 CRISPRa or Adamson 2016 CRISPRi), split it so that entire perturbations — not cells — are held out, train GEARS on the seen perturbations, predict post-perturbation expression for the unseen ones, and score the predictions with MSE on the top differentially expressed genes plus Pearson delta against a control-mean baseline. Everything is retrospective: no new experiments, only re-analysis of published screens.
GEARS pins an older PyTorch Geometric stack, so give it its own environment rather than sharing one with scGPT or scvi-tools. Install PyTorch (matching your CUDA version) first, then the PyG companion wheels, then GEARS. scGen and CPA live in the scvi-tools ecosystem and can share a second environment; pertpy provides the dataset loaders and the differential-expression utilities used for scoring.
# conda / mamba recommended
conda create -n pertpred python=3.10 -y
conda activate pertpred
# 1. PyTorch first — pick the build that matches your driver
pip install torch --index-url https://download.pytorch.org/whl/cu121
# 2. PyTorch Geometric + companion wheels (GEARS' graph backend)
pip install torch_geometric
pip install torch_scatter torch_sparse \
-f https://data.pyg.org/whl/torch-2.1.0+cu121.html
# 3. Core stack
pip install cell-gears scanpy anndata pertpy
# Separate env for the VAE/latent-arithmetic family (scvi-tools pins conflict)
# conda create -n pertvae python=3.10 -y && conda activate pertvae
# pip install scgen cpa-tools scvi-tools scanpyHardware. Norman 2019 (~110k cells, 287 perturbations) trains on a single 16 GB GPU in roughly 1–2 hours for 20 epochs; CPU-only training is possible but an order of magnitude slower. Peak host RAM is dominated by the dense gene × perturbation tensors GEARS builds — budget 32–64 GB. Genome-scale Replogle K562 (>2.5M cells) needs a compute node with 128 GB+ RAM and is usually subsampled before any exploratory run.
AnnData — adata.X holds normalized, log1p-transformed expression (GEARS expects log-space, not raw counts); adata.obs carries the perturbation label; adata.var must contain a gene_name column so perturbation identifiers can be matched to measured genes.adata.obs["condition"] as "GENE+ctrl" for singles, "GENEA+GENEB" for doubles, and the literal string "ctrl" for unperturbed cells. Malformed labels are the single most common load failure.PertData: a gene–gene co-expression graph over the measured transcriptome, and a perturbation–perturbation similarity graph derived from Gene Ontology term overlap. The GO graph is what lets GEARS place a never-seen perturbed gene near functionally related seen genes and thus emit a non-trivial prediction for it.pert_data.set2conditions maps train/val/test to lists of condition strings. This is the object to inspect (and to save alongside results) when you want the split to be reproducible.adata.uns["rank_genes_groups_cov_all"] stores, per condition, the ranked list of differentially expressed genes versus control. All headline metrics are computed on the top 20 of these, not on the full transcriptome.model.predict([["GENE"], ["GENEA", "GENEB"]]) returns a dict keyed by condition, each value a length-n_genes vector of predicted mean post-perturbation expression in log space.Load Norman 2019 through the GEARS dataloader, create a split that holds out entire unseen perturbations, train, predict, and evaluate. The evaluation block is the part worth reading closely: it computes MSE on the top-20 DE genes, Pearson delta, and the control-mean baseline in the same loop so the three numbers are always reported together.
import numpy as np
from scipy.stats import pearsonr
from gears import PertData, GEARS
# 1. Load a published Perturb-seq screen. "norman" = Norman et al. 2019
# (K562, CRISPRa, 100 singles + 131 doubles); "adamson" = Adamson et al.
# 2016 (K562, CRISPRi, single-gene only). PertData downloads the h5ad,
# normalizes/log1p's it, and builds the gene co-expression graph.
pert_data = PertData("./data")
pert_data.load(data_name="norman")
# 2. Hold out entire PERTURBATIONS, not cells. Held-out cells from a
# perturbation the model already trained on measure nothing but denoising.
# split="simulation" mimics the realistic case: some test perturbations
# involve genes never perturbed during training (unseen singles), and some
# doubles combine a seen gene with an unseen one (0/1/2-seen subgroups).
pert_data.prepare_split(split="simulation", seed=1, train_gene_set_size=0.75)
pert_data.get_dataloader(batch_size=32, test_batch_size=128)
print({k: len(v) for k, v in pert_data.set2conditions.items()})
# 3. Train. GEARS embeds each perturbed gene via the GO-derived perturbation
# graph, so an unseen gene inherits signal from its functional neighbours
# instead of falling back to the training mean.
model = GEARS(pert_data, device="cuda")
model.model_initialize(hidden_size=64)
model.train(epochs=20, lr=1e-3)
model.save_model("gears_norman")
# 4. Predict post-perturbation expression for held-out conditions.
# Singles are single-element lists; doubles are two-element lists.
test_conds = pert_data.set2conditions["test"]
preds = model.predict([c.replace("+ctrl", "").split("+")
for c in test_conds if c != "ctrl"])
# 5. Reference quantities from the raw AnnData
adata = pert_data.adata
ctrl_mean = np.asarray(
adata[adata.obs["condition"] == "ctrl"].X.mean(axis=0)
).ravel()
gene_names = adata.var["gene_name"].values
de_map = adata.uns["rank_genes_groups_cov_all"] # condition -> ranked DE genes
# 6. Evaluate. Three numbers, always together:
# (a) MSE on the top-20 DE genes — all genes dilutes the signal to noise
# (b) Pearson delta — cor(pred - ctrl, true - ctrl)
# (c) the control-mean baseline — "predict no change at all"
rows = []
for cond in test_conds:
if cond == "ctrl" or cond not in preds:
continue
true_mean = np.asarray(
adata[adata.obs["condition"] == cond].X.mean(axis=0)
).ravel()
pred_mean = preds[cond]
de20 = de_map[cond][:20]
idx = np.where(np.isin(gene_names, de20))[0]
mse_de = np.mean((pred_mean[idx] - true_mean[idx]) ** 2)
mse_base = np.mean((ctrl_mean[idx] - true_mean[idx]) ** 2)
# Pearson DELTA, not raw Pearson. Raw cor(pred, true) over all genes runs
# ~0.98 for essentially any model, because most genes do not respond to
# the perturbation and both vectors are dominated by the shared baseline
# expression profile. Subtracting the control mean removes that shared
# component and scores only what the model actually had to get right:
# the direction and magnitude of the CHANGE.
r_delta = pearsonr(pred_mean[idx] - ctrl_mean[idx],
true_mean[idx] - ctrl_mean[idx])[0]
r_raw = pearsonr(pred_mean[idx], true_mean[idx])[0]
rows.append((cond, mse_de, mse_base, r_delta, r_raw))
mse_de, mse_base, r_delta, r_raw = (np.nanmean([r[i] for r in rows])
for i in (1, 2, 3, 4))
print(f"MSE (top-20 DE) : {mse_de:.4f}")
print(f"MSE control-mean baseline: {mse_base:.4f} <- must beat this")
print(f"Pearson delta : {r_delta:.3f}")
print(f"Pearson raw (inflated) : {r_raw:.3f}")
# 7. Break the score out by generalization difficulty. A model can look strong
# overall while being useless on the only split that matters.
subgroups = pert_data.subgroup["test_subgroup"] # e.g. unseen_single, combo_seen0/1/2
for grp, conds in subgroups.items():
sub = [r for r in rows if r[0] in set(conds)]
if sub:
print(f"{grp:15s} n={len(sub):3d} "
f"Pearson delta={np.nanmean([r[3] for r in sub]):.3f}")prepare_split behavior can change across GEARS versions; save pert_data.set2conditions to disk with the results so the same held-out perturbations can be re-scored later by a different method.