Pooled CRISPR Screens

Analyse pooled CRISPR screens with scRNA-seq readout.

Overview

Problem. Perturb many genes; read each effect transcriptomically.

Learning goals

• Functional genomics: perturb and observe
• Assigning guides to cells is the hard part

Figures

Tutorial

Analyze pooled CRISPR screens with single-cell RNA-seq readout using a tiered workflow: fast screening → target validation → rigorous differential expression.

When to Use This Skill

Use this skill when you have: - ✅ Pooled CRISPR screens with scRNA-seq (Perturb-seq, CROP-seq, CRISPRi/a) - ✅ 10X Feature Barcoding data (sgRNA captured as feature barcodes) - ✅ Multi-library experiments with biological replicates - ✅ sgRNA-to-cell mapping files (already assigned)

Quick Start (Example Data)

Test this skill with a real CRISPRi Perturb-seq dataset (~10 minutes):

from scripts.load_example_data import load_example_data
data = load_example_data()  # Downloads Papalexi 2021 (~140MB, cached after first run)

adata_list = data['adata_list']  # List of AnnData objects (one per batch)
mapping_files = data['mapping_files']  # sgRNA mapping files

For offline testing: Use load_example_data(dataset='demo') for a small synthetic dataset.

For your own data: Replace with your 10X feature-barcode matrices and sgRNA mapping files (see Inputs).

Installation

Core packages (required):

# Create conda environment
conda create -n crispr-screen python=3.8
conda activate crispr-screen

# Install packages
pip install scanpy==1.9+ anndata==0.8+ pandas numpy scipy
pip install scikit-learn  # For outlier detection
pip install diffxpy  # For differential expression

Visualization packages (required):

pip install matplotlib seaborn

Optional packages:

# For glmGamPoi (requires R ≥4.0)
# Install R, then:
# install.packages("BiocManager")
# BiocManager::install("glmGamPoi")
pip install rpy2  # Python-R interface

# For PDF report generation (optional, recommended)
pip install reportlab

Version requirements:

Inputs

Alternative inputs: CROP-seq pipeline output, custom sgRNA assignment files

See references/crispr_screen_best_practices.md for data format details.

Outputs

Clarification Questions

🚨 ALWAYS ask Question 1 FIRST. Do not ask about screen type, library details, or analysis goals before the user has answered Question 1.

1. Input Files (ASK THIS FIRST):

• Do you have 10X Feature-Barcode matrix files (.h5) to analyze?
  • If uploaded: How many libraries/replicates?
  • Expected format: raw_feature_bc_matrix.h5 from CellRanger
• Do you have sgRNA-to-cell mapping files?
  • If yes: Format? (TSV with barcode-sgRNA pairs). Doublets already filtered?
• Or use example data for testing? (Papalexi 2021 CRISPRi dataset — real data, ~140MB download)

🚨 IF EXAMPLE DATA SELECTED: All parameters are pre-defined (CRISPRi, THP-1 immune cells, 25 target genes). DO NOT ask questions 2-6. Proceed directly to Step 1.

Questions 2-6 are ONLY for users providing their own data:

2. Screen Type: CRISPRi (knockdown), CRISPRa (activation), CRISPRko (knockout), or other?

3. Library Information: How many biological replicates? Non-targeting controls present? Approximate perturbations (10-100, 100-1000, 1000+)? sgRNAs per gene (1-2, 3-5, 5+)?

4. Cell Type: Primary neurons, iPSC-derived, immune cells, cancer cell lines, other? Expected phenotype (cell death, differentiation, activation, subtle)?

5. sgRNA Naming: Do sgRNA IDs include gene names (e.g., "GENE_sgRNA1")? If not, need gene lookup table?

6. Analysis Goals (skip for example data):

• Primary goal?
  • a) Identify all hits with transcriptional phenotypes (standard)
  • b) Compare specific perturbations against each other
  • c) Focus on a known gene set
• Run downstream pathway enrichment?
  • a) Yes (recommended)
  • b) No, just hit identification

Typical Complete Workflow

This skill performs core Perturb-seq analysis: screening, validation, and hit identification. For a complete CRISPR screen workflow:

1. This skill: Screen analysis → validated hits with transcriptional phenotypes (H5AD + hit lists)
2. functional-enrichment-from-degs: Pathway enrichment on hit perturbations → biological interpretation
3. coexpression-network: Network analysis from perturbed cells → regulatory relationships
4. bulk-omics-clustering: Cluster perturbations by signature → functional groups

Quick start: "Analyze my pooled CRISPR screen and identify hits with pathway enrichment"

Why separate skills? Modular design works across screen types (CRISPRi/a/ko) and readouts (scRNA-seq, bulk). See Suggested Next Steps for details.

Standard Workflow

🚨 MANDATORY: USE SCRIPTS EXACTLY AS SHOWN - DO NOT WRITE INLINE CODE 🚨

Step 1 - Load data and map sgRNAs:

For example data:

from scripts.load_example_data import load_example_data
data = load_example_data()  # Papalexi 2021 CRISPRi (~140MB, cached)
adata_list = data['adata_list']
mapping_files = data['mapping_files']

For your own data:

from scripts.load_10x_libraries import load_multiple_libraries
from scripts.map_sgrna_to_cells import map_sgrna_to_adata
from scripts.qc_filtering import apply_qc_filters
from scripts.concatenate_libraries import concatenate_libraries

# Load, map, filter, concatenate
adata_list = load_multiple_libraries(["lib1.h5", "lib2.h5"])
adata_mapped = [map_sgrna_to_adata(ad, mf, sgrna_delimiter="_")
                for ad, mf in zip(adata_list, mapping_files)]
adata_filtered = [apply_qc_filters(ad, min_genes=2000, min_counts=8000, max_mito_pct=0.15)
                  for ad in adata_mapped]
adata = concatenate_libraries(adata_filtered, batch_labels=["lib1", "lib2"])

Step 2 - Preprocess and normalize:

from scripts.gene_name_corrections import correct_gene_names
from scripts.expression_filtering import filter_by_expression
from scripts.normalize_and_scale import normalize_only

# Preprocessing pipeline
adata = correct_gene_names(adata, corrections={}, gene_col='gene')
adata = filter_by_expression(adata, group_key='sgRNA', min_mean_expression=0.5)
adata = normalize_only(adata, target_sum=1e6, exclude_highly_expressed=True)

Step 3 - Run tiered analysis:

from scripts.screen_all_perturbations import screen_all_perturbations, call_hits
from scripts.validate_perturbations import validate_target_effect
from scripts.differential_expression_glmgampoi import run_glmgampoi_batch

# Fast screening → preliminary hits → validation
de_results = screen_all_perturbations(adata, control_group='non-targeting',
                                       gene_col='gene', test_method='t-test',
                                       output_dir='results/initial_screening/')
preliminary_hits = call_hits(de_results, min_de_genes=10)
validation = validate_target_effect(de_results, expected_direction='down',  # 'up' for CRISPRa
                                     min_log2fc=-0.5)
validated_hits = preliminary_hits[
    preliminary_hits['perturbation'].isin(
        validation[validation['target_affected']]['perturbation']
    )
]

# Rigorous DE for top hits (optional, requires R + glmGamPoi)
top_hits = validated_hits.head(50)
glmgampoi_results = run_glmgampoi_batch(adata, perturbations=top_hits['perturbation'].tolist(),
                                         donor_col='batch', output_dir='results/glmgampoi/')

Step 4 - Visualize and export:

from scripts.visualize_perturbations import create_volcano_plots
from scripts.export_results import export_all_results

# Volcano plots for top hits
for gene in top_hits['perturbation'].head(10):
    create_volcano_plots(de_results[gene], perturbation=gene, output_dir='figures/')

# Export all results (CSV, H5AD, pickle, PDF report, hit lists)
export_all_results(adata=adata, perturbation_summary=validated_hits,
                   de_results=de_results, output_dir='results/')

❌ IF YOU DON'T SEE VERIFICATION MESSAGES: You wrote inline code instead of using the functions. Stop and use the provided functions.

NEVER skip directly to writing inline code without trying the script first.

📁 QC Thresholds by Cell Type:

See references/qc_guidelines.md for detailed thresholds.

Common Issues

See references/troubleshooting_guide.md for detailed solutions.

Suggested Next Steps

Related Skills

• functional-enrichment-from-degs - Pathway enrichment on DE genes
• de-results-to-plots - Advanced visualization (heatmaps, custom plots)
• coexpression-network - Gene regulatory network inference
• bulk-omics-clustering - Cluster samples/perturbations by expression

References

See references/crispr_screen_best_practices.md for comprehensive guidelines.

Code preview

scripts/concatenate_libraries.py

"""
Concatenate Multiple Libraries

Merge filtered libraries into a single AnnData object with batch labels.
"""

import scanpy as sc
import anndata as ad
from typing import List, Optional
import warnings


def concatenate_libraries(
    adata_list: List[ad.AnnData],
    batch_labels: Optional[List[str]] = None,
    batch_key: str = 'batch'
) -> ad.AnnData:
    """
    Concatenate multiple AnnData objects (libraries) into one.

    Parameters
    ----------
    adata_list : list of AnnData
        List of filtered AnnData objects to concatenate
    batch_labels : list of str, optional
        Labels for each batch (e.g., ["lib1", "lib2", "lib3"])
        If None, uses numeric labels (0, 1, 2, ...)
    batch_key : str, default='batch'
        Column name for batch labels in obs

    Returns
    -------
    adata_combined : AnnData
        Concatenated AnnData object with batch labels

    Example
    -------
    >>> adata_combined = concatenate_libraries(
    ...     [adata1_filtered, adata2_filtered, adata3_filtered, adata4_filtered],
    ...     batch_labels=["lib1", "lib2", "lib3", "lib4"]
    ... )
    """
    print(f"Concatenating {len(adata_list)} libraries...")

    # Print summary of each library
    for i, adata in enumerate(adata_list):
        label = batch_labels[i] if batch_labels else f"batch_{i}"
        print(f"  {label}: {adata.n_obs} cells x {adata.n_vars} genes")

        # Check for required columns
        if 'gene' in adata.obs.columns:
            n_genes = adata.obs['gene'].nunique()
            print(f"    Unique perturbations: {n_genes}")

    # Concatenate using modern anndata.concat (not deprecated concatenate method)
    if batch_labels:
        # Add batch labels to each AnnData
        for i, adata in enumerate(adata_list):
            adata.obs[batch_key] = batch_labels[i]

        adata_combined = ad.concat(
            adata_list,
            join='outer',
            label=None,  # We already added batch labels above
            keys=None,
            index_unique=None
        )
    else:
        # Simple concatenation without batch labels
        adata_combined = ad.concat(
            adata_list,
            join='outer'
        )

    print(f"\nCombined dataset:")
    print(f"  Total cells: {adata_combined.n_obs}")
    print(f"  Total genes: {adata_combined.n_vars}")

    if 'gene' in adata_combined.obs.columns:
        n_unique_genes = adata_combined.obs['gene'].nunique()

scripts/detect_perturbed_cells.py

"""
Detect Perturbed Cells Using Outlier Detection

For each perturbation, identify cells with significant phenotypic changes using
outlier detection methods in PCA space computed on differentially expressed genes.
"""

import numpy as np
import pandas as pd
import scanpy as sc
import anndata as ad
from sklearn.svm import OneClassSVM
from sklearn.neighbors import LocalOutlierFactor
from sklearn.ensemble import IsolationForest
import diffxpy.api as de
from typing import Dict, List, Optional, Literal


def detect_perturbed_cells(
    adata: ad.AnnData,
    control_group: str = 'non-targeting',
    gene_col: str = 'gene',
    sgrna_col: str = 'sgRNA',
    n_pcs: int = 4,
    method: Literal['LocalOutlierFactor', 'IsolationForest', 'OneClassSVM'] = 'LocalOutlierFactor',
    min_de_genes: int = 10,
    min_outliers: int = 10,
    pval_threshold: float = 0.05,
    save_plots: bool = True,
    plot_dir: str = 'figures/'
) -> Dict:
    """
    Detect perturbed cells with significant phenotypes using outlier detection.

    For each perturbation vs control:
    1. Run differential expression (t-test)
    2. Select top DE genes (p < pval_threshold)
    3. If ≥min_de_genes found, compute PCA on DE genes
    4. Train outlier detector on control cells
    5. Predict outliers in perturbed cells
    6. If ≥min_outliers detected, flag as "hit"

    Parameters
    ----------
    adata : AnnData
        Log-normalized AnnData (not scaled)
    control_group : str, default='non-targeting'
        Label for non-targeting control
    gene_col : str, default='gene'
        Column in obs with gene/perturbation labels
    sgrna_col : str, default='sgRNA'
        Column in obs with sgRNA labels
    n_pcs : int, default=4
        Number of PCs to compute on DE genes
    method : str, default='LocalOutlierFactor'
        Outlier detection method: 'LocalOutlierFactor', 'IsolationForest', 'OneClassSVM'
    min_de_genes : int, default=10
        Minimum DE genes required to proceed with outlier detection
    min_outliers : int, default=10
        Minimum outliers to call perturbation as "hit"
    pval_threshold : float, default=0.05
        P-value threshold for calling DE genes
    save_plots : bool, default=True
        Generate UMAP/PCA plots for each perturbation
    plot_dir : str, default='figures/'
        Directory to save plots

    Returns
    -------
    results : dict
        Dictionary containing:
        - 'outlier_cells': list of cell barcodes flagged as outliers
        - 'perturbation_summary': DataFrame with hit calls per perturbation
        - 'de_results': dict of DE results per perturbation
        - 'adata_dict': dict of per-perturbation AnnData objects with classifications

    Example
    -------
    >>> results = detect_perturbed_cells(
    ...     adata,

scripts/differential_expression.py

"""
Differential Expression Analysis per Perturbation

For each perturbation, perform differential expression analysis comparing
perturbed cells to non-targeting controls.
"""

import os
import pandas as pd
import anndata as ad
import diffxpy.api as de
from typing import Dict, Optional, Literal


def run_de_analysis(
    adata: ad.AnnData,
    control_group: str = 'non-targeting',
    gene_col: str = 'gene',
    test_method: Literal['t-test', 'wilcoxon', 'negative_binomial'] = 't-test',
    use_outliers_only: bool = True,
    outlier_cells: Optional[list] = None,
    top_n_genes: int = 50,
    output_dir: str = 'DEG/rank_test/',
    is_logged: bool = True
) -> Dict[str, pd.DataFrame]:
    """
    Run differential expression analysis for each perturbation vs control.

    Parameters
    ----------
    adata : AnnData
        Log-normalized AnnData (use adata.raw or normalized data)
    control_group : str, default='non-targeting'
        Label for control group
    gene_col : str, default='gene'
        Column in obs with perturbation labels
    test_method : str, default='t-test'
        Statistical test: 't-test', 'wilcoxon', 'negative_binomial'
    use_outliers_only : bool, default=True
        If True and outlier_cells provided, use only outlier cells for DE
    outlier_cells : list, optional
        List of cell barcodes identified as outliers (from detect_perturbed_cells)
    top_n_genes : int, default=50
        Number of top genes to highlight in output
    output_dir : str, default='DEG/rank_test/'
        Directory to save per-perturbation DE results
    is_logged : bool, default=True
        Whether data is log-transformed

    Returns
    -------
    de_results : dict
        Dictionary mapping gene -> DE results DataFrame

    Example
    -------
    >>> de_results = run_de_analysis(
    ...     adata,
    ...     control_group='non-targeting',
    ...     test_method='t-test',
    ...     use_outliers_only=True,
    ...     outlier_cells=results['outlier_cells'],
    ...     top_n_genes=50
    ... )
    """
    os.makedirs(output_dir, exist_ok=True)

    # Subset to outliers + controls if requested
    if use_outliers_only and outlier_cells is not None:
        print(f"Subsetting to outlier cells + controls...")
        control_cell_barcodes = adata[adata.obs[gene_col] == control_group].obs_names.tolist()
        kept_cells = list(set(outlier_cells + control_cell_barcodes))
        adata_selected = adata[adata.obs_names.isin(kept_cells)].copy()
        print(f"  Using {len(outlier_cells)} outlier cells + {len(control_cell_barcodes)} control cells")
    else:
        adata_selected = adata.copy()
        print("Using all cells for DE analysis")

    # Get unique perturbations
    genes = adata_selected.obs[gene_col].unique().tolist()

Companion files