bio-crispr-screens-screen-qc

## Version Compatibility Reference examples tested with: MAGeCK 0.5+ (count + VISPR), MAGeCKFlute 2.0+ (R), pandas 2.2+, numpy 1.26+, scikit-learn 1.4+, matplotlib 3.8+, seaborn 0.13+. Before using code patterns, verify installed versions match. If versions differ: - CLI: `mageck --version` then `mageck count --help`; R: `packageVersion('MAGeCKFlute')` - R: `packageVersion('MAGeCKFlute')` then `?BatchRemove` / `?FluteRRA` If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying. ## CRISPR Screen Quality Control **"Audit my CRISPR screen quality before hit calling"** -> Assess library representation, replicate concordance, depth, drift, and biological signal recovery using DepMap-grade metrics, then decide whether the screen is usable, salvageable, or must be repeated. - Python: `pandas` + `scikit-learn` for Gini, AUC, PCA; `MAGeCKFlute` (R) for one-shot QC dashboard - CLI: `mageck count` writes Gini and mapping stats unconditionally to `<prefix>.countsummary.txt` (`GiniIndex`, `Reads`, `Mapped`, `Percentage`); MAGeCK-VISPR for an interactive dashboard ## QC Stage Hierarchy A pooled screen has six distinct bottlenecks where complexity can collapse. Audit each: | Stage | Metric | Acceptable threshold | Failure consequence | |-------|--------|----------------------|----------------------| | Plasmid pool | Gini, skew, % zero-count guides | Gini <0.1, skew <2 (Joung 2017 states <10), zero <0.5% | Missing guides cannot be screened; dropout indistinguishable from non-coverage | | Day-0 infection | Library coverage, MOI verification | ≥99% guide detection at 500x cells/sgRNA; MOI 0.3 | Founder effects; polyclonality with high MOI | | Selection (puro/blast) | % cells surviving, time-course Gini | 30-40% survival at 5-7 days; Gini drift <0.05 | Selection artifact; fast-growers enriched | | Endpoint | Replicate correlation, depth | Pearson >=0.8 on log-counts (MAGeCK-VISPR floor), Spearman >0.7-0.8, >500 reads/sgRNA (Joung 2017 screening) | Noise dominates; FDR inflates | | Biological signal | CEGv2 PR-AUC, NEGv1 false-positive rate | PR-AUC >0.7 at FDR 5% (community "passing" convention); CEGv2 enrichment in top 1k | Screen lacks essentiality signal; hits not credible | | Copy-number artifact | Amplified-region enrichment, sgRNA-cut-count correlation | No correlation between sgRNA off-target count and depletion | False-positive essentiality at amplicons; ERBB2 in HER2+ etc. | Each metric below quantifies one of these stages. ## Library Representation Metrics **Goal:** Detect dropout, oversaturation, and library bottlenecks at each sequencing stage. **Approach:** Compute per-sample zero-count fraction, low-count fraction (<30 reads, the CRISPRcleanR `ccr.NormfoldChanges` default), and percentile-based skew, then track how these change between plasmid -> Day-0 -> endpoint to localize the bottleneck. ```python import pandas as pd import numpy as np def library_representation(counts_df): '''Per-sample library coverage diagnostics. counts_df: rows = sgRNAs, columns = samples (numeric counts).''' out = pd.DataFrame(index=counts_df.columns) out['n_sgrnas_detected'] = (counts_df > 0).sum() out['pct_zero'] = (counts_df == 0).sum() / len(counts_df) * 100 out['pct_lowcount'] = (counts_df < 30).sum() / len(counts_df) * 100 out['median_count'] = counts_df.median() out['p10_count'] = counts_df.quantile(0.10) out['p90_count'] = counts_df.quantile(0.90) out['skew_ratio'] = out['p90_count'] / out['p10_count'].replace(0, np.nan) return out def stage_specific_thresholds(): '''Stage conventions: Joung 2017 (zero-count, skew) + MAGeCK-VISPR (Gini).''' return { 'plasmid': {'pct_zero_max': 0.5, 'skew_max': 2.0, 'gini_max': 0.10}, # skew 2.0 is a stricter modern convention; Joung 2017 states <10 'day_0': {'pct_zero_max': 1.0, 'skew_max': 2.5, 'gini_max': 0.12}, 'endpoint': {'pct_zero_max': 5.0, 'skew_max': 10.0, 'gini_max': 0.30}, } ``` **Interpretation:** Plasmid pool failing Gini <0.1 indicates synthesis or amplification bias; the screen is unfit for use. Endpoint Gini drifting above 0.30 indicates either heavy biological selection (acceptable for strong-phenotype drug screens) or a bottleneck (must be diagnosed). The Day-0 vs plasmid delta isolates whether the issue arose during infection (cloning is unlikely to lose specific guides between extraction and infection -- the change happens in cells). ## Gini Coefficient **Goal:** Quantify how unevenly reads are distributed across sgRNAs in a single sample. **Approach:** Sort non-zero counts ascending, compute Gini via the cumulative-fraction formula. Compare against stage-specific thresholds. ```python def gini(x): '''Gini coefficient: 0 = perfect equality, 1 = maximal inequality. Uses non-zero counts only; zero-count sgRNAs handled separately by % zero.''' x = np.sort(x[x > 0].astype(float)) if x.size == 0: return np.nan n = x.size cumx = np.cumsum(x) return (n + 1 - 2 * np.sum(cumx) / cumx[-1]) / n ``` **Stage-specific thresholds:** only the plasmid Gini <=0.1 is a published cutoff (MAGeCK-VISPR); the remaining grades are operational convention. | Stage | Excellent | Acceptable | Concerning | Failure | |-------|-----------|------------|------------|---------| | Plasmid pool | <0.10 | <0.15 | 0.15-0.20 | >0.20 | | Day 0 (post-infection) | <0.12 | <0.18 | 0.18-0.25 | >0.25 | | Endpoint (post-selection) | <0.30 | <0.40 | 0.40-0.55 | >0.55 | A Gini that climbs from 0.10 (plasmid) to 0.45 (endpoint) is expected when the screen exerts strong selection (drug, lethal-condition). A Gini that climbs to 0.45 without any biological selection (e.g., a control-vs-control timepoint comparison) indicates technical drift. ## Replicate Concordance **Goal:** Verify that biological/technical replicates agree before testing for between-condition differences. **Approach:** Compute pairwise Pearson on log10(counts+1) (MAGeCK-VISPR convention) and Spearman ρ on raw rank, between every replicate pair within a condition. Flag any pair below the MAGeCK-VISPR floor of 0.8 Pearson on log-scale. ```python def replicate_concordance(counts_df, condition_map): '''condition_map: {condition_name: [sample_col1, sample_col2, ...]}.''' log_counts = np.log10(counts_df + 1) rows = [] for cond, samples in condition_map.items(): if len(samples) < 2: continue for i in range(len(samples)): for j in range(i+1, len(samples)): r_pearson = log_counts[[samples[i], samples[j]]].corr().iloc[0, 1] r_spearman = counts_df[[samples[i], samples[j]]].corr(method='spearman').iloc[0, 1] rows.append({'condition': cond, 'rep1': samples[i], 'rep2': samples[j], 'pearson_log': r_pearson, 'spearman': r_spearman}) return pd.DataFrame(rows) ``` **Thresholds:** 0.8 Pearson is the MAGeCK-VISPR floor; the stricter grades are operational convention. | Metric | Excellent | Acceptable | Failure | |--------|-----------|------------|---------| | Pearson on log10(counts+1) | >0.95 | >0.85 | <0.80 | | Spearman on raw ranks | >0.85 | >0.70 | <0.60 | **When Pearson is high but Spearman is low**, a few outlier sgRNAs are driving correlation (one extreme guide dominates). Inspect the scatterplot; typically caused by PCR jackpotting at a single guide. Hit calling should use a method that ranks (RRA, drugZ) rather than one that fits per-sgRNA fold change directly. ## Essentialome Recovery (CEGv2 PR-AUC) **Goal:** Verify the screen has detectable biological essentiality signal by checking whether known essentials (Hart 2017 CEGv2) drop out faster than known non-essentials (NEGv1). **Approach:** Compute precision-recall AUC where positives are CEGv2 genes and negatives are NEGv1; the screen "passes" if PR-AUC >0.7 (community convention; the CEGv2/NEGv1 sets come from Hart 2017 / Hart 2014). ```python from sklearn.metrics import precision_recall_curve, auc, roc_auc_score def essentialome_recovery(gene_lfc_df, cegv2_set, negv1_set): '''gene_lfc_df: must have ["gene", "lfc"] columns (gene-level mean LFC, negative = depleted). cegv2_set, negv1_set: sets of gene symbols from Hart 2017.''' labeled = gene_lfc_df[gene_lfc_df['gene'].isin(cegv2_set | negv1_set)].copy() labeled['is_essential'] = labeled['gene'].isin(cegv2_set).astype(int) y_score = -labeled['lfc'] # negative LFC = depleted = more essential -> higher score precision, recall, _ = precision_recall_curve(labeled['is_essential'], y_score) return { 'pr_auc': auc(recall, precision), 'roc_auc': roc_auc_score(labeled['is_essential'], y_score), 'n_essential_detected': labeled['is_essential'].sum(), 'n_nonessential_detected': (1 - labeled['is_essential']).sum(), } ``` **Source / threshold:** Hart 2017 *G3* 7:2719 defines CEGv2 (~684 core essentials); NEGv1 (~927 non-essentials) comes from Hart 2014 *Mol Syst Biol* 10:733. Both lists at https://github.com/hart-lab/bagel/blob/master/CEGv2.txt and NEGv1.txt. DepMap convention: PR-AUC >0.7 at FDR 5% is the "passing" threshold; <0.5 means the screen has no essentiality signal and is not interpretable. **When PR-AUC is low despite good Gini and Pearson**: cause is usually one of (a) Cas9 was not selected for before screen start (lots of Cas9-negative cells in the pool diluting signal), (b) puromycin selection truncated too aggressively (over-bottleneck), (c) the timepoint is too early (need 14-21 days for KO + decay + selection to manifest). Each has a different remediation. ## Copy-Number Amplicon Bias Diagnostic **Goal:** Detect the Aguirre 2016 / Munoz 2016 copy-number artifact where sgRNAs targeting amplified loci appear "essential" purely from DNA-damage burden. **Approach:** Bin genes by copy number (if known from matched WGS/SNP-array) and check whether mean LFC correlates with CN. Alternatively, count off-target cut sites per sgRNA and check correlation with depletion -- amplified loci share many identical cut sites. ```python def cn_bias_diagnostic(gene_lfc_df, cn_df): '''cn_df: per-gene copy number (from WGS/SNP-array/matched ASCAT). Tests whether amplified genes show systematically lower LFC.''' merged = gene_lfc_df.merge(cn_df, on='gene') bins = pd.qcut(merged['copy_number'], q=5, duplicates='drop') bin_lfc = merged.groupby(bins, observed=True)['lfc'].agg(['mean', 'median', 'std', 'count']) from scipy.stats import spearmanr rho, p = spearmanr(merged['copy_number'], merged['lfc']) return {'cn_vs_lfc_rho': rho, 'cn_vs_lfc_p': p, 'amplified_mean_lfc': merged[merged['copy_number'] > 4]['lfc'].mean(), 'diploid_mean_lfc': merged[(merged['copy_number'] >= 1.5) & (merged['copy_number'] <= 2.5)]['lfc'].mean(), 'per_bin': bin_lfc} ``` **Interpretation:** A Spearman ρ < -0.1 between copy number and LFC indicates copy-number artifact. The diagnostic threshold is conservative -- Aguirre 2016 showed the effect scales with copy number and with the number of cut sites per sgRNA. Remediation: use CRISPRcleanR, CERES, or Chronos (see [[copy-number-correction]]) before hit calling. ## Sequencing Depth Audit **Goal:** Verify that sequencing depth is sufficient to resolve fold changes at the smallest interesting effect size. **Approach:** Compute reads/sgRNA per sample and the coefficient of variation (CV) of total reads across samples. Compare against Joung 2017's >100 reads/sgRNA for plasmid QC and >500 for screening, or MAGeCK-VISPR's 300x. ```python def depth_audit(counts_df): '''Verify depth: Joung 2017 recommends >100 reads/sgRNA for plasmid QC and >500 for screening; MAGeCK-VISPR uses 300x.''' total = counts_df.sum() n_sgrnas = len(counts_df) depth = total / n_sgrnas cv = total.std() / total.mean() return pd.DataFrame({'total_reads': total, 'reads_per_sgrna': depth, 'depth_grade': np.where(depth < 100, 'FAIL', np.where(depth < 300, 'CAUTION', np.where(depth < 500, 'OK', 'EXCELLENT')))}).assign(across_sample_cv=cv) # 100 = Joung 2017 plasmid-QC floor; 300 = MAGeCK-VISPR; 500 = Joung 2017 screening ``` **CV interpretation:** CV >0.5 across samples in total reads indicates demultiplexing imbalance or library-pooling error; even if individual samples pass depth thresholds, the relative count is then biased. ## MOI Verification **Goal:** Confirm that infection occurred at MOI 0.3-0.5 so that ≤1 sgRNA/cell predominates. **Approach:** From titration plate (control wells with serial-diluted virus), compute infection efficiency, then verify by qPCR of integrated proviral copy number in the screen pool. | MOI | P(≥1 sgRNA/cell) | P(≥2 sgRNAs/cell) | Cells with 2+ guides as fraction of infected |