bio-workflows-crispr-screen-pipeline

## Version Compatibility Reference examples tested with: MAGeCK 0.5.9+, BAGEL2 1.0.5+, drugZ Aug 2019+, JACKS 0.2.0+, Chronos 2.0+, CRISPRcleanR 3.0+ (R), Pertpy 0.6+, PRIDICT2, CRISPResso2 2.2.14+, MAGeCKFlute 2.0+, pandas 2.2+, numpy 1.26+, matplotlib 3.8+. Before using code patterns, verify installed versions match. If versions differ: - CLI: `mageck --version`, `BAGEL.py fc --help`, `drugz -h`, `CRISPResso --version` - Python: `pip show pertpy scanpy anndata` (mageck-vispr via conda `mageck --version`; JACKS/Chronos are GitHub installs — check their repos) - R: `packageVersion('CRISPRcleanR')`, `packageVersion('MAGeCKFlute')` If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying. ## CRISPR Screen Pipeline **"Analyze my pooled or single-cell CRISPR screen end-to-end"** -> Pick the screen design branch, run guide counting, audit six QC stages, apply copy-number and batch correction as needed, run the design-matched hit-calling method, and consolidate across methods for high-confidence hits. This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step. ## The governing principle Every LFC, QC gate, and hit call is computed against a reference that is committed once at library-order time; a wrong-but-silent commitment invalidates the endpoint with no error thrown. 1. **The guide LIBRARY definition (guide->gene map + control classes) is the denominator, the calibrator, and the training reference — committed once.** The library must carry non-targeting controls (NTCs, ~1%, the null distribution) AND CEGv2 reference essentials + NEGv1 non-essentials (the positive/negative calibrators for PR-AUC and BAGEL2/Chronos priors). NTCs calibrate the null/FDR; CEGv2/NEGv1 calibrate PR-AUC — swapping or dropping a class silently breaks FDR or QC. 2. **The baseline choice has a right answer and rescales every hit.** Dropout/enrichment LFC is against a baseline: plasmid pool for the cloning-bottleneck baseline, Day-0/T0 for the biology baseline, and vehicle (NOT Day-0) for drug screens — drug-vs-Day-0 conflates drug effect with normal proliferation. 3. **Copy-number correction MUST precede hit calling in cancer cell lines.** Multiple simultaneous Cas9 cuts at amplified loci trigger a gene-independent DNA-damage/G2 arrest (Aguirre 2016; Munoz 2016; the effect appears in both TP53-mutant and TP53-wild-type lines, though Aguirre 2016 found TP53 status correlates with its magnitude -- separately, Ihry 2018 / Haapaniemi 2018 report p53-dependent toxicity of Cas9 cutting generally), so amplified regions look essential regardless of gene function; calling hits first yields false essentials at ERBB2/MYC/FGFR1. Run CRISPRcleanR/Chronos BEFORE hit calling, or use CRISPRi to bypass the DSB. This is a pipeline step, not a post-hoc interpretation. Verifying `abs(rho(LFC,CN)) < 0.1` afterwards needs a matched CN profile: CRISPRcleanR corrects unsupervised without one, so the check happens in crispr-screens/copy-number-correction (or use Chronos, which consumes CN directly). 4. **CEGv2 essential-gene depletion is the screen's built-in positive control.** If known essentials do not deplete (CEGv2 PR-AUC below ~0.7), the screen failed selection and NO novel hit is trustworthy regardless of its p-value — the seam analog of a spike-in. Normalize -> QC -> (CN correct) -> hit-call, never hit-call first; add batch as an MLE covariate, never pre-corrected with ComBat on counts (distorts the NB mean-variance the caller assumes). ## Made-once commitments | Commitment | Consequence inherited downstream | |------------|----------------------------------| | Guide library (guide->gene map + NTC/CEGv2/NEGv1 control classes) | The counting denominator, the QC calibrator, the hit-calling priors; a missing/misassigned class breaks FDR or PR-AUC | | Baseline (plasmid pool / Day-0 / vehicle) | Every LFC; drug-vs-Day-0 conflates drug effect with proliferation | | Screen type (dropout / enrichment / FACS / drug-modifier) | Which hit-calling method is even valid | | Copy-number profile (cancer lines) | Whether amplicon artifacts are removed before hit calling; residual rho(LFC,CN) is the tell | ## Pipeline Branches by Screen Design ``` Library Design ([[library-design]]) | v FASTQ Files -> mageck count -> count matrix | v Six-Stage QC ([[screen-qc]]) | +---------------------+---------------------+ | | v v Cancer cell line? Non-cancer? Apply CN correction No CN correction needed ([[copy-number-correction]]) | | +---------------------+---------------------+ v Multi-batch? Apply batch covariate ([[batch-correction]]) | v Pick hit-calling method by design ([[hit-calling]]) | +-----------+---------+---------+-----------+-----------+ | | | | | | v v v v v v 2-cond Time Drug Essential Multi- Specialized MAGeCK RRA MAGeCK drugZ BAGEL2 screen (PE/BE/SC/ MLE JACKS or in vivo/ Chronos combinat) | | | | | | +-----------+---------+---------+-----------+-----------+ v Tier-based consensus v Orthogonal validation ``` ## Step 1: Library Design and Pre-Screen Validation Reference [[library-design]] for full library composition. Verify before sequencing: - Plasmid pool Gini <0.1 (Li W et al 2015 MAGeCK-VISPR, Genome Biol 16:281) - >=99% guides detected at >25 reads/guide - Skew (p90/p10) <2 - NTCs comprise ~1% of library; CEGv2 reference essentials + NEGv1 non-essentials included ## Step 2: Guide Counting **Goal:** Turn raw FASTQ into a per-guide count matrix with consistent sample labels. **Approach:** Run mageck count with the library CSV, sample labels in column order, the vector adapter trimmed off the 5' end, and median normalization. ```bash mageck count \ --list-seq library.csv \ --sample-label Plasmid,Day0,Veh_r1,Veh_r2,Drug_r1,Drug_r2 \ --fastq Plasmid.fq.gz Day0.fq.gz Veh_r1.fq.gz Veh_r2.fq.gz Drug_r1.fq.gz Drug_r2.fq.gz \ --norm-method median \ --output-prefix experiment \ --trim-5 5 # integer base-count (or AUTO), NOT an adapter sequence; 5 trims the CACCG scaffold ``` For Cas12a libraries (Inzolia, in4mer): see [[combinatorial-screens]]. For 10X single-cell direct capture: use cellranger-arc or pertpy-aware counting; see [[perturb-seq-analysis]]. ## Step 3: Six-Stage Quality Control **Goal:** Decide whether the screen is analyzable before calling any hits, using six orthogonal QC stages. **Approach:** Load the count matrix, compute per-sample Gini, zero-fraction, and depth plus replicate correlation against the hard gates below. Essential-gene recovery (CEGv2 PR-AUC) is a separate check computed once endpoint-vs-baseline LFCs exist (it needs CEGv2/NEGv1 labels) -- see screen-qc. ```python import pandas as pd import numpy as np counts = pd.read_csv('experiment.count.txt', sep='\t', index_col=0) genes = counts['Gene'] count_matrix = counts.drop('Gene', axis=1) def gini(x): x = np.sort(x[x > 0].astype(float)) if x.size == 0: return np.nan n = x.size cumx = np.cumsum(x) return (n + 1 - 2 * np.sum(cumx) / cumx[-1]) / n per_sample = pd.DataFrame({ 'pct_zero': (count_matrix == 0).sum() / len(count_matrix) * 100, 'gini': count_matrix.apply(gini), 'reads_per_sgrna': count_matrix.sum() / len(count_matrix), }) log_counts = np.log10(count_matrix + 1) pearson = log_counts.corr() print(per_sample) print('Replicate Pearson:', pearson.values[pearson.values < 1].mean()) ``` Hard gates from [[screen-qc]]: - Plasmid Gini <0.1; endpoint <0.3 (or <0.55 for heavy drug screens) - Replicate Pearson on log-counts >0.85 - CEGv2 PR-AUC >0.7 against the Hart 2017 reference essential gene set (community convention, not a threshold defined in that paper) - Reads per sgRNA per sample >=300 (DepMap convention) ## Step 4: Copy-Number Correction (Cancer Cell Lines Only) If screening in a cancer cell line, apply CRISPRcleanR (unsupervised, no CN profile needed) or Chronos (joint with CN profile). Required to remove Aguirre 2016 / Munoz 2016 amplicon artifact. **Goal:** Strip the copy-number amplicon artifact that makes amplified regions look essential in cancer lines. **Approach:** Run CRISPRcleanR unsupervised genome-wide LFC correction (no CN profile needed), then feed the corrected counts downstream; for DepMap-scale panels with matched CN, use Chronos instead. ```r library(CRISPRcleanR) data(KY_Library_v1.0) norm <- ccr.NormfoldChanges('experiment.count.txt', min_reads = 30, EXPname = 'screen', libraryAnnotation = KY_Library_v1.0) # arg 1 is the file PATH gw_lfc <- ccr.logFCs2chromPos(norm$logFCs, KY_Library_v1.0) # $logFCs, not $norm_fold_changes cleaned <- ccr.GWclean(gw_lfc, display = TRUE, label = 'screen') corrected_counts <- ccr.correctCounts('screen', norm$norm_counts, cleaned, KY_Library_v1.0) # (CL, normalised_counts, correctedFCs, libraryAnnotation) # ccr.correctCounts returns an in-memory frame; it does NOT write this file. Persist it, because the # hit callers below read a count TABLE from disk -- the CN-correction commitment in rule 3 is only # honored if that file, not the raw experiment.count.txt, is what MAGeCK / BAGEL2 / drugZ consume. write.table(corrected_counts, 'screen_cleanr_corrected_counts.txt', sep = '\t', quote = FALSE, row.names = FALSE) ``` For DepMap-scale panels with longitudinal data + matched CN, use Chronos. See [[copy-number-correction]]. ## Step 5: Batch Correction (Multi-Batch Screens) For multi-batch screens, add batch as a covariate in MAGeCK MLE rather than pre-correcting with ComBat. See [[batch-correction]] for full decision tree. ## Step 6: Method-Matched Hit Calling **Goal:** Call hits with the method that matches the experimental design, plus at least one orthogonal method for consensus. **Approach:** Pick by design - RRA or BAGEL2 for two-condition essentiality, MLE for time course, drugZ for drug-modifier, JACKS for multi-screen, Chronos for cancer panels - and run two methods so the consensus step has something to reconcile. ### 6a. Two-condition essentiality (MAGeCK RRA or BAGEL2) Cancer cell lines: pass the CRISPRcleanR-corrected count file (`screen_cleanr_corrected_counts.txt` from Step 4) as `--count-table`/`-i` below, NOT the raw `experiment.count.txt` — the CN-correction commitment is only honored if the corrected counts are what the hit caller reads. The `Day0`/`Day14_r*` columns below illustrate a time-course dropout design; they must match the count step's `--sample-label` (the drug-screen count above uses `Plasmid,Day0,Veh_r*,Drug_r*`). ```bash mageck test \ --count-table experiment.count.txt \ --treatment-id Day14_r1,Day14_r2,Day14_r3 \ --control-id Day0 \ --norm-method median \ --output-prefix essentiality_rra ``` ```bash BAGEL.py fc -i experiment.count.txt -o experiment -c Day0 --min-reads 30 # -o is an output LABEL; fc writes experiment.foldchange BAGEL.py bf -i experiment.foldchange -o bayes_factor.txt -e CEGv2.txt -n NEGv1.txt \ -c Day14_r1,Day14_r2,Day14_r3 # add -b -NB 1000 for bootstrapping; -k is not a bf option ``` ### 6b. Time-course / multi-condition (MAGeCK MLE) ```bash mageck mle --count-table experiment.count.txt --design-matrix design.txt \ --output-prefix timecourse_mle --norm-method median ``` ### 6c. Drug-modifier (drugZ) ```bash python drugz.py \ -i experiment.count.txt \ -o drugz_output.txt \ -c Veh_r1,Veh_r2 \ -x Drug_r1,Drug_r2 \ -p 5 ``` drugZ requires vehicle as control, not Day-0. See [[drugz-chemogenomic]]. ### 6d. Multi-screen joint analysis (JACKS) ```bash python run_JACKS.py experiment.count.txt replicatemap.txt guidemap.txt \ --rep_hdr Replicate --sample_hdr Sample --ctrl_sample_hdr Control \ --sgrna_hdr sgRNA --gene_hdr Gene --outprefix jacks_out --apply_w_hp ``` ### 6e. Cancer cell-line panels (Chronos) ```python import chronos # All three inputs must be dicts of DataFrame keyed by library name, not bare DataFrames. model = chronos.Chronos(sequence_map={'screen': sequence_map}, guide_gene_map={'screen': guide_gene_map}, readcounts={'screen': counts_df}) # readcounts=, not reads= model.train(nepochs=301) # nepochs (default 301), not n_steps gene_effects = model.gene_effect # attribute, not a method call # Copy-number correction is a separate post-hoc step (chronos.alternate_CN(gene_effect, copy_number) / a CN matrix), not a constructor arg ``` DepMap quarterly standard; handles CN bias + screen quality + longitudinal jointly. ## Step 7: Tier-Based Consensus **Goal:** Consolidate the per-method calls into confidence tiers. **Approach:** Merge each method's gene-level result, threshold each to a per-method hit flag, and tier by how many methods agree (Tier 1 = all three, Tier 2 = two of three). ```python mageck = pd.read_csv('essentiality_rra.gene_summary.txt', sep='\t')[['id', 'neg|fdr']].rename( columns={'id': 'gene', 'neg|fdr': 'mageck_neg_fdr'}) bagel = pd.read_csv('bayes_factor.txt', sep='\t')[['GENE', 'BF']].rename( columns={'GENE': 'gene', 'BF': 'bagel_bf'}) drugz_df = pd.read_csv('drugz_output.txt', sep='\t')[['GENE', 'fdr_synth']].rename( columns={'GENE': 'gene', 'fdr_synth': 'drugz_synth_fdr'})