GPTomics/bio-crispr-screens-batch-correction

Version Compatibility

Reference examples tested with: pyComBat 0.3.3+ (epigenelabs/pyComBat), MAGeCK 0.5.9+, R/limma 3.58+, sva 3.50+, RUVSeq 1.36+, pandas 2.2+, numpy 1.26+, scikit-learn 1.4+, scipy 1.12+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show pycombat; from combat.pycombat import pycombat
• R: packageVersion('sva'); ?ComBat; packageVersion('RUVSeq'); ?RUVg

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

Batch Correction for CRISPR Screens

"Correct batch effects in my CRISPR screens" -> Diagnose the batch source, decide whether to remove via empirical-Bayes (ComBat), explicit covariate modeling (MAGeCK MLE / Chronos design matrix), control-guide-anchored normalization, or unwanted-variation decomposition (RUV, SVA), then apply only the correction that preserves biological condition signal.

• Python: pyComBat.pycombat for empirical-Bayes correction
• Python: explicit batch covariates in mageck mle --design-matrix
• R: sva::ComBat, RUVSeq::RUVg, limma::removeBatchEffect
• CLI: chronos-cn natively handles screen-batch covariates

Batch Sources in CRISPR Screens

| Source | Mechanism | Detectable by | |--------|-----------|---------------| | Library lot | Different aliquots or PCR amplifications | Gini shift; plasmid-pool sequencing | | Cell passage cohort | Cells passaged through different periods | PCA Day-0 samples clustering by passage | | Infection day | Lentivirus titer drifts; FBS lot changes | PCA Day-0 samples cluster by day | | Cas9 enzyme lot | Cas9 expression heterogeneity | PR-AUC drift across screens | | Sequencing run | Lane bias, flowcell variant, machine | Per-sample read-count distribution | | FBS / culture lot | Fetal bovine serum lot variations confound proliferation | Day-0 vs endpoint differential not present in vehicle | | Tissue-prep batch | In-vivo: animal cohort, surgical day, organ-prep tech | In-vivo screens (see [[in-vivo-screens]]) |

Critical: Batch effects in CRISPR screens often correlate with biology (e.g., the drug arm was processed in batch 2 because that's when the drug arrived). This confounds correction. Always check for confounding before applying ComBat.

Batch Effect Decision Tree

| Diagnostic finding | Recommended correction | |--------------------|------------------------| | PCA shows samples cluster by condition, not batch | No correction needed; biology dominates | | PCA PC1 separates batches, PC2 separates conditions | Apply ComBat with condition as biological_covariate | | Batch fully confounded with condition (e.g. all drug in batch 2, all vehicle in batch 1) | Correction will destroy biology; instead redesign next screen with cross-batch balance OR re-analyze with batch in MAGeCK MLE design matrix | | Day-0 (pre-perturbation) samples cluster by batch | Strong batch effect; ComBat needed | | Endpoint samples cluster by batch but not Day-0 | Selection-driven artifact (FBS lot etc); correct or include batch as covariate | | Replicates within a batch are tight; across-batch much wider | Classic batch effect; ComBat | | Each replicate scatters randomly across PCs | Sample-level noise; no batch correction will help | | Cancer-line panel with multiple batches | Use Chronos (built-in batch modeling) |

Diagnose: PCA + Variance Decomposition

Goal: Quantify what fraction of variance is batch vs condition before correcting.

Approach: Run PCA on log10(counts+1); fit ANOVA decomposing variance into batch and condition components; report variance explained.

import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from scipy import stats

def batch_diagnostic(counts_df, metadata_df, batch_col='batch', condition_col='condition'):
    '''Variance decomposition: report fraction of PC1/PC2 variance attributable to batch vs condition.'''
    log_counts = np.log10(counts_df + 1).T  # samples as rows
    pca = PCA(n_components=5)
    pcs = pca.fit_transform(log_counts)
    out = pd.DataFrame({
        'PC': range(1, 6),
        'var_explained': pca.explained_variance_ratio_,
    })
    pc_df = pd.DataFrame(pcs, columns=[f'PC{i+1}' for i in range(5)], index=counts_df.columns).join(metadata_df)
    for i in range(5):
        pc = pc_df[f'PC{i+1}']
        f_b, p_b = stats.f_oneway(*[pc[pc_df[batch_col] == b] for b in pc_df[batch_col].unique()])
        f_c, p_c = stats.f_oneway(*[pc[pc_df[condition_col] == c] for c in pc_df[condition_col].unique()])
        out.loc[i, 'batch_F'] = f_b
        out.loc[i, 'batch_p'] = p_b
        out.loc[i, 'cond_F'] = f_c
        out.loc[i, 'cond_p'] = p_c
    return out

Interpretation: If PC1 has batch F-stat > condition F-stat by 10x, batch is dominating and correction is warranted. If condition dominates PC1, no correction needed.

ComBat Empirical-Bayes Correction

Goal: Remove batch-specific location and scale shifts while preserving biological condition signal.

Approach: Log-transform counts, fit ComBat with explicit biological_covariate indicating condition (so the model knows which signal to preserve), back-transform.

import numpy as np
from combat.pycombat import pycombat

def combat_correct(counts_df, batch_vector, condition_vector=None):
    '''ComBat on log-counts with optional biological covariate (condition).
    Preserves condition signal while removing batch shifts.'''
    data = np.log2(counts_df.values + 1)
    if condition_vector is not None:
        mod = pd.get_dummies(condition_vector).values.astype(float)
        corrected = pycombat(data, np.array(batch_vector), mod=mod)
    else:
        corrected = pycombat(data, np.array(batch_vector))
    return pd.DataFrame(np.power(2, corrected) - 1,
                         index=counts_df.index, columns=counts_df.columns).clip(lower=0)

Critical caveat: ComBat assumes batch effects are linear shifts of mean and variance in log space. Non-linear effects (e.g., gene-specific batch sensitivity) remain. Always re-check PCA after correction to confirm batches now overlap.

RUV (Remove Unwanted Variation)

Goal: Identify hidden batch sources via control sgRNAs whose true signal is known.

Approach: Designate non-targeting controls as "negative controls" (assumed unchanged); RUV decomposes their variance into unwanted factors, then subtracts these from all data.

library(RUVSeq)
# counts_df: rows = sgRNAs, columns = samples
ntc_indices <- which(rownames(counts_df) %in% ntc_sgrna_names)
seqset <- newSeqExpressionSet(counts = as.matrix(counts_df))
ruv_corrected <- RUVg(seqset, cIdx = ntc_indices, k = 2)  # k = 2 unwanted factors
# Access corrected data
corrected_counts <- normCounts(ruv_corrected)

When to use: RUV preferred over ComBat when batches are not annotated (e.g., unknown technical confounders). Worse than ComBat when batch is known and well-annotated; ComBat is more direct.

SVA (Surrogate Variable Analysis)

Goal: Estimate unknown latent factors that may confound the screen.

Approach: SVA computes surrogate variables that capture variance not explained by known biological factors; these can then be added to the MAGeCK MLE design matrix as covariates.

library(sva)
# counts_df: rows = sgRNAs, columns = samples
mod <- model.matrix(~ condition, data = metadata)
mod0 <- model.matrix(~ 1, data = metadata)
sv_obj <- sva(as.matrix(counts_df), mod, mod0)
n_sv <- sv_obj$n.sv  # number of surrogate variables
# Add to design matrix for MAGeCK MLE
design_mat <- cbind(mod, sv_obj$sv)

Use case: When the screen has clear biological signal (e.g. essentiality recovery passes) but small effect sizes are hidden by noise; SVA-discovered latent factors as covariates can recover them.

Batch as Explicit Covariate (Preferred for MAGeCK MLE / Chronos)

Goal: Model batch and biology in the same regression instead of pre-correcting.

Approach: Add batch indicator columns to the MLE design matrix. The fitted beta for condition is the effect after accounting for batch; no pre-correction needed.

# Design matrix for a screen with 2 batches and 2 conditions
cat > design.txt <<EOF
Samples         baseline    batch2    treatment
Veh_b1_r1       1           0         0
Veh_b1_r2       1           0         0
Drug_b1_r1      1           0         1
Drug_b1_r2      1           0         1
Veh_b2_r1       1           1         0
Veh_b2_r2       1           1         0
Drug_b2_r1      1           1         1
Drug_b2_r2      1           1         1
EOF

mageck mle \
    --count-table counts.txt \
    --design-matrix design.txt \
    --output-prefix batch_aware_mle

Why this is preferred: ComBat shifts counts before testing; the MLE-with-covariates approach correctly propagates uncertainty from the batch term into the condition beta's standard error. ComBat-then-test pretends the corrected counts are noise-free, biasing FDR.

Control-Sgrna Anchored Normalization

Goal: Use non-targeting controls as the per-sample reference so batch shifts cancel.

Approach: Scale each sample so its NTC sgRNAs have a constant median. Subsequent fold changes are relative to NTCs in each sample, automatically batch-controlling.

def ntc_anchored_normalize(counts_df, ntc_sgrna_names, target_median=1000):
    '''Scale each sample so its NTC median is target_median. Subsequent LFC is NTC-anchored.'''
    is_ntc = counts_df.index.isin(ntc_sgrna_names)
    ntc_medians = counts_df.loc[is_ntc].median(axis=0)
    scale_factors = target_median / ntc_medians.replace(0, np.nan)
    return counts_df * scale_factors, scale_factors

Critical: Requires ≥500 NTCs in the library (see [[library-design]]). With fewer, the NTC median is unstable and amplifies noise rather than removing batch.

When NOT to Correct

| Situation | Why correction hurts | |-----------|----------------------| | Batch is fully confounded with condition | Correction destroys biology along with batch; redesign or accept | | Batch effect is smaller than between-replicate noise | Correction adds noise without removing meaningful variance | | Replicates already correlate >0.95 within and across batches | No batch effect to correct | | Single-screen analysis | No "batch" to correct; only replicate noise | | Per-batch sample size <3 | Cannot estimate batch shift reliably; correction is harmful |

Failure Modes

ComBat eliminates biological signal

Trigger: Batch is correlated with condition (e.g., all drug-arm samples were processed week 2; all vehicle-arm samples week 1). Mechanism: ComBat without a mod covariate treats condition variance as batch variance; corrects it away. Symptom: PR-AUC against CEGv2 drops after ComBat correction. Fix: Always supply mod covariate matrix indicating condition; verify by comparing PR-AUC before and after.

RUV adds noise instead of removing it

Trigger: k (number of unwanted factors) set too high. Mechanism: RUV's least-squares decomposition over-fits; "removed" variance includes biology. Symptom: Hits decrease and replicate Pearson drops after correction. Fix: Choose k via cross-validation; default k=1 or 2 for most screens.

Batch-aware MLE collinear design matrix

Trigger: Adding a batch indicator that is fully collinear with another design column (e.g., all of batch 2 is also Day 21). Mechanism: MLE design matrix is singular; betas not estimable. Symptom: MAGeCK MLE errors out or produces NaN betas. Fix: Drop the collinear column; re-design experiment with cross-batch balance.

ComBat after RUV double-corrects

Trigger: Applying multiple corrections sequentially. Mechanism: Both methods remove variance; sequential application removes biology twice. Symptom: All signal gone; counts look uniformly noisy. Fix: Pick one method based on diagnostic; never combine.

Per-batch sample size too small

Trigger: 2 replicates per batch with 3 batches; ComBat estimates batch shift from 2 samples. Mechanism: Insufficient data to estimate batch parameters; high-variance estimates. Symptom: Correction makes some batches worse than uncorrected. Fix: Need ≥3 (preferably 4-6) samples per batch; below this, use covariate modeling instead.

Quantitative Thresholds

| Threshold | Value | Source / Rationale | |-----------|-------|--------------------| | PC1 batch F vs condition F | F_batch > 10x F_cond -> apply correction | Standard variance-decomposition diagnostic | | ComBat min samples per batch | ≥3, ideally 4-6 | Empirical Bayes prior estimation | | RUV k (unwanted factors) | k=1 default; k=2 if multiple known batch sources | Risso 2014; cross-validate | | NTCs needed for NTC-anchored norm | ≥500 in library | Stable median | | Post-correction PCA check | Batches must overlap in PC1/PC2 plot | Visual sanity check | | Post-correction PR-AUC | Should be same or higher than pre | If lower, correction destroyed biology |

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | PR-AUC drops after ComBat | Batch confounded with condition | Add mod covariate; or redesign | | MAGeCK MLE NaN beta after adding batch column | Collinear design matrix | Drop collinear column | | Replicates still cluster by batch after RUV | k too low | Increase k; cross-validate | | Replicates lose internal cohesion after correction | Over-correction | Reduce k or revert | | NTC-anchored norm worse than median | Too few NTCs | Use median; add NTCs to next library | | Sequencing-run-level batch survives ComBat | Non-linear sequencing effect | Pre-normalize with mageck count --norm-method control first |

References

• Johnson WE et al. 2007. Biostatistics 8:118. Original ComBat algorithm.
• Leek JT et al. 2012. Bioinformatics 28:882. SVA package.
• Risso D et al. 2014. Nat Biotechnol 32:896. RUVSeq.
• Pacini C et al. 2021. Cell Syst 12:1132. CRISPR-screen batch effects analysis.
• Hayer KE et al. 2023. Genome Biol 24:1. Benchmark of normalization methods for CRISPR screens.

Related Skills

• crispr-screens/mageck-analysis - MAGeCK MLE with explicit batch covariates
• crispr-screens/screen-qc - Pre-correction PCA diagnostic
• crispr-screens/copy-number-correction - Chronos handles batch + CN jointly
• crispr-screens/library-design - NTC composition for NTC-anchored normalization
• crispr-screens/jacks-analysis - Joint analysis across batches with shared efficacy
• crispr-screens/hit-calling - Post-correction hit calling
• crispr-screens/in-vivo-screens - In-vivo-specific batch sources (animal cohort, tissue prep)