Hit calling

Version Compatibility

Reference examples tested with: MAGeCK 0.5.9+, BAGEL2 1.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, drugz --help
• Python: pip show jacks chronos-cn

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 five 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	Joint efficacy + 2.5x sample-size reduction	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.

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.

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
MAGeCK MLE	BH per condition	`<cond>
BAGEL2	Bootstrap BF; reports BF threshold	BF > 6 ≈ FDR 0.05 (Hart 2017 calibration)
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 roughly to FDR 0.05 in MAGeCK; this calibration was confirmed in Hart 2017 G3 by overlap with CEGv2. 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.

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')

Failure Modes

MAGeCK and BAGEL2 disagree by 200+ hits at FDR 0.05

Trigger: Heavy-selection screen (>40% guides change) or cancer-line CN bias. Mechanism: MAGeCK median normalization breaks; BAGEL2 is robust due to reference-set anchoring. Symptom: MAGeCK hit list inflated; BAGEL2 list closer to expected size. Fix: Run MAGeCK with --norm-method control; apply CN correction; trust BAGEL2 for essentiality.

Chronos and MAGeCK disagree at the top 10 in a cancer line

Trigger: Top hits are at amplified loci. Mechanism: Chronos models CN bias; MAGeCK does not. Symptom: ERBB2 in HER2+, MYC in MYC-amplified, etc. are top hits in MAGeCK but not Chronos. Fix: Apply [[copy-number-correction]] before MAGeCK or switch to Chronos.

drugZ and MAGeCK disagree on small-effect drug-modifier screen

Trigger: Effect size is small; MAGeCK rank-based test is less sensitive than drugZ bidirectional Z. Mechanism: drugZ specifically optimized for small effects in drug screens (Li & Hart 2019); MAGeCK RRA loses sensitivity at small effects. Symptom: drugZ has ~2x more hits at the same FDR threshold; MAGeCK misses real chemogenomic interactions. Fix: Use drugZ as primary for chemogenomic; MAGeCK as confirmatory. See [[drugz-chemogenomic]].

JACKS down-weights efficiency, MAGeCK doesn't, disagreement

Trigger: A gene has one or two strong sgRNAs and 2-3 weak ones; MAGeCK averages them, JACKS down-weights the weak. Mechanism: JACKS variational Bayes correctly identifies low-efficacy guides; MAGeCK aggregates without this prior. Symptom: Gene is JACKS hit but not MAGeCK. Fix: Inspect per-sgRNA LFC; if strong guides are consistent, JACKS is correct. Validate gene orthogonally.

Consensus across 3 methods is empty (no hits)

Trigger: Either no real biology, or each method has different failure mode being triggered. Mechanism: Screen quality is low; signal-to-noise across all methods is poor. Symptom: Tier 1 consensus list is empty. Fix: Re-audit QC. Check Cas9 selection, MOI, timepoint, library positioning. Re-run screen if QC fails.

Quantitative Thresholds

Threshold	Value	Source / Rationale
MAGeCK RRA FDR (gene-level)	<0.05	Li 2014; standard publication
MAGeCK RRA LFC	abs(LFC) >1	2-fold; biological
BAGEL2 Bayes Factor	>6 (≈ FDR 0.05); >12 (≈ FDR 0.005)	Kim & Hart 2021; calibrated against CEGv2
drugZ FDR	<0.05 per direction	Li & Hart 2019
JACKS fdr_log10	<-1 (FDR <0.1); <-2 (FDR <0.01)	Allen 2019
Chronos gene-effect probability	>0.7	DepMap convention
Tier 1 consensus (3 methods)	100% agreement	High confidence; minimal validation needed
Tier 2 consensus (2 of 3)	67% agreement	Arrayed validation required
Tier 3 (1 method only)	Hypothesis; flag for follow-up	Multiple screens or arrayed required
Second-best sgRNA rule	Second-best LFC also passes threshold	Reduces single-guide outliers

Common Errors

Error / symptom	Cause	Solution
All genes significant in MAGeCK RRA	Heavy selection breaks median norm	--norm-method control; or use BAGEL2
BAGEL2 returns no hits despite known essentials	Wrong reference gene set	Verify CEGv2/NEGv1 files match library
drugZ output empty	Used Day 0 as control instead of vehicle	Re-run with vehicle as control
Chronos errors out	Missing CN profile for cell line	Use CRISPRcleanR (unsupervised) instead
Methods disagree by orders of magnitude	Quality issue or design mismatch	Re-audit QC; reconcile via tier consensus
Empty tier 1 consensus	No real biology OR QC failure	Re-audit QC
Single-guide-driven hits	Outlier sgRNA	Apply second-best rule; orthogonal validate

References

• Li W et al. 2014. Genome Biol 15:554. MAGeCK alpha-RRA.
• Li W et al. 2015. Genome Biol 16:281. MAGeCK MLE.
• Kim E & Hart T. 2021. Genome Med 13:2. BAGEL2.
• Li G et al. 2019. Genome Med 11:52. drugZ.
• Allen F et al. 2019. Genome Res 29:464. JACKS.
• Dempster J et al. 2021. Genome Biol 22:343. Chronos.
• Meyers R et al. 2017. Nat Genet 49:1779. CERES.
• Hart T & Moffat J. 2016. BMC Bioinformatics 17:164. BAGEL Bayes factor framework.
• Hart T et al. 2017. G3 7:2719. CEGv2/NEGv1 calibration.

Related Skills

• crispr-screens/mageck-analysis - Full MAGeCK RRA + MLE detail
• crispr-screens/bagel-essentiality - Full BAGEL2 detail
• crispr-screens/drugz-chemogenomic - Full drugZ detail for drug screens
• crispr-screens/jacks-analysis - Full JACKS detail and library calibration
• crispr-screens/copy-number-correction - Chronos, CERES, CRISPRcleanR
• crispr-screens/screen-qc - Quality gates that drive method choice
• crispr-screens/library-design - Library type dictates analysis method
• crispr-screens/combinatorial-screens - GI scoring (synthetic lethality)
• crispr-screens/perturb-seq-analysis - SCEPTRE for single-cell screens
• crispr-screens/batch-correction - Multi-batch normalization upstream of hit calling
• pathway-analysis/gsea - Downstream pathway enrichment
• pathway-analysis/go-enrichment - GO enrichment of hit lists