Batch Effect Correction: The Secrets Behind Those Default Values

· 18 · Tutorials

Batch Effect Correction: The Secrets Behind Those Default Values

You know Harmony's θ default is 2, and scVI's n_latent default is 10—but do you know why? Do you know when adjusting them improves results and when it breaks everything?

This article isn't about "Harmony vs scVI vs Seurat MNN: which is better"—there are already plenty of those. Instead, it covers: what the key parameters mean mathematically when you get a new dataset, why their default values work in most cases, and which boundary conditions cause them to fail.

Building Intuition: What Batch Effect Means Mathematically

The essence of batch effects is that gene expression vectors for the same cell type, measured under different experimental conditions, exhibit systematic shifts in high-dimensional space.

Correction methods generally fall into two categories:

  • Alignment-based (Harmony, Seurat MNN): Pull embeddings of cells from different batches together.
  • Generative model-based (scVI, scANVI, totalVI): Explicitly model "batch as a noise variable" and learn a batch-free latent representation.

The parameter logic differs fundamentally between these two. We'll discuss them separately.

I. Deep Dive into scVI / scANVI Parameters

scVI's core is a Variational Autoencoder (VAE) with the objective function:

$$\mathcal{L} = \underbrace{\mathbb{E}{q(z|x,s)}[\log p(x|z,s)]}{\text{Reconstruction Loss}} - \underbrace{\beta \cdot D_{KL}(q(z|x,s) | p(z))}_{\text{Regularization Loss}}$$

where $x$ is the observed gene count, $s$ is the batch label, and $z$ is the latent space vector. Understanding this formula unlocks all core scVI parameters.

1.1 n_latent: Latent Space Dimensionality

Default: 10

$z \in \mathbb{R}^{n_latent}$ is the "batch-free" cell representation learned by the encoder.

Why is 10 usually sufficient?

In typical scRNA-seq datasets, biological variation in cells is primarily driven by a limited number of cell types/states. Ten dimensions are enough to distinguish 8–12 major cell types in PBMCs and their continuous states. Studies show that for typical datasets of 10k–50k cells, the first 10 PCA components explain > 70% of the variance—the choice of n_latent=10 aligns with this empirical rule.

When should you increase it?

Scenario Suggestion
Dataset > 100k cells, many cell types (e.g., whole-organ atlas) Increase to 20–30
Contains rare cell subpopulations (< 0.5%) requiring fine discrimination Increase to 15–20
Developmental trajectory data (continuous cell states) Increase to 15–20
Only doing batch correction, no downstream subtyping Keep at 10, no need to increase

When does it break?

  • n_latent too large (e.g., 50+): Too many model parameters, prone to overfitting on small datasets (< 5k cells)—training loss looks good, but UMAP still shows separated batch points (the model "hides" batch information in extra dimensions).
  • n_latent too small (e.g., 5): Information bottleneck is too tight, representations of different cell types start overlapping, reducing downstream clustering accuracy.

Quick diagnostic: After training, run model.get_latent_representation(), then perform PCA and look at the explained variance curve for the first n_latent components. If the curve flattens quickly (e.g., near zero after the 6th component), your n_latent can be smaller—the model is wasting capacity.

1.2 n_hidden: Number of Hidden Layer Neurons

Default: 128

The number of neurons in each hidden layer of the encoder/decoder. This controls the capacity of the non-linear mapping.

Why 128 is sufficient: For typical inputs of 2000–5000 highly variable genes, a hidden layer with 128 neurons can capture sufficient non-linear features. Increasing it usually doesn't significantly improve performance on most datasets, only increases training time.

Exception: When the input gene count is very large (e.g., whole-genome transcriptome, 10k+ genes), consider increasing n_hidden=256; otherwise, the first layer's compression step may be too aggressive, losing details.

1.3 dispersion: Dispersion Parameter for Negative Binomial Distribution

Options: "gene" | "gene-batch" | "gene-label" | "gene-cell"
Default: "gene"

scVI models gene counts using a Negative Binomial (NB) or Zero-Inflated Negative Binomial (ZINB) distribution:

$$p(x_{ng} | \mu_{ng}, \theta_g) = \text{NB}(\mu_{ng}, \theta_g)$$

$\theta_g$ (sometimes denoted as dispersion) controls the "spread" of the count distribution for each gene: smaller $\theta_g$ means a more dispersed distribution (heavier tails); as $\theta_g \to \infty$, it approaches a Poisson distribution.

Option Meaning Use Case
"gene" One $\theta_g$ per gene, shared across all cells/batches Most scenarios, default preferred
"gene-batch" One $\theta$ per (gene × batch), allowing different dispersion across batches When sequencing depth varies drastically between batches (e.g., 10x vs Smart-seq2 mix)
"gene-label" One $\theta$ per (gene × cell type) When strong prior knowledge indicates different cell types have vastly different dispersion for the same gene
"gene-cell" One $\theta$ per (gene × cell) Most flexible, but most parameters, requires large datasets (> 50k cells)

Why the default "gene" can break:
When mixing single-end sequencing (10x) and full-length sequencing (Smart-seq2), the count distributions for the same gene differ dramatically—10x produces sparse counts, while Smart-seq2 yields more continuous high counts. Using the same $\theta_g$ to model both batches forces one parameter to fit two completely different distributions. Switching to "gene-batch" can significantly improve ELBO for such datasets.

1.4 gene_likelihood: Probability Distribution for Gene Counts

Options: "zinb" | "nb" | "poisson"
Default: "zinb" (Zero-Inflated Negative Binomial)

ZINB adds a Bernoulli gate on top of NB to model technical zeros:

$$p(x_{ng} = 0) = \pi_{ng} + (1 - \pi_{ng}) \cdot \text{NB}(0 | \mu_{ng}, \theta_g)$$

When to use "nb" instead of "zinb":

  • Newer datasets with higher sequencing depth (> 3000 UMI/cell), fewer zeros
  • Data has been quality-controlled, removing low-quality cells and doublets
  • UMAP shows strange splitting of cell clusters when using "zinb"—this is often the zero-inflation parameter overfitting noise

1.5 β (kl_weight): Weight of KL Divergence

Default: 1.0

This is one of the most overlooked yet impactful parameters.

Revisiting the loss function:

$$\mathcal{L} = \text{Reconstruction Loss} - \beta \cdot D_{KL}(q | p)$$

  • Role of KL term: Pulls the posterior $q(z|x,s)$ towards the standard normal prior $p(z) = \mathcal{N}(0, I)$, forcing the latent space to be smooth, continuous, and interpolatable. This constraint is what "mixes" cells from different batches in the latent space.
  • Role of Reconstruction term: Ensures the decoder can reconstruct the original counts from $z$, preserving biological information.

They are antagonistic:

β Value Effect
β = 0 Pure autoencoder, no regularization, batch effects fully retained in z
β < 1 Light regularization, good reconstruction, weak batch correction, suitable for data with small batch effects
β = 1 (default) Standard VAE, balanced point, suitable for most scenarios
β > 1 (β-VAE) Strong regularization, good batch mixing, but biological details may be compressed and lost
β >> 10 Posterior collapse, all cell z converge, clustering becomes meaningless

Practical tuning suggestions:

  • Start with default kl_weight=1.0, run and check scVI's UMAP
  • If batch separation is still evident: try kl_weight=1.5–2.0
  • If boundaries between cell types are blurred: reduce to kl_weight=0.5–0.8
  • scVI has a kl_weight warm-up mechanism (linearly increases from 0 to the target value over the first few epochs). This warm-up is important; do not disable it.

Note: In scVI's Python API, kl_weight is passed via the train() method using kl_weight or plan_kwargs={"kl_weight": ...}, not during model initialization. The API location varies slightly by version; check the documentation.

1.6 alpha (scANVI-specific): Supervised Signal Strength

scANVI adds cell type label supervision on top of scVI:

$$\mathcal{L}{scANVI} = \mathcal{L}{scVI} + \alpha \cdot \mathcal{L}_{classifier}$$

alpha controls the weight of the classifier loss. Default alpha=1.

  • Increasing alpha: Stronger classifier, more regularized representations for labeled cells, but unlabeled cells rely more on classifier predictions (risk of overfitting).
  • Decreasing alpha: More reliance on the unsupervised part, suitable when the proportion of labeled cells is < 10%.

II. Deep Dive into Harmony2 Parameters

Harmony's core idea is iterative soft clustering + centroid movement:

  1. In PCA space, perform soft K-means clustering on all cells, obtaining membership probabilities for each cell to each cluster.
  2. Calculate the centroid deviation of cells from each batch within each cluster.
  3. Shift each cell's embedding in the opposite direction of its batch centroid (removing batch bias).
  4. Repeat until convergence.

Harmony2 (latest version) introduces soft cosine correction, allowing the correction direction to vary continuously in space, rather than a hard translation per cluster.

2.1 theta (θ): Batch Diversity Penalty

Default: 2 (can be a vector, one value per batch covariate)

$\theta$ controls the penalty strength for batch inhomogeneity when Harmony computes soft cluster centroids:

$$\text{Within-cluster penalty} \propto \theta \cdot \sum_b \frac{n_b}{n} \log \frac{n_b / n}{R_{kb}}$$

where $R_{kb}$ is the membership of cells from batch $b$ to cluster $k$, and $n_b/n$ is the proportion of cells from batch $b$.

Intuition: Larger $\theta$ means Harmony tries harder to make the proportion of cells from each batch within each cluster close to the global proportion, i.e., forces mixing.

θ Value Effect
θ = 0 No penalty, equivalent to no batch correction
θ = 2 (default) Moderate correction, suitable for data with moderate batch effects
θ = 4–6 Strong correction, good batch mixing, but may force different cell types together
θ > 8 Risk of over-correction—biological differences start being removed as batch effects

When to increase θ:

  • Very strong batch differences (different sequencing platforms)
  • Dataset has verified same cell types, but UMAP shows them as two separate clusters

When to decrease θ:

  • Your "batch" is actually a real biological difference (e.g., disease vs normal, different time points)
  • After correction, some cell types that should be separate are merged

θ vector for multiple covariates:

harmony_kwargs = {
    "theta": [2, 1],  # Strong correction for first covariate (sequencing batch), weaker for second (donor)
}

2.2 sigma (σ): Soft Clustering Bandwidth

Default: 0.1

Controls the "fuzziness" of the soft K-means clustering. Smaller $\sigma$ means each cell is more deterministically assigned to one cluster (tending towards hard clustering); larger $\sigma$ means membership distribution is more uniform.

Why 0.1 is usually sufficient: This value ensures that the membership of a cell to the nearest cluster in PCA space is much higher than to others, without becoming a hard assignment.

When it breaks: Data distribution is very uneven (large differences in cell type proportions). With small $\sigma$, rare cell types might be assigned to neighboring large clusters, leading to their "assimilation" after correction. In this case, increase $\sigma$ (e.g., 0.2) to give rare groups more "buffer."

2.3 nclust: Number of Soft Clusters

Default: min(100, N_cells / 30) (approximately)

The number of clusters for Harmony's internal soft K-means. This parameter is often overlooked but has a significant impact.

  • nclust too small: Each "super-cluster" mixes too many different cell types, averaging out real heterogeneity during correction.
  • nclust too large: Exceeds the actual number of biological groups, converges slowly, and centroid estimation becomes unstable when each cluster has too few cells.

Suggestion: At least 5–10 times the number of expected cell types. For a 50k cell dataset, nclust=100–200 is usually reasonable.

2.4 block_size: Proportion of Cells Used per Update

Default: 0.05 (i.e., 5% of cells per iteration)

Harmony does not use all cells to update parameters in each iteration; instead, it samples a block. This is a stochastic gradient approach that:

  • Reduces memory usage
  • Introduces randomness into the convergence process, avoiding local minima

When to adjust:

  • Very small datasets (< 5k cells): block_size=0.05 uses only 250 cells per iteration, too much randomness, unstable convergence. Increase to 0.2–0.5.
  • Very large datasets (> 500k cells): Default 5% is still a lot (25k cells per iteration), can decrease to 0.01–0.02 to speed up each step.

2.5 max_iter_harmony and eps_harmony: Convergence Control

Default: max_iter_harmony=10, eps_harmony=1e-4

Harmony stops when the change in the objective function is less than eps_harmony or the maximum number of iterations is reached.

10 iterations are usually enough, but sometimes not: When batch effects are complex or there are many covariates, observe Harmony's objective function curve (harmony_object.objective_harmony). If it hasn't converged by the 10th iteration (curve still decreasing), increase max_iter_harmony=20.

III. totalVI: Additional Parameters for Multi-modal Data

totalVI processes both RNA and protein (CITE-seq / ASAP-seq) data, adding to scVI:

3.1 Protein Background Prior

CITE-seq protein signals contain non-specific background noise, which totalVI models with a mixture model:

$$p(\text{protein}) = \pi \cdot \mathcal{N}(\mu_{bg}, \sigma_{bg}) + (1-\pi) \cdot \mathcal{N}(\mu_{fg}, \sigma_{fg})$$

  • protein_background_prior_mean (default -10, in log-normal space): Prior mean of the background signal.
  • protein_background_prior_scale (default 0.5): Uncertainty of the background prior.

When these parameters break:
If your protein panel contains many proteins not expressed on target cells (abundant negative controls), the background mean will be systematically high. In this case, protein_background_prior_mean should be increased (e.g., to -5); otherwise, the model will incorrectly classify truly low-expressed proteins as background.

3.2 protein_dispersion

Similar to scVI's dispersion, controls the dispersion modeling for protein counts. Usually keep the default ("protein", one parameter per protein), but when you have multiple measurement batches, "protein-batch" is more appropriate.

IV. General Parameter Tuning Workflow

Regardless of the tool, follow this order:

1. Run with default parameters first
   → Check UMAP: Batch separation? Clear cell type boundaries?
   → Batch separation present → Consider strengthening correction parameters (θ↑ or kl_weight↑)
   → Blurred cell type boundaries → Consider reducing correction strength

2. Tune one parameter, fix others
   → Compare LISI scores (iLISI measures batch mixing, cLISI measures cell type preservation)
   → Look at both scores simultaneously; don't just pursue high iLISI

3. Check biological plausibility
   → Are marker genes highly expressed in the correct clusters?
   → Have rare cell types been "assimilated" and disappeared?

Quick Reference Table

Problem Possible Cause Tuning Direction
Batches still separate in UMAP Insufficient correction strength Harmony: θ↑; scVI: kl_weight↑
Different cell types merged Over-correction Harmony: θ↓; scVI: kl_weight↓
Rare cell types disappear nclust too small / σ too small Harmony: nclust↑, sigma↑
Training doesn't converge Learning rate / warm-up issue Check kl_weight curve, add warm-up
OOM on large datasets block_size too large Harmony: block_size↓; scVI: batch_size↓
ZINB produces strange small clusters Zero-inflation parameter overfitting Switch to gene_likelihood="nb"
scANVI label transfer poor alpha too small alpha↑, or increase proportion of labeled cells

Summary

Choosing batch effect correction parameters is fundamentally about finding a balance between two opposing goals:

  • Batch mixing (integration) vs Biological preservation
  • Model capacity vs Generalization

Default values are reasonable starting points validated on many datasets, but they are not universal. Understanding what each parameter controls mathematically allows you to know where to go when defaults fail—rather than blindly grid searching.

Related tool documentation:

scRNA-seq batch effects scVI Harmony parameter tuning single cell