Bio Crispr Screens Bagel Essentiality

Version Compatibility

Reference examples tested with: BAGEL2 2.0 (hart-lab/bagel, build 115), pandas 2.2+, numpy 1.26+, scipy 1.12+, matplotlib 3.8+.

Before using code patterns, verify installed versions match. If versions differ:

• CLI: BAGEL.py fc --help; BAGEL.py bf --help; BAGEL.py pr --help
• Python: BAGEL2 is distributed via git clone (no canonical PyPI release); confirm BAGEL.py version after checkout.

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

BAGEL2 Essentiality Analysis

"Identify essential genes from my CRISPR fitness screen using BAGEL2" -> Compute per-sgRNA fold changes from counts, derive per-gene log-likelihood ratios against reference essential and non-essential gene sets, sum to Bayes Factor, and apply BF threshold calibrated by precision-recall against the reference.

• CLI: BAGEL.py fc to compute fold changes
• CLI: BAGEL.py bf to compute Bayes Factors
• CLI: BAGEL.py pr for precision-recall curves
• Reference sets: CEGv2 (essentials) and NEGv1 (non-essentials); both at https://github.com/hart-lab/bagel

The BAGEL2 Bayesian Framework (under the hood)

Why this matters for postdoc-level use: BAGEL2 uses a Bayes-factor classifier trained on known essential and non-essential genes. The chain:

1. For each sgRNA, compute log-fold-change (LFC) treatment vs control.
2. For each gene, look up per-sgRNA LFCs.
3. For each sgRNA, compute the log-likelihood ratio: log( P(LFC | gene is essential) / P(LFC | gene is non-essential) ). The numerator and denominator are KDEs (kernel density estimates) of LFC distributions from CEGv2 and NEGv1 reference sgRNAs.
4. Sum per-gene log-likelihood ratios across all sgRNAs targeting the gene -> per-gene Bayes Factor.
5. Resampling for the confidence interval (default: 10-fold cross-validation; -b switches to bootstrapping with -NB, default 1000); BF >6 corresponds to ~90% posterior probability (Hart 2017 G3); ~5% FDR by BAGEL convention.

Critical BAGEL2 improvements over BAGEL1:

• Linear extrapolation: BAGEL1 truncated the LLR at the edges of its KDE; BAGEL2 fits a linear regression in the stable region and extrapolates, giving wider dynamic range. This recovers tumor suppressors (highly positive LFC) that BAGEL1 missed.
• Multi-target correction: For sgRNAs targeting multiple genomic loci (off-targets), BAGEL2 discards them and regresses out their BF contribution, but only when -m/--filter-multi-target is given together with --align-info. The original BAGEL counted off-target hits as essentiality signal.
• Tumor suppressor sensitivity: BAGEL2 correctly identifies positive selection (enrichment) genes -- not possible in BAGEL1.

Calibration to CEGv2 / NEGv1

Why these reference sets matter: BAGEL2's discriminative power depends on KDEs of LFCs from known essential vs known non-essential genes. CEGv2 (Hart 2017) is 684 core essential genes shared across cell lines; NEGv1 (Hart 2014) is 927 non-essential genes verified across multiple screens. These act as positive and negative controls within every screen.

Reference set integrity:

• CEGv2: pan-cancer essentials -- common dropouts across most cancer cell lines
• NEGv1: confidently non-essential -- genes without expression or genes with verified neutral status

Critical pitfall: Using a custom essentiality reference (e.g., a single-cell-line CRISPR screen) instead of CEGv2 biases the BAGEL2 model toward that line's specific biology. Always use the standardized references unless there is a specific reason for custom training.

Compute Per-Sample Fold Changes

Goal: Generate per-sgRNA fold-change matrix as input for Bayes-factor calculation.

Approach: Take normalized counts, compute log-fold-change vs a control (Day 0 or plasmid baseline) per sgRNA.

# BAGEL2 installation: distributed via git clone (no canonical PyPI release).
git clone https://github.com/hart-lab/bagel
cd bagel

# Inputs:
# counts.txt: tab-separated with columns: sgRNA, GENE, Sample1, Sample2, ...
# Control column(s): typically Day 0 or plasmid sample(s)
# Treatment column(s): screen endpoint

BAGEL.py fc \
    -i counts.txt \
    -o foldchange \                        # NOTE: -o is a LABEL for fc; writes foldchange.foldchange
    -c Plasmid \                           # control sample (or Day 0)
    --min-reads 30                         # default is 0; 30 is a common convention
# Output: foldchange.foldchange (per-sgRNA LFCs) and foldchange.normed_readcount

Compute Bayes Factors

Goal: Score per-gene essentiality as a Bayes Factor.

Approach: Run BAGEL.py bf with the fold-change matrix and reference gene sets. Resampling defaults to 10-fold cross-validation; add -b -NB N to bootstrap instead.

BAGEL.py bf \
    -i foldchange.foldchange \
    -o bayes_factor.txt \
    -e CEGv2.txt \                         # essentials reference (CEGv2)
    -n NEGv1.txt \                          # non-essentials reference
    -c Sample1,Sample2,Sample3 \            # treatment samples to score
    -b -NB 1000                            # opt into bootstrapping (default is 10-fold cross-validation)
# Output: bayes_factor.txt - per-gene Bayes Factor + CI

Output columns:

Interpretation rule: BF >6 corresponds to ~90% posterior probability of essentiality against CEGv2 (Hart 2017; FDR ≤3% in that calibration, with ~5% a looser BAGEL convention); higher BF = stronger evidence the gene is essential. BAGEL2 also reports negative BFs which can indicate tumor suppressors (positive selection).

Precision-Recall Curve

Goal: Empirically select BF threshold for a given precision/recall tradeoff.

Approach: Run BAGEL.py pr to compute precision and recall at every BF level against CEGv2; pick the BF that gives desired precision.

BAGEL.py pr \
    -i bayes_factor.txt \
    -o precision_recall.txt \
    -e CEGv2.txt \
    -n NEGv1.txt
# Output: precision_recall.txt - precision/recall at each BF threshold

Practical BF ladder (regenerate precision and recall per screen with BAGEL.py pr):

Pick threshold based on application: For exploratory hit calling, BF >0 with low precision is acceptable; for clinical-grade essentiality calls, BF >12 or higher.

Interpret BAGEL2 Results

Goal: Stratify genes into essential, non-essential, and tumor-suppressor categories.

Approach: Apply BF threshold to classify; flag negative BF as candidate tumor suppressors.

import pandas as pd

def interpret_bagel(bf_path, bf_essential=6, bf_tumor_suppressor=-6):
    '''Classify genes from BAGEL2 BF output.'''
    df = pd.read_csv(bf_path, sep='\t')
    df['call'] = 'neutral'
    df.loc[df['BF'] > bf_essential, 'call'] = 'essential'
    df.loc[df['BF'] < bf_tumor_suppressor, 'call'] = 'tumor_suppressor'
    return df.sort_values('BF', ascending=False)

Tumor suppressor identification: Genes with significantly negative BF (e.g., <-6) are enriched in the screen, indicating fitness advantage from their loss. This is biologically distinct from "non-essential" and may indicate tumor-suppressor function. BAGEL1 could not detect this; BAGEL2's linear extrapolation enables it.

Bayesian Reasoning Per Sgrna

Why this matters: BAGEL2 computes per-sgRNA contributions; a gene with 4 sgRNAs each contributing +5 to BF gets +20 total. A gene with 3 sgRNAs contributing +5 and 1 sgRNA contributing -3 (off-target or low-efficacy) gets +12 net.

# Per-sgRNA contributions for diagnosis
# Output table: each sgRNA's LLR contribution to gene-level BF
# Useful for identifying low-efficacy guides

Critical: When per-sgRNA contributions are very heterogeneous (one sgRNA dominates BF), the gene is "guide-of-one"; verify with JACKS efficiency analysis or apply the second-best-sgRNA rule from [[hit-calling]].

Comparing BAGEL2, MAGeCK, drugZ

Reconciliation: BF >6 ≈ 90% posterior probability (Hart 2017 G3); commonly treated as roughly MAGeCK FDR 0.05 by convention. BAGEL2 hits absent from MAGeCK suggest weak signal that BAGEL2's reference anchoring detects but MAGeCK's null-based test misses; verify by inspecting per-sgRNA contributions.