bio-crispr-screens-hit-calling

## Version Compatibility Reference examples tested with: MAGeCK 0.5.9+, BAGEL2 2.0, drugZ Aug 2019+, JACKS 0.2.0+, Chronos 2.0+ (DepMap), CERES 1.0+, pandas 2.2+, numpy 1.26+, scipy 1.12+, statsmodels 0.14+. Before using code patterns, verify installed versions match. If versions differ: - CLI: `mageck --version`, `BAGEL.py version`, `python drugz.py --help` - Python: `pip show crispr_chronos` (JACKS installs from GitHub, not PyPI) If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying. ## Hit Calling Decision Tree **"Identify significant hits in my CRISPR screen"** -> Choose the analysis method that matches the experimental design, statistical assumptions, and quality grade of the screen. Reconcile across methods when high-stakes hits must be validated. The primary hit-calling methods cover non-overlapping niches; the decision is not "which is best" but "which matches the design." | Design / question | Primary method | Why | Secondary check | |--------------------|----------------|-----|------------------| | Two-condition essentiality, one cell line, no CN concerns | MAGeCK RRA | Robust, fast, gold-standard for ranked analysis | BAGEL2 (Bayes factor on same data) | | Time course (3+ timepoints) | MAGeCK MLE | RRA cannot model multi-condition | JACKS (efficacy-aware) | | Multi-cell-line panel (cancer dependency) | Chronos | Models CN bias + screen quality jointly | MAGeCK MLE per line + meta-analysis | | Drug screen (vehicle vs drug) | drugZ | Bidirectional Z; vehicle-anchored | MAGeCK MLE with dose covariate | | Multi-screen joint, same library | JACKS | Shared efficacy; enables ~2.5x smaller screens | MAGeCK MLE; results should converge | | Essentiality classification with reference sets | BAGEL2 | Bayes factor with CEGv2/NEGv1 calibration | MAGeCK RRA | | Combinatorial / paired guide | MAGeCK MLE with GI scoring | Models interaction term; see [[combinatorial-screens]] | Custom GI scoring | | Single-cell perturbation (Perturb-seq) | SCEPTRE | NB GLM + permutation; see [[perturb-seq-analysis]] | Mixscape pre-filter | | Cancer-line copy-number screen | Chronos (preferred) or CERES | Joint CN-bias + gene-effect modeling; see [[copy-number-correction]] | CRISPRcleanR pre-hoc + MAGeCK | ## Statistical Models Compared | Method | Year | Statistical model | Tests | Best for | Fails when | |--------|------|-------------------|-------|----------|------------| | MAGeCK RRA | 2014 | NB per-sgRNA -> alpha-RRA per gene | Two-sided | General two-condition | >40% guides change (median norm breaks); time course; cancer-line CN | | MAGeCK MLE | 2015 | NB GLM with design matrix; per-gene beta | Wald per condition | Multi-condition / time course | Cell-line specific essentiality; CN bias | | BAGEL2 | 2021 | Bayes factor from log-likelihood ratio | Essential vs non-essential | Essentiality classification | Non-essentiality screens; drug screens | | drugZ | 2019 | Bidirectional Z-score on guide-level LFC | Sensitizer vs suppressor | Drug-modifier / chemogenomic | Essentiality (no biological prior); time-course | | JACKS | 2019 | Variational Bayes: LFC = gene * efficacy | Per-gene posterior | Multi-screen joint, library calibration | Single screen; cross-chemistry | | Chronos | 2021 | Cell-population dynamics ODE + NB | Gene effect adjusted for screen quality | Cancer-line panels, longitudinal | Single screen; non-cancer applications | | CERES | 2017 | Nonlinear model decoupling CN-bias from gene effect | Per-gene effect | Cancer-line panel with CN profile | Superseded by Chronos at DepMap | ## RRA vs MLE Within MAGeCK | Property | RRA (`mageck test`) | MLE (`mageck mle`) | |----------|----------------------|---------------------| | Conditions supported | 2 | Multiple (design matrix) | | Statistical test | Robust rank aggregation | Wald on beta from NB GLM | | Output | neg/pos score, FDR per direction | beta per condition | | sgRNA efficiency | Not modeled (optional fixed input) | Modeled via `--sgrna-efficiency` | | Outlier robustness | High (rank-based) | Lower (likelihood-based) | | Best for | Standard 2-condition screen | Time course, drug screen, multi-cell-line, paired | | Speed | Fast | Slow (per-gene optimization) | ## Algorithmic Taxonomy: Why Each Was Built | Method | Designed to solve | |--------|--------------------| | MAGeCK RRA | First robust statistical framework for CRISPR-screen ranking; alpha-RRA borrowed from RRA in microarray meta-analysis | | MAGeCK MLE | Extend MAGeCK to multi-condition; explicit beta scores allow direct LFC interpretation | | BAGEL2 | Reference-set-anchored Bayesian classification; precision-recall calibrated; tumor-suppressor sensitivity (BAGEL1 was uni-directional) | | drugZ | Drug-modifier screens have low effect sizes and need bidirectional sensitivity; STARS/MAGeCK miss synthetic-lethal hits | | JACKS | Sample-size reduction via library-shared efficacy; library calibration as side product | | Chronos | DepMap-scale (1000+ cell lines, billions of cell-divisions) needs population-dynamics model; CN bias + screen quality first-class | | CERES | First to formally decouple CN from gene effect at DepMap scale; superseded but historically important | ## Run All Five on the Same Data (Consensus Strategy) **Goal:** For high-stakes hits (drug-target nomination, paper-level claims), require agreement across 2-3 orthogonal methods. **Approach:** Run MAGeCK + BAGEL2 + (drugZ or JACKS) on the same count matrix; rank by each; classify hits as called by 1, 2, or 3 methods. ```python import pandas as pd def consensus_hits(mageck_path, bagel_path, drugz_path, mageck_fdr_thresh=0.05, bagel_bf_thresh=5, drugz_fdr_thresh=0.05): '''Build consensus across MAGeCK / BAGEL2 / drugZ on the same screen. Each hit gets a count of supporting methods.''' mageck = pd.read_csv(mageck_path, sep='\t')[['id', 'neg|fdr']].rename(columns={'id': 'gene', 'neg|fdr': 'mageck_neg_fdr'}) bagel = pd.read_csv(bagel_path, sep='\t')[['GENE', 'BF']].rename(columns={'GENE': 'gene', 'BF': 'bagel_bf'}) drugz = pd.read_csv(drugz_path, sep='\t')[['GENE', 'fdr_synth']].rename(columns={'GENE': 'gene', 'fdr_synth': 'drugz_synth_fdr'}) merged = mageck.merge(bagel, on='gene', how='outer').merge(drugz, on='gene', how='outer') merged['mageck_hit'] = merged['mageck_neg_fdr'] < mageck_fdr_thresh merged['bagel_hit'] = merged['bagel_bf'] > bagel_bf_thresh merged['drugz_hit'] = merged['drugz_synth_fdr'] < drugz_fdr_thresh merged['consensus_count'] = (merged[['mageck_hit', 'bagel_hit', 'drugz_hit']].astype(int)).sum(axis=1) return merged.sort_values('consensus_count', ascending=False) ``` **Confidence tiers:** | Tier | Definition | Validation requirement | |------|------------|-------------------------| | Tier 1 (high) | Called by 3/3 methods | Arrayed validation; orthogonal modality (CRISPRi if originally Cas9) | | Tier 2 (medium) | Called by 2/3 methods | Arrayed validation in matched line | | Tier 3 (exploratory) | Called by 1/3 methods | Treat as hypothesis; further screens before publication | ## Reconciliation: When Two Methods Disagree | Pattern | Likely cause | Action | |---------|--------------|--------| | MAGeCK significant, BAGEL2 not | BAGEL2 trained on CEGv2/NEGv1; gene is essential but not in reference | Trust MAGeCK; flag for follow-up | | BAGEL2 significant, MAGeCK not | BAGEL2 has tumor-suppressor sensitivity MAGeCK lacks | Investigate sgrna_summary for one weak guide | | MAGeCK significant, JACKS not | JACKS down-weighted one outlier guide | Trust JACKS if guides agree; outlier may be off-target | | Chronos and MAGeCK disagree on cancer line | Chronos accounts for CN; MAGeCK does not | Trust Chronos; apply [[copy-number-correction]] | | drugZ significant, MAGeCK not on drug screen | drugZ bidirectional Z is more sensitive | Trust drugZ for chemogenomic; MAGeCK may miss small effects | | MAGeCK MLE significant, MAGeCK RRA not in 2-condition | Beta-score effect size is significant but rank-based not | Trust MLE if guides consistent; RRA may be over-conservative | | All methods disagree | Either no real biology or all methods are mis-applied | Stop. Re-audit QC; check chemistry / library / design matrix | ## Second-Best sgRNA Conservative Rule **Goal:** Reduce false positives from single outlier sgRNAs by requiring the second-most-extreme guide per gene to also be a hit. **Approach:** For each gene, sort sgRNAs by LFC; require the second-best LFC to exceed a threshold. Rejects genes that depend on one extreme guide. ```python def second_best_lfc(sgrna_lfc_df, genes_series, direction='neg'): '''Return per-gene LFC of the second-best sgRNA in the direction of interest. For dropout (direction="neg"), second-most-negative LFC.''' results = [] for gene in genes_series.unique(): gene_lfc = sgrna_lfc_df[genes_series == gene].sort_values() if direction == 'neg': second = gene_lfc.iloc[1] if len(gene_lfc) >= 2 else gene_lfc.iloc[0] else: second = gene_lfc.iloc[-2] if len(gene_lfc) >= 2 else gene_lfc.iloc[-1] results.append({'gene': gene, 'second_best_lfc': second}) return pd.DataFrame(results) ``` **Rule:** A high-confidence hit has second-best LFC also passing the threshold. A guide-of-one hit has only one extreme guide and should be flagged for orthogonal validation. This rule predates JACKS and is implicit in MAGeCK RRA but explicit elsewhere. ## Multiple-Testing Correction Conventions | Method | Native correction | Cross-method comparison | |--------|---------------------|---------------------------| | MAGeCK RRA | BH per direction | `neg|fdr`, `pos|fdr` | | MAGeCK MLE | BH per condition | `<cond>|fdr` | | BAGEL2 | Bootstrap BF; reports BF threshold | BF > 6 ≈ 90% posterior (Hart 2017); ~5% FDR by convention | | drugZ | BH per direction | `fdr_synth`, `fdr_supp` | | JACKS | Posterior probability + BH | `fdr_log10` (log10 FDR) | | Chronos | DepMap gene-effect probability | `effect_probability` | **Reconciliation:** BF >6 in BAGEL2 corresponds to ~90% posterior probability (Hart 2017 G3, by overlap with CEGv2) and is commonly used as a stringent cutoff roughly comparable to MAGeCK FDR 0.05. Treat that equivalence as an approximate convention, not an exact calibration. drugZ FDR is per-direction; the `fdr_synth` and `fdr_supp` columns are independent BH corrections. ## Order of Operations ``` 1. Library design (see library-design) <- design quality dictates hit calling 2. Plasmid pool sequencing <- baseline; non-negotiable 3. Run screen at MOI 0.3, 500x coverage 4. Sequence endpoint 5. Run mageck count <- generates raw + normalized counts 6. Screen QC (see screen-qc) <- gates downstream method choice 7. Copy-number correction if cancer line <- CRISPRcleanR or Chronos; see copy-number-correction 8. Batch correction if multi-batch <- see batch-correction 9. Hit calling (this skill) <- choose method by design 10. Consensus across 2-3 methods <- for high-stakes hits 11. Orthogonal validation <- arrayed; different chemistry 12. Pathway analysis <- see pathway-analysis/gsea ``` ## Custom z-score Hit Calling (when standard tools don't fit) **Goal:** Compute gene-level z-scores when neither MAGeCK nor BAGEL2 fits the experimental design. **Approach:** RPM-normalize, compute per-sgRNA log2 fold-changes, aggregate to gene level, derive z-score from the null distribution of non-targeting controls (cleanest) or all genes (assumes <40% changing), apply BH correction. ```python import pandas as pd import numpy as np from scipy import stats from statsmodels.stats.multitest import multipletests def custom_zscore_hit_calling(counts_df, ctrl_cols, treat_cols, genes_series, ntc_genes=None): '''Z-score gene-level hit calling. If ntc_genes provided, null derived from NTCs only; otherwise from all genes (assumes <40% changing).''' def rpm(df): return df.div(df.sum(axis=0), axis=1) * 1e6 ctrl_rpm = rpm(counts_df[ctrl_cols]) treat_rpm = rpm(counts_df[treat_cols]) lfc_per_sgrna = np.log2((treat_rpm.mean(axis=1) + 1) / (ctrl_rpm.mean(axis=1) + 1)) gene_lfc = pd.DataFrame({'gene': genes_series, 'lfc': lfc_per_sgrna}).groupby('gene')['lfc'].agg(['mean', 'std', 'count']) gene_lfc.columns = ['mean_lfc', 'std_lfc', 'n_sgrnas'] if ntc_genes is not None: null = gene_lfc.loc[gene_lfc.index.isin(ntc_genes), 'mean_lfc'] null_mean, null_std = null.median(), null.std() else: null_mean = gene_lfc['mean_lfc'].median() null_std = gene_lfc['mean_lfc'].std() gene_lfc['z'] = (gene_lfc['mean_lfc'] - null_mean) / null_std gene_lfc['p'] = 2 * stats.norm.sf(np.abs(gene_lfc['z'])) gene_lfc['fdr'] = multipletests(gene_lfc['p'], method='fdr_bh')[1] return gene_lfc.sort_values('z') ```