ComBat-seq

Unlike standard ComBat which operates on continuous log-transformed data, ComBat-seq directly models raw counts using a negative binomial (NB) distribution. This preserves the mean-variance relationship inherent to sequencing data. The model estimates batch parameters using maximum likelihood estimation within a generalized linear model (GLM) framework, avoiding the need for artificial pseudo-counts or log transformations that distort variance structure.

A defining feature of ComBat-seq is its output: adjusted integer counts rather than continuous normalized values. This is critical for downstream analyses that require discrete data:

• Differential expression tools like DESeq2 and edgeR expect count input
• Splicing analysis methods rely on exact count distributions
• Variant calling pipelines require integer allele depths The algorithm achieves this by sampling from the estimated NB distribution after batch effect removal.

ComBat-seq uses a two-step estimation procedure to prevent removal of biological signal:

• Step 1: Fit a full model including both batch and biological covariates to estimate batch parameters
• Step 2: Adjust counts using only the batch component while holding biological effects constant This ensures that genuine condition-specific expression differences are retained while technical artifacts are removed. The model explicitly separates treatment effects from batch effects in the design matrix.

RNA-seq data often contains excess zeros due to genes with low or absent expression. ComBat-seq handles this through the NB distribution's natural accommodation of zero counts without requiring zero-inflation models. The dispersion parameter captures both biological variability and technical noise, preventing overcorrection of stochastic zeros. For extremely sparse datasets, the method remains robust because shrinkage is applied at the gene level rather than globally.

ComBat-seq extends beyond simple two-batch correction to handle complex experimental designs:

• Multiple batches: Simultaneously correct for 3+ processing groups
• Continuous covariates: Adjust for variables like RNA integrity number (RIN) or age
• Interaction terms: Model batch effects that differ across biological conditions The design matrix formulation allows flexible specification of additive and interactive batch effects, making it suitable for large multi-center studies.

ComBat-seq is implemented in the sva (Surrogate Variable Analysis) R/Bioconductor package, ensuring seamless integration with standard genomics workflows:

• Accepts SummarizedExperiment and matrix input objects
• Compatible with DESeq2, edgeR, and limma pipelines
• Returns adjusted counts in the same format as input The function signature ComBat_seq(counts, batch, group = NULL, covar_mod = NULL) mirrors the familiar ComBat interface, minimizing the learning curve for experienced bioinformaticians.

A discrete probability distribution that forms the statistical backbone of ComBat-seq. Unlike the Poisson distribution, the negative binomial explicitly models overdispersion—the phenomenon where the variance of RNA-seq counts exceeds the mean due to biological and technical variability.

• Parameters: Mean (μ) and dispersion (φ)
• Why it matters: ComBat-seq preserves this distributional assumption, unlike log-transform methods
• Tools using it: DESeq2, edgeR, and ComBat-seq all rely on negative binomial GLMs

A statistical technique that borrows information across all genes to stabilize estimates for individual genes with limited data. ComBat-seq applies this principle when estimating batch effect parameters, preventing overfitting for genes with low counts.

• Mechanism: Shrinks gene-specific estimates toward a common prior distribution
• Benefit: Reduces variance of parameter estimates without introducing substantial bias
• Applications: Dispersion estimation in DESeq2, fold change moderation in limma, and batch parameter estimation in ComBat-seq

An unsupervised dimensionality reduction technique essential for visualizing batch effects before and after correction. PCA projects high-dimensional gene expression data onto principal components that capture the greatest variance.

• Pre-correction use: Identify whether samples cluster by batch rather than biology
• Post-correction use: Verify that batch-driven separation has been removed while preserving biological signal
• Interpretation: If PC1 separates batches instead of conditions, correction is needed