bio-crispr-screens-mageck-analysis

## 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: sgRNA id, sequence, gene (in that order) --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 AUTO # length(s) to trim from 5' end; AUTO detects it (int or comma-list also accepted) # 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,Sequence,Gene BRCA1_1,ATGGATTTATCTGCTCTTCG,BRCA1 BRCA1_2,CAGCAGATACTTGATGCATC,BRCA1 NTC_0001,GACGCATCGAATCAATAGCC,NonTargeting_0001 ``` ## 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` listing NTC sgRNAs | **Diagnostic:** Run with `median` first. If essentialome PR-AUC against CEGv2 is high and Gini is in range, accept. If PR-AUC <0.5 despite reasonable Gini, retry with `control` -- this isolates whether median normalization was masking the signal. ## MAGeCK Test (RRA for Two-Condition) **Goal:** Identify genes significantly enriched or depleted between treatment and control. **Approach:** Compute per-sgRNA NB-model p-values, rank, apply alpha-RRA to get per-gene scores, permute for gene-level FDR. Outputs separate columns for positive and negative selection. ```bash mageck test \ --count-table screen.count.txt \ --treatment-id Drug_r1,Drug_r2 \ --control-id Veh_r1,Veh_r2 \ --norm-method median \ --output-prefix drug_vs_veh \ --gene-lfc-method median \ # alternative: mean (less robust) --sort-criteria pos # rank for positive selection (choices: neg, pos; default neg) # Outputs: drug_vs_veh.gene_summary.txt (gene-level), drug_vs_veh.sgrna_summary.txt ``` **Gene-summary columns:** | Column | Meaning | |--------|---------| | `id` | Gene symbol | | `num` | Number of sgRNAs in library for this gene | | `neg|score`, `pos|score` | alpha-RRA score, negative- or positive-selection direction | | `neg|p-value`, `pos|p-value` | Permutation p-value | | `neg|fdr`, `pos|fdr` | BH-corrected FDR | | `neg|rank`, `pos|rank` | Rank by score (1 = top hit in that direction) | | `neg|lfc`, `pos|lfc` | Median log-fold-change of sgRNAs in that direction | **Interpretation rule:** A gene is a strong essential if `neg|fdr < 0.05` AND `neg|lfc < -1`. A gene is a strong resistance hit if `pos|fdr < 0.05` AND `pos|lfc > 1`. Genes with `fdr < 0.05` but `|lfc| < 0.5` are statistically significant but biologically weak -- worth flagging for orthogonal validation. ## MAGeCK MLE (Multi-Condition) **Goal:** Estimate gene effects across complex experimental designs with multiple conditions, time points, batches, or covariates. **Approach:** Specify a design matrix mapping samples to conditions, then run `mageck mle` which fits a NB GLM with sgRNA-efficiency term and outputs per-gene per-condition beta scores. ```bash # Design matrix: design.txt (tab-separated) # Samples must match sample-label in mageck count output # baseline column must be present and all 1 cat > design.txt <<EOF Samples baseline day7 day14 day21 Day0 1 0 0 0 Day7_r1 1 1 0 0 Day7_r2 1 1 0 0 Day14_r1 1 0 1 0 Day14_r2 1 0 1 0 Day21_r1 1 0 0 1 Day21_r2 1 0 0 1 EOF mageck mle \ --count-table screen.count.txt \ --design-matrix design.txt \ --output-prefix timecourse_mle \ --norm-method median \ --sgrna-efficiency efficiency.txt \ # OPTIONAL: from JACKS or library design --sgrna-eff-name-column 0 \ # 0-based; 0 is the default --sgrna-eff-score-column 1 \ # 0-based; 1 is the default --max-sgrnapergene-permutation 40 # skip genes with more sgRNAs than this (default 40) # Outputs: timecourse_mle.gene_summary.txt with beta scores per condition + Wald p-values ``` **Gene-summary columns (MLE):** | Column | Meaning | |--------|---------| | `Gene` | Gene symbol | | `sgRNA` | Number of sgRNAs | | `<condition>|beta` | Effect-size estimate (log-fold-change relative to baseline) | | `<condition>|z` | Wald z-statistic | | `<condition>|p-value` | Two-sided p-value | | `<condition>|fdr` | BH-corrected FDR | | `<condition>|wald-fdr` | Wald-statistic-based FDR (alternative) | **Interpretation rule:** Beta scores are log2-fold-changes; a beta of -1 in condition day21 means sgRNAs are 2-fold depleted in day-21 relative to baseline. NaN betas indicate convergence failure (typically a gene with too few non-zero counts in that condition); exclude from interpretation. ## Sample MAGeCK Test for Drug Screen with sgRNA Efficiency ```bash mageck test \ --count-table screen.count.txt \ --treatment-id Drug_r1,Drug_r2,Drug_r3 \ --control-id Veh_r1,Veh_r2,Veh_r3 \ --norm-method control \ # NTCs as normalization reference --control-sgrna ntcs.txt \ # one NTC sgRNA per line --gene-lfc-method median \ --variance-estimation-samples Day0_r1,Day0_r2 \ # estimate variance from these samples --output-prefix drug_screen_normalized # --sgrna-efficiency / --sgrna-eff-name-column / --sgrna-eff-score-column are `mageck mle` options, # not `mageck test` options; pass them to mle if guide-efficacy weighting is needed. ``` **Reading order:** Run `mageck count` -> screen-qc skill for QC -> `mageck test` or `mageck mle` -> `MAGeCKFlute` for visualization -> downstream pathway analysis. ## Time-Course Analysis **Goal:** Identify genes with consistent depletion or enrichment across multiple time points. **Approach:** Run `mageck mle` with a design matrix where each timepoint is a separate column, then test for monotonic trends in per-condition beta scores. ```python import pandas as pd import numpy as np def time_course_consistency(mle_results, conditions=['day7', 'day14', 'day21']): '''Identify genes with monotonic beta trends across timepoints. Returns genes where all betas same sign and trend is monotone.''' beta_cols = [f'{c}|beta' for c in conditions] fdr_cols = [f'{c}|fdr' for c in conditions] df = mle_results[['Gene'] + beta_cols + fdr_cols].copy() df['all_negative'] = (df[beta_cols] < 0).all(axis=1) df['all_positive'] = (df[beta_cols] > 0).all(axis=1) df['monotone'] = df[beta_cols].apply(lambda x: (np.diff(x) <= 0).all() or (np.diff(x) >= 0).all(), axis=1) df['any_sig'] = (df[fdr_cols] < 0.05).any(axis=1) return df[df['monotone'] & df['any_sig']].sort_values(beta_cols[-1]) ``` ## Visualizing Results **Goal:** Generate publication-grade volcano plot and rank plot of MAGeCK output. **Approach:** Load gene_summary.txt, plot `-log10(fdr)` vs LFC, color by significance, annotate top hits. ```python import matplotlib.pyplot as plt import numpy as np def volcano(gene_summary_path, direction='neg', fdr_threshold=0.05, lfc_threshold=1.0): '''direction: "neg" for dropout, "pos" for enrichment.''' df = pd.read_csv(gene_summary_path, sep='\t') lfc_col = f'{direction}|lfc' fdr_col = f'{direction}|fdr' fig, ax = plt.subplots(figsize=(9, 7)) sig = (df[fdr_col] < fdr_threshold) & (np.abs(df[lfc_col]) > lfc_threshold) ax.scatter(df.loc[~sig, lfc_col], -np.log10(df.loc[~sig, fdr_col].clip(lower=1e-10)), c='lightgray', alpha=0.4, s=10)