bio-applied-screen-qc-normalization

name	bio-applied-screen-qc-normalization
description	Compute Gini index and replicate LFC correlation on CRISPR sgRNA count matrices; apply DESeq2-style median-ratio normalization before MAGeCK/CRISPRcleanR. Use when QC'ing a CRISPR screen count table or flagging copy-number-biased dropout.
tool_type	python
primary_tool	MAGeCK

CRISPR Screen QC, Normalization, and Copy-Number Bias Detection

When to Use

• You have an sgRNA-by-sample count matrix (from mageck count or a similar tool) and need to check it before calling hits.
• You need to verify library representation evenness (Gini index) and flag samples with too many dropped-out guides.
• You need to normalize raw counts across plasmid/T0 and treatment replicates before differential testing (MAGeCK RRA, CRISPRcleanR, or a custom LFC calculation).
• You need to check replicate concordance (Pearson/Spearman on log-fold-changes) before trusting a screen's hit list.
• You suspect copy-number amplified regions are inflating dropout signal in a cancer cell-line screen.

Version Compatibility

• Python ≥3.10, numpy ≥1.24, pandas ≥2.0, scipy ≥1.10
• MAGeCK ≥0.5.9 (for mageck count/mageck test upstream/downstream of this QC)
• CRISPRcleanR ≥3.0 (R package) if doing copy-number correction on real cancer cell-line data instead of the segment-flagging heuristic below

Prerequisites

• pip install numpy pandas scipy matplotlib
• An sgRNA count matrix: guides as rows, samples (plasmid/T0 + replicates) as columns, e.g. from bio-crispr-screens-mageck-analysis or mageck count
• Familiarity with log-fold-change (LFC) as the basic screen readout: log2((treatment + pseudocount) / (control + pseudocount))

Goal: Decide whether a CRISPR screen's count matrix is good enough to call hits from, and normalize it if so.

Approach: Compute per-sample Gini index and zero-count fraction to catch under-representation; apply median-ratio (DESeq2-style) normalization to remove library-size differences; check replicate LFC correlation; scan for genomic segments with copy-number-biased dropout.

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


def gini_index(counts: np.ndarray) -> float:
    """Gini index of guide count distribution. Perfect evenness = 0, one guide all reads = 1.
    Good screens: Gini < 0.1. Problematic: Gini > 0.2.
    """
    sorted_counts = np.sort(counts)
    n = len(sorted_counts)
    return (2 * np.sum(np.arange(1, n + 1) * sorted_counts) /
            (n * np.sum(sorted_counts)) - (n + 1) / n)


def screen_qc_report(count_df: pd.DataFrame) -> dict:
    """Generate QC metrics for a CRISPR screen count matrix.

    Args:
        count_df: DataFrame with guides as rows, samples as columns

    Returns:
        dict with QC metrics per sample (gini, zero_fraction, median_count, ...)
    """
    metrics = {}
    for col in count_df.columns:
        counts = count_df[col].values
        metrics[col] = {
            'total_reads': int(counts.sum()),
            'guides_detected': int((counts > 0).sum()),
            'zero_count_guides': int((counts == 0).sum()),
            'zero_fraction': float((counts == 0).mean()),
            'gini': float(gini_index(counts[counts > ])),
            : (np.median(counts)),
            : (np.mean(counts)),
        }
     metrics


np.random.seed()
n_guides = 
count_matrix = pd.DataFrame({
    : np.random.negative_binomial(, , n_guides),
    : np.random.negative_binomial(, , n_guides),
    : np.random.negative_binomial(, , n_guides),
})

qc = screen_qc_report(count_matrix)
 sample, m  qc.items():
    (
          )

Goal: Remove library-size differences between samples so raw counts are comparable.

Approach: Median-ratio normalization (the DESeq2 size-factor principle) — divide each guide's counts by a per-guide geometric-mean pseudo-reference, take the median of those ratios per sample as the size factor, then divide raw counts by it.

def median_ratio_normalize(count_df: pd.DataFrame) -> pd.DataFrame:
    """Normalize counts using median-ratio method (DESeq2-style).

    1. Compute geometric mean per guide across samples (pseudo-reference)
    2. Divide each count by the pseudo-reference
    3. Size factor = median of ratios for each sample
    4. Divide raw counts by size factor
    """
    log_means = np.log(count_df + 1).mean(axis=1)
    geo_means = np.exp(log_means)
    ratios = count_df.div(geo_means, axis=0)
    size_factors = ratios.replace([np.inf, -np.inf], np.nan).median(axis=0)
    return count_df.div(size_factors, axis=1)


norm_counts = median_ratio_normalize(count_matrix)
print("Size factors:", (count_matrix.sum() / norm_counts.sum()).round(3).to_dict())

Goal: Confirm two treatment replicates agree before trusting the screen's hit calls.

Approach: Compute per-guide LFC vs. the plasmid/T0 control for each replicate, then correlate the two LFC vectors (Pearson for linear agreement, Spearman for rank agreement). Good screens: Pearson r > 0.8.

def replicate_correlation(count_df: pd.DataFrame, control: str,
                           rep1: str, rep2: str) -> dict:
    """Calculate LFC correlation between two replicates vs a shared control column."""
    pseudo = 0.5
    lfc1 = np.log2((count_df[rep1] + pseudo) / (count_df[control] + pseudo))
    lfc2 = np.log2((count_df[rep2] + pseudo) / (count_df[control] + pseudo))

    pearson_r, pearson_p = stats.pearsonr(lfc1, lfc2)
    spearman_r, spearman_p = stats.spearmanr(lfc1, lfc2)

    return {
        'pearson_r': pearson_r, 'pearson_p': pearson_p,
        'spearman_rho': spearman_r, 'spearman_p': spearman_p,
    }


corr = replicate_correlation(norm_counts, 'plasmid', 'replicate_1', 'replicate_2')
print(f"Replicate concordance: Pearson r={corr['pearson_r']:.3f}, "
      f"Spearman rho={corr['spearman_rho']:.3f}")

Goal: Flag genomic regions where copy-number amplification — not gene essentiality — is driving guide dropout, before calling hits.

Approach: Bin guides into fixed-size genomic segments, compute median LFC per segment, and flag segments whose median deviates more than 2 MADs from the genome-wide median. For real cancer cell-line data, prefer CRISPRcleanR's ccr.GWclean() (circular binary segmentation), which is the validated production tool for this step.

def detect_cn_bias(lfc: pd.Series, guide_locations: pd.DataFrame,
                    segment_size: int = 50) -> pd.DataFrame:
    """Flag genomic segments with copy-number-biased LFC (heuristic; use
    CRISPRcleanR ccr.GWclean() for production cancer cell-line screens)."""
    guide_locations = guide_locations.copy()
    guide_locations['lfc'] = lfc.values
    guide_locations['segment'] = guide_locations['position'] // segment_size

    seg_stats = guide_locations.groupby(['chr', 'segment'])['lfc'].agg(
        ['median', 'count']
    ).reset_index()

    genome_median = lfc.median()
    mad = np.median(np.abs(lfc - genome_median))
    seg_stats['cn_biased'] = np.abs(seg_stats['median'] - genome_median) > 2 * mad
    return seg_stats


guide_locs = pd.DataFrame({
    'chr': np.random.choice(['chr1', 'chr2', 'chr3'], n_guides),
    'position': np.random.randint(0, 10000, n_guides),
})
lfc = np.log2((norm_counts['replicate_1'] + 0.5) / (norm_counts['plasmid'] + 0.5))
cn_result = detect_cn_bias(lfc, guide_locs)
print(f"CN-biased segments: {cn_result['cn_biased'].sum()} / {len(cn_result)}")

Key QC Thresholds

Metric	Good	Acceptable	Poor
Gini index	< 0.1	0.1–0.2	> 0.2
Zero-count guides	< 1%	1–5%	> 5%
Replicate Pearson r	> 0.8	0.6–0.8	< 0.6
Reads per guide (median)	> 300	100–300	< 100
Essential gene depletion	> 5σ	3–5σ	< 3σ

Pitfalls

• Low MOI violations: If MOI > 0.3, multiple guides per cell corrupt phenotype-guide assignment.
• Plasmid library bias: Always compare to a plasmid/T0 control, not just between treatment groups.
• Batch effects: Screen date and cell passage number can dominate biological signal; include batch as a covariate in the model.
• Copy-number confounding: Amplified loci show guide dropout even without essentiality — always check CN bias before calling hits, especially in cancer cell lines.
• Normalizing after filtering zero-count guides inconsistently: apply the same guide filter across all samples before computing size factors, or size factors become incomparable.

See Also

• bio-crispr-screens-mageck-analysis — upstream mageck count/mageck test hit calling
• bio-crispr-screens-hit-calling — downstream statistical hit identification
• bio-crispr-screens-batch-correction — modeling batch effects across screen replicates
• bio-crispr-screens-library-design — designing sgRNA libraries that avoid these QC failure modes