# Bio Applied Mageck Gene Essentiality

> Run MAGeCK count/test on pooled CRISPR sgRNA screens, scoring gene essentiality via RRA, FDR, and log2 fold-change. Use when analyzing CRISPR screen FASTQ/count data, calling essential or drug-resistance genes, or benchmarking vs DepMap.

- Skill: `pavel-kravchenko/bio-applied-mageck-gene-essentiality` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-mageck-gene-essentiality`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-mageck-gene-essentiality/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

| Screen type | Selection pressure | Read-out |
|---|---|---|
| Negative selection (viability) | None (growth over time) | sgRNA depletion → essential genes |
| Positive selection | Drug/toxin | sgRNA enrichment → resistance genes |
| Phenotypic (FACS) | Cell sorting | Enrichment/depletion in gate |

**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%.

```bash
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:

```python
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.

```bash
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
```

```python
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.

```python
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.

```python
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:

```r
# 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