Batch Handling in RNA-Seq Count Models

Why Count Models Keep Normalization and Batch Separate

RNA-seq data are recorded as nonnegative integer counts, but counts from different samples may not be directly comparable. A sample with greater sequencing depth can have larger counts even when its underlying expression profile is similar to that of another sample. Earlier material on count scaling motivates representing this sample-level exposure difference with a normalization factor. A known technical batch, by contrast, is a sample characteristic that may be associated with expression after accounting for the normalization factor.

These are distinct roles in a count model. A normalization factor rescales the expected count through an offset. A batch covariate estimates a batch-associated difference in the expected count while holding the other model terms fixed. Neither operation means that the observed count should be replaced by a continuously adjusted value before fitting the count likelihood.

Definition: Let $Y_{ij}$ be the observed RNA-seq count for gene $i$ in sample $j$. Let $s_j>0$ be the sample-specific normalization factor, treated as known for the gene-wise model, and let $x_j$ and $b_j$ denote biological and batch covariates. A negative-binomial generalized linear model can be written $$ Y_{ij}\sim \operatorname{NB}(\mu_{ij},\alpha_i),\qquad \operatorname{Var}(Y_{ij})=\mu_{ij}+\alpha_i\mu_{ij}^2, $$ $$ \log(\mu_{ij})=\log(s_j)+\beta_{i0}+x_j^\mathsf{T}\beta_i+b_j^\mathsf{T}\gamma_i. $$ Here $\mu_{ij}$ is the expected count, $\alpha_i\geq 0$ is the gene-specific dispersion, and $\log(s_j)$ is an offset: its coefficient is fixed at 1 rather than estimated. The vectors $x_j$ and $b_j$ encode biological and batch covariates; their coefficients are estimated from the data.

The negative-binomial variance allows counts to vary more than a Poisson model, whose variance equals its mean. The model is applied gene by gene; it does not require counts from different genes to be independent for the purpose of specifying each gene's likelihood. It does assume an appropriate mean structure, a suitable dispersion for the gene, and that samples are independent conditional on the modeled covariates, unless the design calls for an additional dependence structure. The normalization factors must be positive and must represent sample-scale differences rather than the batch effect one hopes to estimate.

For a categorical batch with two levels, choose one level as the reference. The indicator $B_j$ is 0 for the reference batch and 1 for the other batch, so its coefficient $\gamma_i$ estimates the log mean-count ratio between batches, conditional on the other terms. With several levels, use reference-level indicator columns or an equivalent full-rank coding. The intercept and reference level must not be represented by redundant columns.

What the Offset Does

The offset enters on the log-mean scale because a sample with normalization factor $s_j$ has expected count proportional to that factor, all else equal. For a single gene with no other covariates, the model says $\mu_{ij}=s_j\exp(\beta_{i0})$. Thus dividing the expected count by $s_j$ gives the model's relative expression scale, $\exp(\beta_{i0})$. The fitted count itself remains a count response, and its likelihood is evaluated at the observed integer $Y_{ij}$.

Worked Example: Applying a Normalization Offset

Suppose a hypothetical gene has observed counts 80 and 150 in two samples. Their normalization factors are $s_1=0.8$ and $s_2=1.5$. Compute the simple scaled counts and compare them with the offset interpretation.

For descriptive scaling, divide each observed count by its sample factor: $$ 80/0.8=100,\qquad 150/1.5=100. $$ The two scaled values are both 100 on the relative scale. This does not imply the original observations were 100 counts; the observed count data remain 80 and 150.

In a count model with no biological or batch covariates, if $\exp(\beta_{i0})=100$, the expected counts are $0.8(100)=80$ and $1.5(100)=150$. Equivalently, the log means include $\log(0.8)+\log(100)$ and $\log(1.5)+\log(100)$. The offset accounts for their different sample scales without altering either observed count.

Normalization factors are often estimated from the count data, then supplied to the model as fixed offsets. Their construction and adequacy matter: factors that fail to capture relevant sample-scale differences can leave systematic variation in the model. A factor is not a batch label, and a batch label is not a substitute for a factor. If every sample in one batch has a different sequencing scale, the model may need both terms: the factor handles the scale specified by normalization, while the batch coefficient estimates any remaining batch-associated shift.

Representing Known Batch in the Design

A known batch is represented by columns in the model's design matrix, alongside the biological terms needed to answer the analysis question. For a two-level condition and two-level batch, one possible model for gene $i$ is $$ \log(\mu_{ij})=\log(s_j)+\beta_{i0}+\beta_{iC}C_j+\gamma_{iB}B_j, $$ where $C_j=1$ for treatment and 0 for control, and $B_j=1$ for batch B and 0 for batch A. The condition effect is $\beta_{iC}$; the batch effect is $\gamma_{iB}$. These effects are additive on the log-mean scale, meaning multiplicative on the expected-count scale.

Worked Example: Interpreting a Fitted Batch-Adjusted Model

Consider a hypothetical gene with fitted coefficients $\beta_{i0}=\log(20)$, $\beta_{iC}=\log(1.5)$, and $\gamma_{iB}=\log(2)$. For a sample with $s_j=1.2$, calculate the fitted expected counts for control in batch A, treatment in batch A, control in batch B, and treatment in batch B.

For control in batch A, $C_j=B_j=0$: $$ \hat\mu_{ij}=1.2\exp(\log(20))=1.2(20)=24. $$ For treatment in batch A, $C_j=1$ and $B_j=0$: $$ \hat\mu_{ij}=1.2\exp\{\log(20)+\log(1.5)\} =1.2(20)(1.5)=36. $$ For control in batch B, $C_j=0$ and $B_j=1$: $$ \hat\mu_{ij}=1.2\exp\{\log(20)+\log(2)\} =1.2(20)(2)=48. $$ For treatment in batch B, both indicators equal 1: $$ \hat\mu_{ij}=1.2\exp\{\log(20)+\log(1.5)+\log(2)\} =1.2(20)(1.5)(2)=72. $$

Conditional on batch and the normalization factor, the fitted treatment-to-control mean-count ratio is $\exp(\beta_{iC})=1.5$, a 50% higher expected count for treatment in this model. Conditional on condition and normalization, the batch B to batch A ratio is $\exp(\gamma_{iB})=2$. These are model-based associations, not evidence that either condition or batch causally changed expression.

This coefficient interpretation relies on the stated model: there is no condition-by-batch interaction, and the batch effect is the same log ratio across the modeled conditions. If scientific aims or the design require different condition effects by batch, an interaction term may be considered, but it needs data that support estimating it. Including every possible term by default is not a remedy for a poorly supported design.

Worked Example: Checking Whether Batch and Condition Can Be Distinguished

Suppose a hypothetical experiment has two samples: one control sample in batch A and one treatment sample in batch B. With an intercept, a treatment indicator $C_j$, and a batch-B indicator $B_j$, the design rows are $$ \begin{array}{c|ccc} &1&C_j&B_j\\ \text{control, batch A}&1&0&0\\ \text{treatment, batch B}&1&1&1 \end{array} $$

The condition and batch columns are identical: $C_j=B_j$ for both samples. Therefore, the design matrix does not have three independent columns. In the linear predictor, $\beta_{iC}C_j+\gamma_{iB}B_j=(\beta_{iC}+\gamma_{iB})C_j$, so the data can inform the sum of these coefficients but cannot determine them separately. Changing normalization factors does not resolve this problem, because an offset does not create a sample with condition and batch separated.

In a better-balanced hypothetical design with samples in all four condition-by-batch combinations, the rows include $(1,0,0)$, $(1,1,0)$, $(1,0,1)$, and $(1,1,1)$. The condition and batch columns are no longer identical, so their main effects can in principle be distinguished. Estimability still depends on the design and on having adequate observations; the matrix should be checked before interpreting fitted coefficients.

Why Not Adjust Counts First?

A common mistake is to subtract an estimated batch effect from each raw count, or otherwise adjust counts to produce a “batch-corrected” expression matrix, then fit a negative-binomial count model to those adjusted values. This confuses two operations. A count likelihood describes the distribution of the observed count response; subtracting or continuously rescaling counts can produce non-integers, negative values, or a variance structure that is no longer the assumed negative-binomial one. The corrected values are not the original count observations required by that likelihood.

For inference about a condition effect with known batch, the direct approach is to retain the observed counts as responses and include batch in the model design, together with the normalization offset. Linear adjustment of transformed expression values can be useful for other purposes, such as visual summaries, but it is not equivalent to fitting the count model and should not be presented as if adjusted values were raw counts.

Model Checks and Limitations

A sensible fit begins with the design and offset, not only the final coefficient table. Check that sample labels match the count columns, each batch level is coded as intended, and the design matrix is full rank for the effects to be estimated. Confirm that the normalization factors are positive, attached to the correct samples, and on a consistent scale. Inspect whether normalized expression patterns still show a strong batch structure and whether the modeled mean and dispersion provide a plausible description of the counts.

Mis-coded or omitted batch: Verify batch labels against the sample records. An incorrect indicator or a missing batch term can leave technical structure in the model residuals or associate technical variation with a biological coefficient.
Redundant design columns: A rank-deficient model cannot uniquely estimate all listed coefficients. Check the design before interpreting condition and batch effects separately.
Inadequate normalization: A batch coefficient does not automatically repair a poor sample-scale factor. Examine normalization diagnostics and whether the offset has been supplied to the model as intended.
Overly simple mean structure: A single batch coefficient assumes a common batch log ratio across conditions. If that assumption is scientifically implausible, assess whether the design can support a more suitable model.

The model also has limits. It estimates batch-associated mean shifts under its chosen design and distributional assumptions; it does not guarantee that every technical effect has been removed. In a confounded design, the biological and batch effects may not be separately estimable. With few samples, dispersion and coefficient estimates may be uncertain. Finally, model-based condition comparisons depend on the specified normalization and covariates; they should be interpreted as conditional associations for the modeled samples and gene, not as universal or causal effects.

Key takeaway: Keep the integer count as the response, put the log normalization factor in the model as an offset, and represent known batch with design-matrix covariates. The offset sets the sample scale; the batch coefficient estimates a conditional batch association.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A gene has count 90 in a sample with normalization factor 1.5. What is its simply scaled value, and what count should remain the response in the count model?
In a model, the fitted batch coefficient is $\log(0.8)$. Interpret the batch B to reference-batch expected-count ratio, conditional on the other modeled terms.
Every control sample is in batch A and every treatment sample is in batch B. Can a model with separate condition and batch main effects identify both effects from these data? Why?
Why is subtracting a fitted batch value from each raw count not the same as including batch in a negative-binomial model?

Solutions. (1) The simply scaled value is $90/1.5=60$. The response remains the observed count 90; the factor enters as an offset. (2) The ratio is $\exp\{\log(0.8)\}=0.8$: batch B has 0.8 times the expected count, or 20% lower expected count, than the reference batch conditional on the other terms. (3) No. The condition and batch indicators are identical, so their coefficients cannot be separated; only their combined contribution is informed by the data. (4) Subtraction changes the observed response and can yield non-integer or negative values, which do not have the stated negative-binomial count interpretation. Including batch as a covariate retains the observed counts and estimates its association within the count model.