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-2026, including recent critical evaluations of foundation models.
Perturbation modeling is a rally between model and lab — played across four very different courts. The model serves a prediction, the wet lab returns the measured response, and the gap between them is the loss, the benchmark, and the next drill. 四大满贯，四种场地，同一个问题：your prediction only counts if it holds on every surface.
「读懂发球，才接得住回球。」
Every perturbation experiment is a rally: the model serves a predicted response, the cell returns the measured one. And the tour never stays on one surface — hard court, clay, grass, and night sessions each change how the ball behaves, the way new cell types, platforms, and conditions change how a perturbation reads out. Six moments from the practice courts, re-read as the concepts this module's methods and benchmarks formalize.
On the practice court you pick the serve — flat, slice, topspin — and the drill is only as informative as the serves you bother to practice.
Perturbation design. The intervention is the serve: gene knockout, drug, dose, combination. CRISPR screens are the cataloged serve drills — thousands of serves, one court. See the task ↓
A rally player calls the return before the ball crosses the net — from the serve's spin, speed, and the opponent's stance, not from watching it land.
Response prediction. Given control state + serve, call the perturbed state: control + Δ = perturbed. The whole task is calling the Δ before the experiment returns it. See the formulation ↓
On clay at Roland Garros you meet the grinder: serves that look wild but all land in the same corner. Stand in that corner and you look like a genius predictor — until the tour moves to a different surface.
Systematic variation & the mean baseline. Any perturbation shifts the whole transcriptome in one direction, so predicting the perturbed-group mean scores like SOTA without learning anything per-gene. Always check what the corner-stander scores. See the critique ↓
Drilling a hundred known serves on hard court is not the Championships. The first morning on grass, the ball skids through low and late — a serve you have never seen, and your read has to travel.
Unseen-perturbation generalization. The hard split is not random cells but new genes, new combinations, new doses. Most methods' scores drop sharply here; this is the split that counts. See benchmarks ↓
A highlight-reel champion with no footwork loses to the club player who simply returns every ball deep. Reputation is not a return game.
Foundation models vs simple baselines. On standard perturbation benchmarks, large pretrained models often fail to beat train-mean or linear baselines. Scale is not mechanism; check the fundamentals first. See the evidence ↓
Good players debrief every rally — which reads were wrong, by how much, on which serve — and choose tomorrow's drill from the errors, not from the highlight reel.
Honest metrics + active learning. Debrief with centroid-referenced metrics (top-DEG Δ, E-distance), then spend the next experiment where the error is largest: optimal intervention design closes the loop. See the loop ↓
每一杆都是实验，每个回合都在校准 — every serve is a prediction, every return a measurement; the rally is the loop.
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.
Tabular-attention foundation model trained on 149M cells that performs true in-context learning at the cell-set level, reframing perturbation prediction as prompting rather than fine-tuning: give the model a labelled prompt cell-set at inference time and it predicts the response to a new, unseen perturbation with no weight updates.
Both State and STACK are 2026-era transformer approaches, but they represent opposite adaptation paradigms: State fine-tunes a set-based transformer on the target task, while STACK uses in-context learning with no weight updates. For the pretraining objectives, tokenisation choices, and scaling laws behind these and other foundation models, see the Pretraining & Foundation Models module.
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.
CellOracle and CIPHER predict perturbation responses through a fixed or inferred regulatory graph. For methods that instead discover the underlying causal graph from perturbation data (DCDI, ENCO, and related structure-learning approaches), see the Causal Inference module — the two problems share data (Perturb-seq screens) but answer different questions: this module asks "what will happen," causal inference asks "what regulates what."
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 + 3 datasets from other perturbation/readout modalities, 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).
Systematic benchmark pitting five foundation models plus two additional deep-learning methods against deliberately simple linear/additive baselines for single- and double-gene-perturbation transcriptome prediction. None of the deep models beat the baselines.
Arc Institute's recurring benchmark competition for perturbation prediction, aiming to establish rigorous, CASP-style standards for the field. The inaugural 2025 challenge (H1 hESC context generalization) closed with results in early 2026; the 2026 edition raises the bar to zero-shot CRISPRi-response prediction across six held-out cell lines with no released training set.
Standardized benchmarking framework from Altos Labs — unified data loaders, biologically meaningful train/test splits, a modular model zoo, and rank-based metrics designed to complement (and expose the limits of) standard RMSE/R² fit measures.
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) and VAEs (SAMS-VAE, CPA) provide alternative frameworks. Each method claims superiority on different metrics and datasets.
GEARS (dual knowledge-graph GNN) reaches Nature Biotechnology. scFoundation and scGPT are adapted for perturbation tasks.
100M+ cell datasets emerging — Tahoe-100M, X-Atlas/Orion, and Mixscale all released. Physics-inspired methods (CIPHER) and biological priors (scLAMBDA with LLM embeddings, code) show up. Field recognizes need for proper baselines and evaluation.
The Virtual Cell Challenge's inaugural 2025 competition wraps up — the first CASP-style blind evaluation for perturbation prediction — and reveals that even the winning approach needed hand-engineered statistical features layered on deep learning; the harder, zero-shot 2026 edition launches immediately after. STACK (bioRxiv) reframes perturbation prediction as in-context learning: no fine-tuning, just a labelled prompt cell-set at inference time. Geneformer-scaling (Nat Comput Sci) brings 4-bit QLoRA quantization to genome-scale virtual knockout screens, cutting a 316M-parameter screen from $25k/33 days to <$5k/6 days. The field converges on rank-based metrics and mandatory linear baselines following the Ahlmann-Eltze benchmark.
| 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 | 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) |
| State | 2025 | Transformer (Set-based) | 100M+ cells pretrain | State Embedding + State Transition, bidirectional attention over cell populations |
| STACK | 2026 | Transformer (ICL) | 149M cells pretrain | True in-context learning — no fine-tuning, no weight updates |
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:
For Zero-Shot / Few-Shot Transfer (no fine-tuning):
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 — pin the exact version so it matches the PyG wheel below
pip install torch==2.1.0 --index-url https://download.pytorch.org/whl/cu121
# 2. PyTorch Geometric + companion wheels (GEARS' graph backend).
# The wheel URL is locked to torch 2.1.0+cu121 — if you install a
# different torch/CUDA combo above, update this URL to match
# (see https://data.pyg.org/whl/ for the available combinations).
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.