bio-crispr-screens-perturb-seq-analysis

What this skill does

## 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

| Method | Year | Architecture | Readout | MOI | Single-cell sgRNA detection |
|--------|------|--------------|---------|-----|------------------------------|
| Perturb-seq (Dixit 2016, *Cell*) | 2016 | sgRNA expressed in cassette; direct PCR capture | scRNA-seq | High (4-50 sgRNAs/cell) | Yes via amplicon-PCR pre-sequencing |
| CROP-seq (Datlinger 2017, *Nat Methods*) | 2017 | sgRNA expressed via tRNA-spacer-sgRNA-pgk-puro; barcoded in 3'UTR of fluorescent marker; captured in 3' droplet | scRNA-seq | Low (1-2 sgRNAs/cell) | Native via 10X 3' chemistry |
| Perturb-CITE-seq (Frangieh 2021, *Nat Genet*) | 2021 | Adds surface-protein hashtag oligos to CROP-seq | scRNA-seq + ADT (protein) | Low | CROP-seq architecture |
| ECCITE-seq (Mimitou 2019, *Nat Methods*) | 2019 | Surface-protein hashtag with sgRNA-marked cells | scRNA-seq + ADT | Low | Hash + sgRNA |
| Perturb-ATAC (Rubin 2019, *Cell*) | 2019 | scATAC-seq readout | scATAC | Low | sgRNA capture via separate library prep |
| Perturb-multiome (10X) | 2021+ | scRNA + scATAC simultaneously | scRNA + ATAC | Low | Direct capture from sgRNA cassette |
| Replogle GW Perturb-seq (2022, *Cell*) | 2022 | Multiplexed CRISPRi with sgRNA barcoding | scRNA-seq | 1 sgRNA/cell | Direct capture |

**Decision rule:** For genome-wide CRISPRi screens, Replogle's CRISPRi + 10X 3' direct-capture protocol is the gold standard (2.5M cells across 2,058 perturbations 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.

```python
# 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. 30-60% of cells may receive sgRNA but fail to undergo gene knockdown ("escapers"). Including escapers dilutes the perturbation effect; Mixscape (Papalexi 2021) 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.

```python
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'].isin(['KO'])].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

```r
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:
results <- run_sceptre_low_moi(
    expression_matrix    = gene_expr,            # genes x cells
    grouping_var         = pert_indicator,       # length(cells); 0 = NTC, 1 = perturbed
    technical_factors    = covariates_df,        # n_genes, n_umi, batch
    response_grouping_var = 'gene_id',
    n_permutations       = 1000
)
# 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 po