Context-Dependent Prediction of CRISPR Perturbation Responses in Melanoma

A rigorous benchmark study of perturbation × condition single-cell transcriptional response prediction on the Frangieh et al. (2021) melanoma Perturb-CITE-seq dataset.

Data: Frangieh et al. 2021, Nature Genetics — Perturb-CITE-seq (melanoma), via scPerturb

---

TL;DR — main finding

Evaluation design, and especially the choice of metric, determines the scientific conclusion.

• On absolute DE-gene expression profiles every model scores highly (R² DE 0.91–0.98), but this metric is weakly discriminative: a no-effect control baseline (predict the condition-matched control mean) already reaches R² DE = 0.969, above the linear additive model (0.909). Absolute scores are dominated by baseline expression levels, not by the perturbation effect.

• Evaluating the effect itself (the delta, change from control) flips the picture. On held-out perturbation × condition combinations:

  model	delta Pearson DE	delta R² DE
  control-only	n/a (predicts zero change)	−0.28 ± 0.15
  linear additive	0.05 ± 0.07	−3.39 ± 0.77
  random forest	0.52 ± 0.11	0.05 ± 0.15
  CPA (off-the-shelf, 3 seeds)	0.24 ± 0.21	−0.62 ± 0.72

• A simple random forest is the only model that reliably captures the effect direction, and its advantage is specific to unseen perturbation × condition combinations (it is concentrated in IFN-γ-context hits such as B2M, CD274, AEBP1, ACTA2, CTSD). A specialized deep model (CPA) was not better here: high absolute scores but a low and highly seed-unstable delta (−0.00 to 0.37 across seeds), underlining that specialized models must be assessed across seeds and against strong baselines under effect-level metrics.

The full write-up is in reports/report_final.md.

---

Repository structure

crispr_perturb_project/
├── README.md
├── environment.yml          # local conda env (Python 3.11) for the classical pipeline
├── data/                    # FrangiehIzar2021_RNA.h5ad (not committed; see "Data")
├── src/
│   ├── data_loading.py      # load + QC + tidy perturbation/condition fields
│   ├── splits.py            # the three OOD splits
│   ├── evaluate.py          # unified metrics (absolute + delta Pearson/R² on DE genes)
│   ├── models.py            # control-only, linear additive, random forest
│   ├── run_benchmark.py     # one model × split benchmark table
│   ├── run_multiseed.py     # multi-seed benchmark with mean ± std
│   ├── hyperparam_search.py # random-forest hyperparameter selection on validation pairs
│   ├── plot_results.py      # benchmark figures + B2M case study + absolute-vs-delta scatter
│   ├── plot_umap.py         # UMAP data-exploration figures
│   ├── analyze_improvements.py
│   └── run_cpa.py           # CPA (cpa-tools), same evaluation footing (needs a GPU)
├── notebooks/
│   ├── cpa_colab.ipynb      # CPA on Colab/GPU (self-contained)
│   └── CPA_Colab_README.md
├── scripts/
│   ├── run_cpa.sbatch       # SLURM job for CPA on an HPC cluster
│   └── CLUSTER_SETUP.md
├── reports/                 # the written report
└── results/                 # benchmark tables and figures

---

Setup

conda env create -f environment.yml      # creates the `crispr-pert` env (Python 3.11)
conda activate crispr-pert

The classical pipeline (control / linear / random forest + all evaluation and figures) runs on a laptop. CPA needs a GPU and a separate install — see Running CPA.

Data

The harmonized Frangieh Perturb-CITE-seq RNA matrix (218,331 cells × 23,712 genes) comes from scPerturb on Zenodo:

mkdir -p data
wget -O data/FrangiehIzar2021_RNA.h5ad \
  "https://zenodo.org/record/10044268/files/FrangiehIzar2021_RNA.h5ad?download=1"

Key columns: perturbation (249 target genes incl. control), perturbation_2 (condition: Control / IFNγ / Co-culture). For tractable benchmarking each run uses a reproducible 30,000-cell subsample and the top 1,000 highly variable genes.

---

Reproduce the benchmark

Canonical multi-seed run (3 models incl. the control baseline, cell resampling, the min_eval_cells filter, a fixed pre-registered holdout set, and the delta metrics):

conda run -n crispr-pert python -m src.run_multiseed \
  --seeds 0 1 2 3 4 --vary-sample-seed \
  --path data/FrangiehIzar2021_RNA.h5ad --sample-n 30000 \
  --n-hvg 1000 --rf-trees 50 --rf-max-depth 8 --min-eval-cells 15 \
  --models control linear_additive random_forest \
  --holdout-pairs "B2M@@IFNγ" "CD274@@IFNγ" "AEBP1@@IFNγ" "ACTA2@@IFNγ" "CTSD@@IFNγ" \
  --holdout-perts ACTA2 AEBP1 CD274 B2M APOD \
  --outdir results/multiseed_final

Random-forest hyperparameter selection (on validation pairs, test pairs untouched):

conda run -n crispr-pert python -m src.hyperparam_search \
  --path data/FrangiehIzar2021_RNA.h5ad --sample-n 30000 --sample-seed 0 \
  --n-hvg 1000 --test-holdout-n 5 --val-holdout-n 5 \
  --n-estimators-grid 10 30 50 --max-depth-grid 8 16 None \
  --outdir results/hyperparam

Figures (benchmark bars, delta-Pearson panel, B2M case study, absolute-vs-delta scatter):

conda run -n crispr-pert python -m src.plot_results \
  --summary-path results/multiseed_final/benchmark_summary_mean_std.csv \
  --figdir results/figures_final \
  --path data/FrangiehIzar2021_RNA.h5ad --sample-n 30000 --sample-seed 0 --seed 0 \
  --n-hvg 1000 --rf-trees 50 --rf-max-depth 8 \
  --holdout-pairs "B2M@@IFNγ" "CD274@@IFNγ" "AEBP1@@IFNγ" "ACTA2@@IFNγ" "CTSD@@IFNγ" \
  --case-perturbation B2M --case-condition IFNγ

conda run -n crispr-pert python -m src.plot_umap \
  --path data/FrangiehIzar2021_RNA.h5ad --sample-n 0 --figdir results/figures_final

Per-group improvement analysis:

conda run -n crispr-pert python -m src.analyze_improvements \
  --scores-path results/multiseed_final/benchmark_group_scores_all_seeds.csv \
  --split holdout_pert_condition --metric r2_de --top-n 10 \
  --outdir results/analysis_final

Canonical outputs: results/multiseed_final/, results/figures_final/ (incl. benchmark_delta_pearson_de.* and delta_scatter.*), results/analysis_final/.

---

Running CPA

CPA needs a GPU and a separate environment (it is not in environment.yml). src/run_cpa.py reuses the same data loading / HVG / evaluation, so its row is apples-to-apples with the classical benchmark and automatically reports the delta metrics.

• Cluster (SLURM): see scripts/CLUSTER_SETUP.md and scripts/run_cpa.sbatch.
• Colab / rented GPU: see notebooks/CPA_Colab_README.md and notebooks/cpa_colab.ipynb.

Install gotchas (cpa-tools 0.8.8): use Python 3.10/3.11 (not 3.12), do not pin an old pyyaml, and pin pyarrow==12.0.1 (ray 2.9.3 needs the removed PyExtensionType).

Example run on a GPU box:

python -m src.run_cpa \
  --path data/FrangiehIzar2021_RNA.h5ad --sample-n 30000 --max-epochs 300 --seed 0 \
  --outdir results/cpa_seed0
# default --holdout-pairs is the pre-registered B2M/CD274/AEBP1/ACTA2/CTSD @ IFNγ

---

Methods at a glance

Three OOD splits (src/splits.py): random (leaky reference), holdout_perturbation (unseen perturbation), holdout_pert_condition (unseen perturbation × condition — the on-target test).

Metrics (src/evaluate.py): Pearson and R² on the top differentially expressed genes, reported both as absolute (conventional, but inflated — kept as a diagnostic) and as delta (change from the condition-matched control — the effect itself, used for all conclusions). Cross-split comparability is enforced with --min-eval-cells 15; uncertainty is estimated by resampling the subsample across 5 seeds.

Models (src/models.py): control-only baseline (no-effect floor), linear additive baseline, tuned random forest; plus CPA (src/run_cpa.py) as the deep compositional model.

Compute

Classical models + all evaluation/figures run on a laptop (arm64 macOS tested). CPA was run on a rented GPU (RTX 4090) / can run on an HPC cluster. The full Perturb-CITE-seq matrix is stored CSC; data_loading.load_data converts to CSR (per-cell ops are pathologically slow on CSC).