Bio Applied Mageck Gene Essentiality

CRISPR Screen Analysis with MAGeCK

When to Use

• Turning pooled CRISPR knockout/knockdown screen FASTQ reads into an sgRNA count table (mageck count)
• Calling gene-level essentiality or drug-resistance hits from a count table (mageck test, RRA algorithm)
• Filtering hits by FDR/LFC, or comparing screen results to DepMap CERES/Chronos essentiality scores
• Building volcano plots, rank plots, or Gini-index QC for a CRISPR screen
• Designing screen coverage/replicates before sequencing, or debugging a high unmapped-read fraction

Version Compatibility

• MAGeCK ≥ 0.5.9 (RRA + MLE modules), MAGeCK-VISPR ≥ 0.5.6 for interactive reports
• Python ≥ 3.10 with numpy, pandas, scipy ≥ 1.10, matplotlib
• R ≥ 4.3 with Bioconductor MAGeCKFlute ≥ 2.0 for downstream visualization

Prerequisites

• pip install numpy pandas scipy matplotlib (analysis) and MAGeCK installed via conda install -c bioconda mageck (counting/testing)
• Optional: BiocManager::install("MAGeCKFlute") in R for rank/scatter plots on MAGeCK output
• Familiarity with FASTQ format and basic hypothesis testing/FDR (see bio-read-qc-fastp-workflow, bio-experimental-design-multiple-testing)

Pooled CRISPR Screen Design

A pooled screen introduces thousands of sgRNAs into a cell population — each cell ideally receives one guide via lentiviral transduction at low MOI (≈0.3, single-integration). After selection pressure (drug, growth, immune challenge), deep sequencing of sgRNA abundance reveals which gene knockouts drive the phenotype.

Design rules: genome-wide libraries (Brunello: 4 sgRNAs/gene, ~77K guides; GeCKOv2: 6/gene) need ≥300× coverage (cells per sgRNA) at every timepoint; include 500-1000 non-targeting controls (NTCs) for normalization and false-positive calibration; use ≥2 replicates (3 recommended, count correlation r > 0.95 expected). Multiple sgRNAs per gene let MAGeCK's RRA algorithm distinguish a true on-target effect from one noisy guide.

sgRNA Counting: FASTQ → Count Table

Goal: convert raw screen FASTQ files into a per-sample sgRNA read-count table. Approach: mageck count trims reads to the sgRNA length, exact/near-exact-matches them against the library, and tallies counts per sample; unmapped fraction should stay below ~20%.

mageck count \
    --list-seq library.csv \              # sgRNA library: id,sequence,gene
    --name experiment_name \
    --sample-label "plasmid,d14_rep1,d14_rep2,d14_rep3" \
    --fastq plasmid.fastq.gz d14_rep1.fastq.gz d14_rep2.fastq.gz d14_rep3.fastq.gz \
    --sgrna-len 20 \                      # sgRNA length to extract from read
    --norm-method median                  # normalization: median, total, or none

This produces experiment_name.count.txt (rows = sgRNAs, columns = samples). The snippet below simulates an equivalent table for testing downstream code without running MAGeCK:

import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import nbinom

np.random.seed(42)

# 1000 genes x 4 sgRNAs + 200 non-targeting controls (NTCs), 3 replicates at day 14
N_GENES, SGRNAS_GENE, N_NTC = 1000, 4, 200
N_ESSENTIAL, N_RESISTANCE = 150, 50   # depleted vs enriched gene categories

gene_ids = [f'GENE_{i:04d}' for i in range(N_GENES)]
essential_genes = set(gene_ids[:N_ESSENTIAL])
resistance_genes = set(gene_ids[N_ESSENTIAL:N_ESSENTIAL + N_RESISTANCE])

sgrna_ids, gene_labels = [], []
for gene in gene_ids:
    for j in range(SGRNAS_GENE):
        sgrna_ids.append(f'{gene}_sg{j + 1}')
        gene_labels.append(gene)
for k in range(N_NTC):
    sgrna_ids.append(f'NTC_{k + 1:03d}')
    gene_labels.append('NonTargeting')

n_sgrna = len(sgrna_ids)

# Plasmid (input) counts: overdispersed negative binomial, mean 500
plasmid_mu, plasmid_size = 500, 5
plasmid_counts = nbinom.rvs(n=plasmid_size, p=plasmid_size / (plasmid_size + plasmid_mu), size=n_sgrna)
plasmid_counts = np.maximum(plasmid_counts, 1)

# True log2 fold-change by gene category
lfc_true = np.zeros(n_sgrna)
for i, gene in enumerate(gene_labels):
    if gene in essential_genes:
        lfc_true[i] = np.random.normal(-3.2, 0.5)   # strong depletion
    elif gene in resistance_genes:
        lfc_true[i] = np.random.normal(+2.5, 0.6)   # enrichment
    else:
        lfc_true[i] = np.random.normal(0.0, 0.3)    # non-essential noise


def sim_day14(plasmid, lfc, n_reps=3, dispersion=5):
    """Simulate replicate day-14 counts from plasmid input and a true LFC."""
    reps = []
    for _ in range(n_reps):
        mu = plasmid * 2 ** lfc
        counts = nbinom.rvs(n=dispersion, p=dispersion / (dispersion + mu))
        reps.append(np.maximum(counts, 0))
    return np.array(reps).T   # (n_sgrnas, n_reps)


day14_counts = sim_day14(plasmid_counts, lfc_true)
count_df = pd.DataFrame({
    'sgRNA': sgrna_ids, 'Gene': gene_labels, 'Plasmid': plasmid_counts,
    'Day14_R1': day14_counts[:, 0], 'Day14_R2': day14_counts[:, 1], 'Day14_R3': day14_counts[:, 2],
})
print(count_df.head(8).to_string(index=False))

MAGeCK Test: Robust Rank Aggregation (RRA)

Goal: call gene-level hits from an sgRNA count table by aggregating guide-level fold-changes. Approach: run mageck test, which computes per-sgRNA LFC against a negative-binomial mean-variance model, ranks all guides library-wide, then scores each gene by how extreme its guides sit in that ranking (alpha-RRA), calibrated by permutation.

mageck test \
    -k experiment_name.count.txt \
    -t d14_rep1,d14_rep2,d14_rep3 \        # treatment/late-timepoint samples
    -c plasmid \                          # control/baseline sample
    -n mageck_output \
    --gene-lfc-method median \
    --norm-method median
# writes mageck_output.gene_summary.txt and mageck_output.sgrna_summary.txt

def alpha_rra_score(gene_lfc, null_lfc_pool, n_permutations=1000, seed=0):
    """Educational analogue of MAGeCK's alpha-RRA statistic for one gene.

    Ranks a gene's sgRNA LFCs as percentiles against the library-wide LFC
    pool, then takes rho = min_i Beta_cdf(rank_i; i+1, k-i) over the sorted
    percentiles (an order-statistic extremity score), and calibrates a
    p-value by permuting random guide sets of the same size. Real MAGeCK
    uses a faster analytic approximation of this permutation null.
    """
    rng = np.random.default_rng(seed)
    pool = np.sort(np.asarray(null_lfc_pool))
    ranks = np.sort(np.searchsorted(pool, gene_lfc) / len(pool))
    k = len(ranks)
    rho = min(stats.beta.cdf(ranks[i], i + 1, k - i) for i in range(k))

    null_rhos = np.empty(n_permutations)
    for p in range(n_permutations):
        perm = np.sort(rng.uniform(size=k))
        null_rhos[p] = min(stats.beta.cdf(perm[i], i + 1, k - i) for i in range(k))
    p_value = float(np.mean(null_rhos <= rho))
    return rho, p_value

Hit Calling and Interpretation

Goal: turn mageck_output.gene_summary.txt into a clean list of significant hits. Approach: filter on the RRA output's neg|fdr/pos|fdr and neg|lfc/pos|lfc columns (MAGeCK reports separate negative- and positive-selection statistics per gene), then cross-check known essential-gene categories (ribosome, proteasome, spliceosome) and DepMap Chronos scores as an orthogonal validation.

def call_hits(gene_summary_path, fdr_neg=0.05, fdr_pos=0.05, lfc_neg=-1.0, lfc_pos=1.0):
    """Load a MAGeCK RRA gene_summary.txt and flag significant hits.

    Expected columns: id, num, neg|score, neg|p-value, neg|fdr, neg|lfc,
    pos|score, pos|p-value, pos|fdr, pos|lfc.
    """
    df = pd.read_csv(gene_summary_path, sep='\t')
    df['is_depleted_hit'] = (df['neg|fdr'] < fdr_neg) & (df['neg|lfc'] < lfc_neg)
    df['is_enriched_hit'] = (df['pos|fdr'] < fdr_pos) & (df['pos|lfc'] > lfc_pos)
    return df

Note: DepMap's default essentiality model is Chronos (replaced CERES starting ~21Q3), which jointly corrects for copy-number amplification and cross-library batch effects — treat CERES as legacy when comparing hits.

Visualization

Goal: QC the library and communicate hits. Approach: Gini index checks read-count uniformity of the plasmid pool (want < 0.2); a volcano plot (gene LFC vs -log10 FDR) highlights top hits.

import matplotlib.pyplot as plt


def gini_index(counts):
    """Gini coefficient of sgRNA read-count uniformity (0=even, 1=skewed).
    Plasmid-pool libraries should score < 0.2; higher values flag uneven
    cloning/representation before any biological selection is applied.
    """
    x = np.sort(np.asarray(counts, dtype=float))
    n = len(x)
    cum = np.cumsum(x)
    return (n + 1 - 2 * np.sum(cum) / cum[-1]) / n


def volcano_plot(gene_df, lfc_col='neg|lfc', fdr_col='neg|fdr', label_col='id', top_n=10):
    """Plot gene LFC vs -log10(FDR) from a MAGeCK gene_summary table, labeling top_n hits."""
    fig, ax = plt.subplots(figsize=(6, 5))
    fdr_floor = gene_df[fdr_col].clip(lower=1e-300)
    ax.scatter(gene_df[lfc_col], -np.log10(fdr_floor), s=8, alpha=0.4, color='grey')
    hits = gene_df.nsmallest(top_n, fdr_col)
    ax.scatter(hits[lfc_col], -np.log10(hits[fdr_col].clip(lower=1e-300)), color='red', s=20)
    for _, row in hits.iterrows():
        ax.annotate(row[label_col], (row[lfc_col], -np.log10(max(row[fdr_col], 1e-300))), fontsize=7)
    ax.axhline(-np.log10(0.05), ls='--', color='black', lw=0.8)
    ax.set_xlabel('Gene LFC')
    ax.set_ylabel('-log10(FDR)')
    return fig

For polished downstream figures, MAGeCKFlute (Bioconductor) reads MAGeCK output directly:

# BiocManager::install("MAGeCKFlute")
library(MAGeCKFlute)

gdata <- ReadRRA("mageck_output.gene_summary.txt")
sdata <- ReadsgRRA("mageck_output.sgrna_summary.txt")

# Rank plot highlighting top depleted/enriched genes by RRA score
p1 <- RankView(gdata, top = 5, bottom = 5)
ggsave("mageck_rank_plot.png", p1, width = 6, height = 5)

# 9-square scatter classification (useful for two-condition comparison screens)
p2 <- ScatterView(gdata, x = "Score", y = "Score", label = "id")

Pitfalls

• Low coverage inflates noise: <300x cells/sgRNA lets sampling noise dominate over selection signal, widening confidence intervals and causing both false positives and false negatives
• NTCs drifting: if non-targeting controls show large LFCs, suspect a normalization or coverage problem before trusting any gene hit
• Copy-number artifacts: amplified loci get cut more often (more Cas9 toxicity per sgRNA), making non-essential genes in amplified regions look like false negative-selection hits — correct with CRISPRcleanR before hit calling
• High unmapped fraction (>20%): usually a wrong --sgrna-len/trim offset or a library CSV that doesn't match the sequenced pool — spot-check a few raw reads against the library file
• Multiple testing: MAGeCK's neg|fdr/pos|fdr are already BH-corrected per gene set; don't re-filter on raw p-values across thousands of genes
• CERES vs Chronos: don't assume DepMap scores are still CERES — Chronos has been the default since ~21Q3 and handles copy-number/batch effects differently

See Also

• bio-applied-screen-qc-normalization
• bio-applied-genetic-engineering-in-silico
• bio-applied-statistics-for-bioinformatics
• bio-applied-cancer-transcriptomics