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."

Statistical Models Compared

RRA vs MLE Within MAGeCK

Algorithmic Taxonomy: Why Each Was Built

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:

Reconciliation: When Two Methods Disagree

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

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 fdrsynth and fdrsupp 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 (Colic et al. 2019); MAGeCK RRA loses sensitivity at small effects. Symptom: At matched FDR, drugZ calls small-effect chemogenomic interactions (e.g. DDR genes) that MAGeCK RRA misses, with stronger expected-pathway enrichment. 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

Common Errors

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.
• Colic M 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