bio-crispr-screens-mageck-analysis

What this skill does

## Version Compatibility

Reference examples tested with: MAGeCK 0.5.9+, MAGeCKFlute 2.0+ (R/Bioconductor), MAGeCK-VISPR 0.5.6+, pandas 2.2+, numpy 1.26+, matplotlib 3.8+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `mageck --version`, `mageck count --help`, `mageck test --help`, `mageck mle --help`
- R: `packageVersion('MAGeCKFlute')`, `?FluteRRA`, `?FluteMLE`

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

## MAGeCK CRISPR Screen Analysis

**"Run MAGeCK on my pooled CRISPR screen"** -> Count sgRNAs from FASTQ, normalize across samples, and rank genes by enrichment or depletion using either the robust rank aggregation (RRA) test for two-condition designs or the maximum-likelihood (MLE) model with explicit design matrix for multi-condition / time-course / drug screens.

- CLI: `mageck count` -> `mageck test` for two-condition RRA
- CLI: `mageck mle` for multi-condition / time-course / multi-cell-line MLE
- R: `MAGeCKFlute::FluteRRA()` / `FluteMLE()` for downstream visualization and pathway analysis
- Python: `mageck-vispr` for interactive QC + result dashboard

## RRA vs MLE Decision Tree

| Experimental design | Recommended | Why |
|---------------------|-------------|-----|
| Two conditions (e.g. drug vs vehicle, treated vs untreated), single cell line, no covariates | `mageck test` (RRA) | RRA is more robust to outlier sgRNAs; faster; default for most published screens |
| Time series (Day 0 -> Day 7 -> Day 14 -> Day 21) | `mageck mle` | RRA cannot model multiple timepoints jointly; MLE estimates per-condition beta scores |
| Multi-cell-line panel (e.g. DepMap-style 5-50 lines) | `mageck mle` with cell-line covariate, or Chronos | MLE handles >2 conditions; Chronos preferred at DepMap scale |
| Paired samples (each replicate matched donor/cell prep) | `mageck mle` with paired design | RRA does not support pairing |
| Combinatorial (treatment x cell line x time) | `mageck mle` with full factorial design | RRA only handles 1 factor |
| Drug screen (vehicle vs drug, multiple doses) | `mageck mle` with dose covariate OR drugZ (preferred for chemogenomic) | drugZ optimized for chemogenomic; see [[drugz-chemogenomic]] |
| Essentiality (Day 0 -> endpoint, single condition) | `mageck test` for simple dropout; `mageck mle` if multi-cell-line | RRA suffices; or BAGEL2 for Bayesian essentiality; see [[bagel-essentiality]] |
| Multi-batch / multi-screen joint analysis | `mageck mle` with batch covariate, JACKS for guide efficacy, or Chronos | See [[batch-correction]] |

**Fails when:**
- Using RRA for time-series and treating each timepoint as a separate test: false-discovery inflation from un-modeled temporal correlation. Use MLE.
- Using MLE without explicit design matrix in a screen with severe batch effects: beta scores include batch variance. Add batch covariate.
- Using MAGeCK for drug screens without vehicle control: comparing drug vs Day 0 conflates drug effect with proliferation. Compare drug vs vehicle, not Day 0. See [[drugz-chemogenomic]].

## The RRA Algorithm (under the hood)

**Why this matters for postdoc-level use:** RRA is not "non-parametric"; it tests whether per-gene sgRNA p-value ranks are more clustered toward the extremes than expected under the uniform null. The chain:

1. `mageck count` normalizes raw counts (median by default) and outputs `*.count_normalized.txt`.
2. `mageck test` computes a per-sgRNA p-value under a NB model fitted to per-gene variance (mean-variance trend stabilized via empirical Bayes shrinkage; same family as edgeR/DESeq2 but simpler dispersion estimation).
3. Per-sgRNA p-values are ranked; ranks are normalized to percentile (rho).
4. For each gene with k sgRNAs, the alpha-RRA score is the minimum over `Pr(beta(i, k-i+1) <= rho_i)` for i in 1..k -- the probability of observing the i-th smallest rank in a uniform sample.
5. Gene-level p-value is computed by permuting sgRNA-to-gene assignment to derive an empirical null over alpha-RRA scores.
6. Multiple-testing correction is Benjamini-Hochberg.

**Critical assumption:** RRA assumes most sgRNAs are non-changing (used to estimate the NB dispersion). If >40% of sgRNAs change, median normalization fails and the dispersion estimate is biased. Symptom: every gene appears significant. Fix: use control-sgRNA normalization (`--norm-method control`) with non-targeting controls as the reference.

## The MLE Model (under the hood)

**Why this matters:** MLE assumes the per-sgRNA log-count is Negative-Binomial with mean determined by a linear combination of condition betas plus an sgRNA-efficiency term (if supplied). The chain:

1. Define a design matrix where rows are samples and columns are conditions; entries are 1 if the sample belongs to the condition, 0 otherwise. A baseline column (all 1s) is required.
2. Per-gene, the model is `log(count_ij) = baseline_i + sum_c (beta_gc * design_jc) + log(sgrna_efficiency_i) + log(size_factor_j)` where i is sgRNA, j is sample, c is condition.
3. The optimizer finds the per-gene beta_gc maximizing the NB likelihood; per-gene Wald-statistic gives a p-value per condition.
4. Each beta represents the log-fold-change in that condition relative to baseline, accounting for sgRNA efficiency.

**Critical assumption:** Beta scores are interpretable only when the design matrix is correctly specified. A common error is omitting batch as a covariate -- the resulting betas absorb batch variance. The fix is to add batch columns; the betas then estimate biology after batch adjustment.

**Convergence failure:** MLE NaN beta scores indicate optimizer divergence -- usually because a gene has too few sgRNAs with non-zero counts in the relevant conditions. The output includes such genes with NaN; do not interpret them as zero effect.

## Count sgRNAs from FASTQ

**Goal:** Quantify sgRNA representation from raw sequencing data.

**Approach:** Run `mageck count` to map FASTQ reads to the sgRNA library reference, producing a normalized count matrix and QC summary. Set `--norm-method median` (default; robust to outliers) or `--norm-method control` (when many sgRNAs change).

```bash
mageck count \
    --list-seq library.csv \                       # sgRNA library: header sgRNA,Gene,Sequence
    --sample-label Plasmid,Day0,Veh_r1,Veh_r2,Drug_r1,Drug_r2 \
    --fastq Plasmid.fq.gz Day0.fq.gz Veh_r1.fq.gz Veh_r2.fq.gz Drug_r1.fq.gz Drug_r2.fq.gz \
    --norm-method median \                         # see normalization decision below
    --output-prefix screen \
    --trim-5 CACCG                                 # 5'-adapter (lentiGuide-Puro convention)

# Outputs:
#   screen.count.txt           raw counts
#   screen.count_normalized.txt normalized counts (median-scaled)
#   screen.countsummary.txt    per-sample QC: Gini, reads, mapping rate, % zero-count
#   screen.log                  per-FASTQ mapping stats
```

**Library file format (tab- or comma-separated; first row is header):**

```
sgRNA,Gene,Sequence
BRCA1_1,BRCA1,ATGGATTTATCTGCTCTTCG
BRCA1_2,BRCA1,CAGCAGATACTTGATGCATC
NTC_0001,NonTargeting_0001,GACGCATCGAATCAATAGCC
```

## Normalization Decision

| `--norm-method` | When to use | Mechanism | Fails when |
|------------------|-------------|-----------|------------|
| `median` (default) | Standard screen, <40% guides change | Scale each sample to a common median of all guides | Heavy selection (>40% guides change) inflates median, biasing scaling |
| `total` | Equal sequencing depth assumed, no outliers | Scale to total reads | Sensitive to PCR jackpots / outlier high-count guides |
| `none` | Already-normalized inputs (DEPRECATED for raw FASTQ counting) | Skips scaling | Almost never appropriate; only for already-normalized inputs |
| `control` | Heavy selection screens; library has ≥500 NTCs | Scale each sample so non-targeting controls have constant median | Requires `--control-sgrna ntcs.txt`