# Bio Applied Screen Qc Normalization

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

- Skill: `pavel-kravchenko/bio-applied-screen-qc-normalization` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-screen-qc-normalization`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-screen-qc-normalization/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-screen-qc-normalization

---


# 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.

```python
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 > 0])),
            'median_count': float(np.median(counts)),
            'mean_count': float(np.mean(counts)),
        }
    return metrics


np.random.seed(42)
n_guides = 5000
count_matrix = pd.DataFrame({
    'plasmid': np.random.negative_binomial(5, 0.01, n_guides),
    'replicate_1': np.random.negative_binomial(4, 0.01, n_guides),
    'replicate_2': np.random.negative_binomial(4, 0.01, n_guides),
})

qc = screen_qc_report(count_matrix)
for sample, m in qc.items():
    print(f"{sample}: Gini={m['gini']:.3f}, zero={m['zero_fraction']:.1%}, "
          f"median={m['median_count']:.0f}")
```

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

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

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

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