Bio Crispr Screens Perturb Seq Analysis

Version Compatibility

Reference examples tested with: Pertpy 0.6+, SCEPTRE 0.10+ (R / katsevich-lab/sceptre), Mixscape via Seurat 4.3+ or Pertpy, scanpy 1.10+, anndata 0.10+, pandas 2.2+, numpy 1.26+, scipy 1.12+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show pertpy scanpy anndata
• R: packageVersion('sceptre'); ?sceptre; ?Seurat::PrepLDA

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Single-Cell Perturb-Seq Analysis

"Analyze a single-cell pooled CRISPR perturbation screen" -> Assign sgRNAs to cells, filter unperturbed escapers, normalize counts, fit per-gene differential expression conditioned on perturbation, and rank perturbations by their molecular effect.

• Python: pertpy unified framework for Mixscape + SCEPTRE-via-R + differential expression
• R: sceptre for low-MOI NB GLM + permutation testing
• Python/R: Seurat::MixscapeLDA and downstream

Experimental Architecture Comparison

Decision rule: For genome-wide CRISPRi screens, Replogle's CRISPRi + 10X 3' direct-capture protocol is the gold standard (>2.5M cells; the genome-scale K562 screen targeted ~9,866 expressed genes in Replogle 2022). For protein readout, Perturb-CITE-seq. For chromatin, Perturb-multiome. For low-throughput pilot, original Dixit Perturb-seq.

MOI and sgRNA Assignment

The central technical challenge: Each cell must receive exactly one sgRNA (otherwise the perturbation is undefined). At MOI 0.3, ~26% of cells get ≥1 sgRNA, but 4% get ≥2; the cells with multiple sgRNAs must be filtered or analyzed as combinatorial perturbations.

Assignment workflow:

1. Detect sgRNA reads per cell: From the sgRNA library prep (direct capture or 3'UTR barcode), count reads per sgRNA per cell.
2. Threshold: Most pipelines use 10+ reads of one sgRNA to assign that perturbation.
3. Multiplets: Cells with 2+ sgRNAs at >10 reads each are either multi-perturbed (analyzable as combinatorial) or doublets.
4. Doublet detection: Use scDblFinder, Scrublet, or AMULET (multiome) to identify doublets independently from sgRNA assignment.

Goal: Assign a single perturbation identity (or 'multiplet'/'none') to every cell from the sgRNA counts matrix.

Approach: Threshold per-cell sgRNA reads at ≥10 (Pertpy convention); cells exceeding the threshold for exactly one sgRNA are assigned that perturbation; cells with multiple sgRNAs above threshold are flagged as multiplets for filtering or combinatorial analysis.

# sgRNA assignment via threshold counting
def assign_sgrna(adata, sgrna_counts_layer='sgrna_counts', threshold=10):
    '''Per-cell sgRNA assignment. Returns single assignment or 'multiplet'/'none'.'''
    import numpy as np
    counts = adata.layers[sgrna_counts_layer]  # cells x sgRNAs
    above_thresh = counts >= threshold
    n_sgrna_per_cell = above_thresh.sum(axis=1)
    assignments = np.where(
        n_sgrna_per_cell == 0, 'none',
        np.where(n_sgrna_per_cell == 1,
                  [adata.var_names[i] for i in counts.argmax(axis=1)],
                  'multiplet'))
    adata.obs['sgrna_assignment'] = assignments
    return adata

Escaper Cell Filtering (Mixscape)

Why this matters: Not all sgRNA-positive cells actually edit. The escaper fraction is guide- and gene-dependent: Papalexi 2021 measured ~25% escapers for IFNGR2, perturbation rates of 39-92% across four IRF1 guides (i.e. 8-61% escapers), and no detectable perturbation at all for 15 genes. Including escapers dilutes the perturbation effect; Mixscape identifies and filters them.

Mixscape algorithm: For each perturbed cell, compute a "perturbation signature" = (its expression) - (mean of K nearest non-targeting-control cells). This signature isolates the perturbation effect from cell-state variation. Cells with perturbation signature similar to NTC distribution are escapers.

import pertpy as pt
import scanpy as sc

# adata is a scRNA-seq AnnData with 'sgrna_assignment' column
# Pertpy 0.6+ Mixscape API (verify against installed pertpy with help(pt.tl.Mixscape))
mixscape = pt.tl.Mixscape()
mixscape.perturbation_signature(
    adata=adata,
    pert_key='sgrna_assignment',     # .obs column with sgRNA target per cell
    control='NTC',                   # name of non-targeting control in pert_key column
    n_neighbors=20,                  # K neighbors for KNN-NTC subtraction
)
# Writes .layers['X_pert'] with perturbation-signature-corrected expression

# Filter escapers: classify perturbed cells as KO (true perturbation) or NP (non-perturbed/escaper)
mixscape.mixscape(
    adata=adata,
    pert_key='sgrna_assignment',     # pert_key (not 'labels' in modern pertpy)
    control='NTC',
    new_class_name='mixscape_class', # .obs column to write
)
# Defaults to layer='X_pert' (output of perturbation_signature)

# Keep only KO cells for downstream analysis
adata_ko = adata[adata.obs['mixscape_class_global'].isin(['KO'])   # mixscape_class holds '<gene> KO'; the bare label is in mixscape_class_global].copy()
print(f'KO cells: {adata_ko.n_obs} ({adata_ko.n_obs/adata.n_obs:.1%} of perturbed)')

Critical: Mixscape can fail when the perturbation has weak phenotype; empirically Mixscape detects perturbations with log-fold-change <-0.5 (depletion) reliably, but weaker effects collapse into the NTC distribution. For genome-wide screens, run Mixscape per perturbation; for low-effect perturbations, trust the assignment without filtering.

SCEPTRE for Low-MOI Differential Expression

Why this matters: Standard differential-expression tools (DESeq2, MAST) assume Gaussian-mixture distribution and fail at single-cell scale with sparse, zero-inflated data. SCEPTRE (Katsevich Lab, 2021; low-MOI variant Barry 2024 Genome Biol) uses a negative-binomial GLM with conditional resampling:

1. Per gene, fit NB GLM: log(expr_g) ~ pert_indicator + technical_factors
2. Compute z-score for the perturbation coefficient
3. Resample the pert_indicator (conditional on counts) 500-1000 times; compute permutation null
4. Get FDR via permutation; not parametric

library(sceptre)

# Input: sce object or sparse matrix + metadata
# Required: gene_expression_matrix, perturbation_indicator (binary per cell per pert),
#           technical_factors (batch, n_genes, etc.)

# For each gene + perturbation pair:
# Current sceptre API is a pipeline of composable steps:
sceptre_object <- import_data(response_matrix, grna_matrix, grna_target_data_frame,
                              moi = 'low', extra_covariates = covariates_df)
sceptre_object <- set_analysis_parameters(sceptre_object, discovery_pairs = pairs_df)
sceptre_object <- assign_grnas(sceptre_object)
sceptre_object <- run_qc(sceptre_object)
sceptre_object <- run_calibration_check(sceptre_object)
sceptre_object <- run_discovery_analysis(sceptre_object)
results <- get_result(sceptre_object, analysis = 'run_discovery_analysis')
# Output: per-gene-per-pert p-value, log-fold-change, FDR

Advantage over MAST: SCEPTRE's permutation NB GLM is the only method that maintains calibrated FDR in pooled-screen scRNA-seq (Barry 2024 benchmark). MAST and Wilcoxon are over-confident due to data sparsity.

Pertpy Unified Framework

Pertpy (https://pertpy.readthedocs.io) integrates Mixscape, distance-based perturbation comparison, EdgeR/PyDESeq2/WilcoxonTest DE, and factor models in a single AnnData-based interface. For SCEPTRE specifically, invoke the R sceptre package separately (Pertpy does not wrap it).

import pertpy as pt
import scanpy as sc

# Load data
mdata = pt.dt.papalexi_2021()      # returns a MuData object
adata = mdata['rna']  # built-in example from Mixscape paper

# Standard scRNA-seq preprocessing (scanpy)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)

# Mixscape escaper filtering (writes .layers['X_pert'] and .obs['mixscape_class'])
ms = pt.tl.Mixscape()
ms.perturbation_signature(adata, pert_key='perturbation', control='NT', n_neighbors=20)
ms.mixscape(adata, pert_key='perturbation', control='NT')

# Filter to KO cells
adata_ko = adata[adata.obs['mixscape_class'].isin(['KO', 'NT'])].copy()

# Pseudobulk differential expression via pertpy (PyDESeq2 backend)
de = pt.tl.PyDESeq2(adata_ko, design='~perturbation')
de.fit()
results_df = de.test_contrasts(contrast=('perturbation', 'GENE_X', 'NT'))

# For calibrated SCEPTRE on low-MOI single-cell data, use R sceptre directly
# (Barry 2024 Genome Biol; not bundled in pertpy)

Genome-Wide Perturb-Seq (Replogle 2022)

Replogle 2022 Cell 185:2559 demonstrated genome-wide Perturb-seq:

• 2.5M cells total; the genome-scale K562 screen targeted ~9,866 expressed genes (with a 2,057-gene essential subset)

• CRISPRi via dCas9-KRAB
• Native 10X 3' direct-capture for sgRNA
• Median >100 cells per perturbation as screened (Replogle 2022)
• Cluster-based analysis of perturbed cells reveals gene-program organization

Scaling principles:

• Cells per perturbation: 500-1,000 minimum for stable DE
• 10X channels: 10-30 channels at 5,000-10,000 cells each
• Cost: ~$50-100K for genome-scale

# Replogle-style genome-wide design
# Each cell -> 1 library element (low MOI)
# Each gene -> 1 dual-sgRNA CRISPRi element (2 distinct sgRNAs per element)
# Replogle 2022 retained >2.5M cells at a median >100 cells per perturbation
# Total: ~9,900 expressed genes x 1 element = ~9,900 elements

Factor-Based Analysis

For complex perturbation responses, decompose the per-cell perturbation effect into shared latent factors:

import pertpy as pt

# FR-Perturb ("Factorize-Recover") decomposes perturbation effects into shared factors.
# It is NOT part of pertpy: it is a standalone CLI from douglasyao/FR-Perturb
# (Yao et al. 2023 Nat Biotechnol). Run it outside Python:
#   python run_FR_Perturb.py --input <expression> --perturbations <matrix> --out <prefix>

Multiomic Perturb-seq (RNA + ATAC)

For chromatin readout: Use 10X Multiome with CRISPRi/a; sgRNA assignment via the same scATAC-seq library.

# Multiome: scRNA + scATAC + sgRNA
# Use ArchR or Signac for ATAC integration
# Use Pertpy for RNA-side DE
import muon as mu
mdata = mu.MuData({'rna': adata_rna, 'atac': adata_atac})
# Joint differential analysis across modalities