PaperlayerView exact revision on arXiv

Source-linked AI summary

Geometric coherence of single-cell CRISPR perturbations reveals regulatory architecture and predicts cellular stress

Prashant C. Raju

arXiv:2604.16642v3q-bio.QMq-bio.CBq-bio.GNstat.AP

TL;DR

Single-cell CRISPR screens measure perturbation magnitude but not whether cells move coherently through state space. This paper introduces Shesha perturbation stability (S_p), a geometric coherence metric, and shows that stability is robustly related to effect magnitude while revealing regulatory architecture and predicting stress-associated outcomes.

  • Problem

    Existing genome-editing evaluations assess sequence fidelity but provide limited evidence about whether edited cells occupy stable, coherent cellular states.

  • Method

    Shesha perturbation stability (S_p) measures mean cosine similarity between individual cell shift vectors and the mean perturbation direction.

  • Results

    Across five CRISPR datasets spanning over 2,200 perturbations, stability correlated strongly with effect magnitude, with Spearman ρ=0.75–0.97.

  • Takeaways & Limitations

    Geometric stability provides a complementary axis for evaluating perturbation screens, assessing hit reproducibility, and characterizing regulatory architecture shaping cellular responses.

  • Takeaways & Limitations

    The association between geometric instability and UPR activation remains correlative, and the study does not demonstrate that incoherence causes UPR activation.

Abstract

from arXiv · show

Genome engineering has achieved sequence-level precision, yet predicting the transcriptomic state a cell will occupy after perturbation remains open. Single-cell CRISPR screens measure how far cells move, but effect magnitude ignores whether the cells move together. We introduce Shesha perturbation stability ($S_p$), which quantifies directional coherence as the mean cosine similarity between individual cell shift vectors and the mean perturbation direction. Across five CRISPR datasets (2,200+ perturbations), stability correlates with magnitude (Spearman $ρ= 0.75$--$0.97$), but discordant cases expose regulatory architecture: pleiotropic regulators such as CEBPA pay a ``geometric tax,'' producing large but incoherent shifts, while lineage-specific factors such as KLF1 produce coordinated responses. $S_p$ and Song et al.'s perturbation-response score (PS) share partial overlap ($ρ_{\text{partial}} = +0.51$ after controlling for magnitude), but $S_p$ provides significant incremental prediction of UPR pathway activation beyond both PS and magnitude ($p < 10^{-18}$). In a split-half reproducibility assay, $S_p$ predicts directional reproducibility beyond magnitude ($ρ_{\text{partial}} = +0.384$) while PS does not ($ρ_{\text{partial}} = -0.193$), with the advantage consistent across all magnitude strata and both datasets. Geometric instability is independently associated with UPR activation across four datasets. $S_p$ is implemented in the open-source shesha-geometry Python package.

Introduction

Single-cell CRISPR editing achieves sequence-level precision without reliably predicting the resulting cellular state. The Shesha perturbation stability metric addresses this gap by measuring whether perturbed cells move coherently, revealing regulatory architecture beyond effect magnitude.

  • Sequence-level editing precision has outpaced the ability to predict cellular consequences and unintended cell fates.
  • Existing evaluation frameworks measure edit syntax, including indel rates, off-target cleavage, and sequence fidelity, but not whether the resulting state is stable.
  • Single-cell CRISPR screens typically summarize responses by effect magnitude, overlooking how cells move relative to one another.
  • Shesha perturbation stability (S_p) measures the mean cosine similarity between individual cell shift vectors and the mean perturbation direction.
  • Across five single-cell CRISPR datasets comprising over 2,200 perturbations, stability correlated with effect magnitude at Spearman ρ=0.75–0.97.Pleiotropic regulators such as CEBPA produced large but incoherent shifts, whereas lineage-specific factors such as KLF1 produced tightly coordinated responses.

Materials and Methods … Quantifying geometric stability of perturbations

The study analyzed five single-cell CRISPR perturbation datasets using independent preprocessing and multi-stage control matching. It quantified perturbation stability as the directional coherence of cell-level PCA shifts relative to controls, alongside effect magnitude.

  • Datasets: Five datasets spanning CRISPRa, CRISPRi, pooled screening, and genome-scale CRISPRi were analyzed across K562, BMDC, and THP-1 cells.The datasets included 2,200+ perturbations overall, with Replogle et al. contributing 1,832 after filtering.
  • Preprocessing: Each dataset was preprocessed independently with quality filtering, library-size normalization, log transformation, 2,000 highly variable genes, and 50-component PCA.Cells with fewer than 100 detected genes were removed; normalization targets differed for Replogle and Adamson.
  • Control group identification: Controls were identified through case-insensitive term matching, delimiter-aware regular expressions for short tokens, and substring matching for embedded keywords.Dataset-specific rules assigned Replogle labels containing “non-targeting” or beginning with “chr” to control.
  • Quantifying geometric stability of perturbations: For each perturbation, each cell’s shift vector was computed as its PCA coordinates minus the control centroid.The mean perturbation direction was defined as the average of the individual cell shift vectors.
  • Quantifying geometric stability of perturbations: Shesha perturbation stability (S_p) is the mean cosine similarity between individual cell shift vectors and the mean perturbation direction.This definition measures directional coherence among perturbed cells.
  • Quantifying geometric stability of perturbations: S_p near 1 indicates coordinated movement, whereas S_p near 0 indicates scattering across expression space.Effect magnitude (M_p) is the Euclidean norm of the mean shift vector, M_p = ∥d̄∥.
  • Quantifying geometric stability of perturbations: Perturbations with fewer than 50 cells, or fewer than 10 cells in Dixit, were excluded, while all 8 Adamson perturbations were retained.The metric was robust to Euclidean, Mahalanobis, and k-NN distance choices, with wide bootstrap confidence intervals for Adamson.

Discordance

Discordance measures whether a perturbation’s stability is lower or higher than expected from its magnitude, using a nonlinear LOESS fit. The resulting biological conclusions remain robust across residual methods, although high-magnitude gene rankings can shift.

  • Discordance: Discordance is defined as the residual from the magnitude–stability relationship.It identifies perturbations where magnitude and stability diverge.
  • Discordance: LOESS with bandwidth fraction = 0.3 models the nonlinear relationship at low magnitudes.This avoids relying on ordinary least-squares regression where signal-to-noise constrains coherence.
  • Discordance: Positive discordance indicates lower-than-predicted stability, whereas negative discordance indicates higher-than-predicted stability.Positive values fall below the magnitude–stability curve; negative values lie above it.
  • Discordance: Biological conclusions are robust to residual-method choice, although specific gene rankings shift for high-magnitude perturbations.Linear and rank-based residual methods are compared in the SI Appendix.

Functional diversity of differentially expressed genes … Pathway-level stress scoring

The analyses assess functional diversity, perturbation-response scoring, pseudobulk heterogeneity, directional reproducibility, and pathway-level cellular stress using complementary transcriptomic and geometric measures. Stress analyses additionally control for effect magnitude and examine canonical markers and within-stratum differences.

  • Functional diversity of differentially expressed genes: Top-k differentially expressed genes (k ∈ {25, 50, 100}) were ranked by absolute log-fold change and tested for GO Biological Process enrichment with g:Profiler.The analysis used the Norman 2019 dataset and assessed sensitivity across three gene-list sizes.
  • Perturbation-response score comparison: Song et al.’s perturbation-response score was computed using Euclidean distance, Mahalanobis distance, and a Python port of scMAGeCK constrained optimization.The three approaches increased in fidelity from PCA-space distance to covariance-adjusted distance and constrained optimization.
  • Comparison with pseudobulk heterogeneity metrics: Pseudobulk heterogeneity was represented by η2, the proportion of PCA variance explained by perturbation identity across the first 50 principal components.Spearman correlations between η2 and LOESS-residual discordance were computed for each dataset.
  • Split-half reproducibility assay: A split-half assay tested directional reproducibility by correlating mean shift vectors from two randomly partitioned cell halves across 50 independent splits.The assay covered Replogle (n = 1,832) and Norman (n = 236) perturbations.
  • Pathway-level stress scoring: Cellular stress was scored using four MSigDB Hallmark pathways: Unfolded Protein Response, p53 Pathway, Apoptosis, and Reactive Oxygen Species.Scores were computed with scanpy.tl.score_genes using curated gene lists.
  • Pathway-level stress scoring: Pathway scores used 72–78 UPR, 32–51 p53, 53–73 Apoptosis, and 45–55 ROS genes per dataset, with correlations controlling for effect magnitude and Benjamini-Hochberg correction.Both raw Spearman and partial correlations were reported across 16 pathway-dataset combinations.
  • Pathway-level stress scoring: Stress analyses also assessed DDIT3, ATF4, XBP1, and HSPA5, tested quadrant depletion, stratified Replogle perturbations by GO molecular function, and compared HSPA5 within magnitude quartiles.The HSPA5 comparison used a one-sided Mann-Whitney U test between high- and low-discordance groups.

scGPT validation · Statistical analysis

The magnitude–stability relationship was tested beyond linear PCA using scGPT Whole Human embeddings across three datasets, with statistical uncertainty and cross-dataset effects assessed using prespecified resampling and mixed-effects procedures.

  • scGPT validation: scGPT validation recomputed stability and magnitude in Whole Human pretrained embeddings for Norman, Dixit, and Replogle datasets.The analysis tested whether the relationship depended on linear PCA embeddings.
  • scGPT validation: The scGPT validation used raw counts rather than log-normalized counts as embedding input.Embedding used embed_data() with deterministic settings.
  • scGPT validation: Stability and magnitude were computed identically in the scGPT and PCA pipelines.This preserved the metric definitions while changing the embedding representation.
  • Statistical analysis: Confidence intervals used bootstrap resampling with 10,000 iterations, seed 320, and the percentile method.The same procedure defined uncertainty intervals for the reported analyses.
  • Statistical analysis: Cross-dataset generalization was assessed using a linear mixed-effects model with dataset as a random effect.Fixed effects were magnitude, spread, and sample size.
  • Statistical analysis: Partial correlations between stability and stress markers controlled for effect magnitude, with p-values two-sided unless otherwise noted.The full mixed-effects specification was provided in the SI Appendix.

Results

Across five CRISPR datasets, perturbation stability strongly tracked effect magnitude while discordant perturbations revealed distinct regulatory architectures. Stability also predicted directional reproducibility and cellular stress beyond magnitude and perturbation-response score, with geometric instability associated with UPR activation across datasets.

  • Magnitude and stability: ρ=0.968 across 2,254 perturbations after dataset-wise z-scoring, confirming a robust cross-dataset magnitude–stability relationship.Within individual datasets, Spearman correlations ranged from ρ=0.746 to ρ=0.985 across CRISPRa and CRISPRi screens.
  • Regulatory architecture: CEBPA perturbations showed large but discordant shifts, whereas KLF1 and its combinations consistently produced highly concordant responses.In Norman, CEBPA-family perturbations occupied the discordant extreme at magnitudes Mp=8–10, while KLF1 was the most concordant single-gene perturbation under both residual methods.
  • Perturbation-response score: ρ=−0.203 after controlling for magnitude and PS, showing that Sp significantly predicted UPR pathway activation beyond both predictors.Sp and real PS were moderately correlated after magnitude control (partial ρ=+0.507), but each retained unique information.
  • Cellular stress: partial ρ=−0.214 in Replogle and partial ρ=−0.231 in Dixit, while stability showed a sign-consistent negative association with UPR score across all four testable datasets.The Papalexi association was directionally consistent but underpowered, and the Norman estimate was not significant.

Discussion

Across datasets and embedding methods, geometric stability complements perturbation magnitude by revealing regulatory architecture, reproducibility, and stress-related biological costs. The framework supports hit prioritization while requiring caution about correlation, metric scope, embedding choices, and residual-method sensitivity.

  • Regulatory architecture: Pleiotropic regulators produce large but incoherent responses, whereas lineage-specific factors produce coherent responses aligned with developmental trajectories.This magnitude–stability decoupling motivates the term “geometric tax.”
  • Regulatory architecture: CEBP-family perturbations engage more than twice as many GO Biological Process categories as KLF1 combinations (p=0.013), but functional diversity does not explain discordance globally (ρ=−0.119, p=0.11).The result supports broad transcriptional scope as one mechanism without making it a universal explanation.
  • Population-level signal: ρ=−0.276 (p<10^-33): Sp-discordance and eta2 are negatively correlated in Replogle, consistent with geometric incoherence diluting aggregate transcriptomic signal.eta2 is the proportion of PCA variance explained by perturbation identity.
  • Biological costs: Geometric instability is associated with UPR activation across four datasets, but the relationship remains correlational, with effects explaining only 4–5% of residual variance.The study does not demonstrate that geometric incoherence causes UPR activation.
  • Applications: At equivalent magnitude, higher Sp predicts directionally reproducible phenotypes across all magnitude strata, providing a secondary criterion for prioritizing CRISPR screen hits.The analysis frames geometric stability as a practical ranking criterion when candidate effects have similar sizes.
  • Limitations: The magnitude–stability relationship persists in scGPT embeddings, but PCA remains primary; nonlinear manifold methods may reveal additional structure, and discordance rankings depend on residual method.Top discordant genes differ substantially across residual methods for high-magnitude perturbations.

Extended Methods

The analysis used independently preprocessed single-cell CRISPR datasets, standardized control assignment, and 50-dimensional PCA embeddings. Shesha perturbation stability was computed from cell-level perturbation shifts relative to dataset-specific control centroids, with magnitude-adjusted discordance defined by regression residuals.

  • Dataset preprocessing: Five datasets were accessed through pertpy 1.0.4 and preprocessed independently to prevent batch effects.The pipeline included quality filtering, library-size normalization, and log transformation.
  • Dataset preprocessing: Downstream stability and magnitude computations were performed on a 50-dimensional PCA embedding.This embedding was used after dataset-specific preprocessing.
  • Control assignment: Controls were identified using case-insensitive exact matches, delimiter-aware regex, substring matching, and dataset-specific label handling.The protocol accommodated heterogeneous labeling conventions across datasets.
  • Stability metric: Shesha perturbation stability S_p measures coherence between individual cell shifts from the control centroid and the mean perturbation direction.For perturbation p, cell shifts were defined in PCA coordinates relative to the same-dataset control centroid.
  • Stability metric: Positive discordance indicates lower stability than predicted by magnitude, whereas negative discordance indicates higher stability than predicted.Discordance is the standardized residual from ordinary least-squares regression of S_p on M_p.

Robustness Analyses

Robustness analyses show that the magnitude–stability relationship persists across distance metrics, PCA dimensionalities, random seeds, and perturbation subsets. Observed stability also far exceeds a hyperspherical null model, supporting a biological coherence signal.

  • Distance metrics: All tested distance metrics produced comparable or stronger magnitude–stability correlations than the standard Euclidean metric.Whitened Mahalanobis-scaled and k-NN matched methods were equal or higher across all four datasets, potentially reducing batch and control heterogeneity noise.
  • PCA dimensionality: Spearman ρ increased monotonically with PCA dimensionality, from 0.949 to 0.969 for Norman and 0.781 to 0.869 for Dixit.The default 50 components yielded ρ=0.959 for Norman and ρ=0.844 for Dixit, within the upper range across tested dimensionalities.
  • Random seeds: 0.99997 and 0.99963 were the mean cross-seed stability-ranking correlations for Norman and Dixit across 15 random seeds.Magnitude–stability Spearman ρ varied by less than 0.001 across seeds: Norman 0.959–0.960 and Dixit 0.844–0.844.
  • Perturbation influence: 0.0019 and 0.0106 were the maximum absolute changes in ρ after removing any single perturbation for Norman and Dixit, respectively.The leave-one-out ranges were [0.959, 0.961] for Norman and [0.840, 0.854] for Dixit, indicating that no individual perturbation drove the overall correlation.
  • Null model: 0.05–0.85 were the observed stability values, far exceeding the S_p≈0 expected under uniformly sampled 50-dimensional hyperspherical shift vectors.The null variance is inversely proportional to n_p, supporting a biological rather than statistical coherence signal.

Mixed-Effects Model

A linear mixed-effects model tested whether stability generalizes across datasets after accounting for effect magnitude, within-perturbation variance, and cell count. Magnitude explained 11.2 times more variance than sample size, while near-zero dataset random-effect variance supported cross-dataset generalization.

  • Mixed-Effects Model: The model predicted S_p from effect magnitude, within-perturbation expression variance, and cell count, with a dataset-level random intercept and residual error.All predictors were z-scored within each dataset before fitting.
  • Mixed-Effects Model: 11.2 times more variance in stability was explained by magnitude than by sample size (|β_1/β_3| = 11.2).This indicates that the magnitude–stability relationship was not driven by differential cell sampling.
  • Mixed-Effects Model: Near-zero dataset random-effect variance indicated that the magnitude–stability relationship generalized across datasets rather than being driven by one screen.Full results were reported in Table S2.

Extended Stress Marker Analysis

Extended stress-marker analysis found marker-specific and context-dependent associations with perturbation stability. DDIT3 correlations varied by modality and cell type, whereas HSPA5 showed the most consistent negative partial correlations in CRISPRi datasets.

  • Analysis framework: Stability–stress associations were evaluated using raw and magnitude-controlled Spearman correlations across four markers and three datasets.Median stability and stress expression defined high-stability/high-stress quadrants, whose counts were tested against independence expectations with a one-sided binomial test.
  • DDIT3: DDIT3 raw correlations were −0.365 in Dixit CRISPRi and +0.382 in Replogle CRISPRi, indicating contributions from cell type and experimental design.The sign also differed between CRISPRa and CRISPRi; after magnitude control, all three DDIT3 partial correlations were negative, but only Replogle survived with a small effect.
  • HSPA5: HSPA5 showed negative partial correlations in both CRISPRi datasets, Dixit and Replogle, but a null partial correlation in CRISPRa Norman.This made HSPA5 the most consistent marker across modalities.

Extended Replogle Analysis

The Replogle 2022 genome-scale CRISPRi screen independently validates the observed discordance pattern at unprecedented scale: pleiotropic regulators show lower-than-expected coherence, whereas narrowly acting factors show higher-than-expected coherence.

  • Extended Replogle Analysis: 1,832 perturbations and 310,385 cells provide independent validation of the discordance pattern in the Replogle 2022 genome-scale CRISPRi screen.The analysis filtered perturbations for ≥50 cells.
  • Extended Replogle Analysis: Pleiotropic regulators GATA1, CHMP3, and AQR cluster below the regression line, while narrowly acting factors LSG1, ISG20L2, and KRI1 cluster above.This reproduces the pattern observed in Norman.

scGPT Validation Protocol

The scGPT validation used the official Whole Human checkpoint on raw-count inputs with specified embedding, reproducibility, and geometric-analysis settings. Across datasets, the rank order was preserved, although scGPT correlations were slightly lower, especially for Replogle.

  • Embedding protocol: The official scGPT Whole Human checkpoint generated cell embeddings from raw counts using embed_data() with gene_col="index", batch_size=64, and use_fast_transformer=False.The checkpoint comprised best_model.pt, vocab.json, and args.json; inputs were not log-normalized before embedding.
  • Reproducibility: All scGPT computations used deterministic mode with Python, NumPy, and PyTorch seeds set to 320.The PyTorch seed was applied on both CPU and CUDA.
  • Geometric analysis: Stability used cosine distances and magnitude used Euclidean distances, computed with perturbation_name and control labels; confidence intervals used 10,000 bootstrap resamples with seed 320.The analyses used shesha.bio.compute_stability() and compute_magnitude() with percentile bootstrap intervals.
  • Dataset scope: scGPT analysis covered Norman 2019, Dixit 2016, and Replogle 2022, while Adamson 2016 was omitted from main-text figures and Papalexi 2021 was excluded.Adamson 2016 had n=8 and was retained in the output CSV; Papalexi 2021 was excluded because of pertpy loader incompatibility with current mudata versions.
  • Validation result: The dataset rank order remained Norman > Replogle > Dixit, with the largest scGPT correlation decrease in Replogle (Δρ= −0.119).The passage attributes the lower scGPT correlations to nonlinear embedding resolving manifold structure that PCA collapses.

Combinatorial Analysis

In the Norman 2019 dataset, combinatorial perturbations were more stable than single-gene perturbations even after accounting for magnitude, with lineage-specific combinations outperforming pleiotropic-factor combinations. Nonlinear correction also changed discordance rankings, reducing high-magnitude artifacts and highlighting subunits of essential multiprotein complexes.

  • Combinatorial versus single-gene perturbations: 0.460 versus 0.306: combinatorial perturbations had significantly higher mean stability than single-gene perturbations in Norman 2019.The dataset included 131 combinatorial and 105 single-gene perturbations (Mann-Whitney U=2,637, p=4.1 × 10^-16).
  • Combinatorial versus single-gene perturbations: 0.973 and 0.919: magnitude-stability correlations remained strong within single-gene and combinatorial categories, respectively.Their regression slopes were 0.089 and 0.058, showing that higher combinatorial stability was not simply due to larger effect magnitudes.
  • Regulatory-factor combinations: Lineage-specific combinations such as KLF1+SET and KLF1+TGFBR2 were more stable than pleiotropic combinations such as CEBPA+JUN and CEBPA+CEBPB.This reproduced the discordance pattern observed in the single-gene analysis.
  • Discordance-method robustness: ≤1/10: linear discordance rankings overlapped minimally with rank-based and LOESS methods, whereas rank-based and LOESS rankings agreed substantially better.Overall correlations between rank-based and LOESS methods were ρ=0.940 and 0.925; linear residuals overestimated discordance for high-magnitude perturbations.
  • Discordance-method robustness: 5 of the top 8: CEBPA variants occupied most top discordant positions under linear residuals but largely disappeared after LOESS correction.CEBPA+JUN fell from rank 1 to rank 172, while KLF1 rankings were nearly unchanged at 236 and 235.
  • Biological interpretation of discordance: Five of ten LOESS-discordant genes were spliceosome components, while three were ESCRT or membrane-remodeling factors and one was a proteasome subunit.All encoded subunits of large, essential multiprotein complexes; CHMP3 was the only original linear top-three gene retaining a top-five rank.
Loading 2604.16642v3…