Why Model Batch?
Pseudobulk analysis aggregates cells into sample-level gene counts, so technical and processing differences between samples can affect the counts being compared. A batch may represent a processing run, sequencing lane, or other grouping that was not the biological condition of interest. If condition groups are distributed unevenly across batches, some apparent expression difference can reflect batch rather than condition. Including batch as a sample-level predictor lets the model estimate the condition contrast while accounting for systematic differences between batches.
Adjustment only works when the study design contains information to distinguish batch from condition. If every control sample is in one batch and every treated sample is in another, the data cannot determine which factor explains a difference. A model formula cannot repair that lack of overlap.
The Pseudobulk Model With Batch
For gene \(g\) in pseudobulk sample \(s\), let \(Y_{gs}\) be the count and let \(L_s\) be the sample's library-size offset on the log scale. A standard negative-binomial model is
Here \(\mu_{gs}\) is the expected count, \(\phi_g\) is the gene-specific dispersion, \(C_s\) indicates condition, and \(B_{ks}\) are indicator variables for the batch levels. One batch is the reference level, so its indicator is omitted. The variance is \(\mu_{gs}+\phi_g\mu_{gs}^2\). The library-size offset accounts for differences in sequencing depth; it is not a freely estimated batch effect. The expression above uses a condition indicator coded 0 for control and 1 for treated.
The model assumes that the negative-binomial mean and variance are suitable for the pseudobulk counts, that samples are the independent observational units (or that their dependence is handled by the model), and that the specified mean structure is reasonable. In particular, a model with batch main effects but no condition-by-batch interaction assumes a common condition log-fold change across batches. It also requires that the design matrix has enough independent information to estimate the requested coefficients. Dispersion estimation and library-size offsets follow the earlier material on pseudobulk models; the central issue here is how batch and condition are represented together.
For a batch with several levels, treatment coding compares each non-reference batch with the reference, at the same condition and offset. Alternative contrast coding changes the meanings of individual batch coefficients, but does not change the fitted comparisons or the condition contrast when the same model space is used. The condition coefficient is the target contrast only when that coefficient is identifiable under the chosen design.
Worked Examples of Adjustment and Identifiability
Worked Example: A Balanced Design Gives an Adjusted Fold Change
Suppose a hypothetical study has two batches and two conditions, with two independent pseudobulk samples per batch-condition combination. All samples have the same library-size offset. For one gene, the counts are:
| Batch | Control counts | Treated counts |
|---|---|---|
| A | 100, 120 | 150, 180 |
| B | 200, 240 | 300, 360 |
The within-batch treated-to-control ratios of the group means are \(165/110=1.5\) in batch A and \(330/220=1.5\) in batch B. The batch B to batch A ratios are \(220/110=2\) for controls and \(330/165=2\) for treated samples. Thus these data have exactly the common-ratio pattern represented by a log-linear model with batch and condition main effects.
For an explicit arithmetic illustration, consider the Poisson log-linear version of the model, with equal offsets. In this balanced example, the fitted mean for each batch-condition combination is its observed group mean. The condition coefficient is therefore
The adjusted fold change is \(\exp(0.4055)=1.500\): the expected count is 1.5 times as large in treated samples as in controls at a common batch and offset. The batch coefficient for B relative to A is \(\log(2)=0.6931\). This Poisson calculation makes the adjustment transparent; an actual pseudobulk analysis commonly uses a negative-binomial model and estimates dispersion as well. The numerical calculation here is not a claim that every negative-binomial fit will give exactly these coefficients.
Worked Example: Exact Confounding Prevents Separate Effects
Now suppose a hypothetical experiment has three control samples, all processed in batch A, and three treated samples, all processed in batch B. Let \(C_s=1\) for treated samples and \(B_s=1\) for batch B. Every sample then has \(B_s=C_s\). With an intercept, the condition and batch columns of the design matrix are identical.
For any gene, the linear predictor contains \(\beta_{Cg}C_s+\gamma_gB_s\). Since \(B_s=C_s\), this is \((\beta_{Cg}+\gamma_g)C_s\). The data can inform the sum, but not the two coefficients separately. For any number \(a\), replacing \(\beta_{Cg}\) by \(\beta_{Cg}+a\) and \(\gamma_g\) by \(\gamma_g-a\) leaves every fitted mean unchanged. Thus the condition effect is not identifiable independently of batch. A fitting program may report a dropped coefficient or a rank-deficiency warning; assigning a value to the condition coefficient does not make it estimable.
Worked Example: Near-Confounding Makes Adjustment Unstable
Consider a hypothetical design with four control samples in batch A, three treated samples in batch B, and one treated sample in batch A. There are no control samples in batch B. Unlike the previous example, the condition and batch columns are not identical: the treated sample in batch A supplies some overlap. The design can therefore have full rank, and a batch-adjusted condition contrast may be identifiable.
However, the comparison relies heavily on that single treated sample in batch A. If its gene count is unusual, or if it is removed for a justified reason, the estimated contrast may change substantially or become non-identifiable. The arithmetic check is about allocation, not the gene's count values: the small overlap is one sample out of eight, so there is little direct support for separating condition from batch. This is near-confounding, not exact confounding. Full rank is necessary for estimation, but does not guarantee a precise or robust estimate.
Why the Design Matrix Is the Key Check
The design matrix has one row per pseudobulk sample and one column per model coefficient. Its columns encode the intercept, condition, batch indicators, and any other included predictors. A coefficient is separately estimable only if its column contributes information that cannot be represented as a linear combination of the other columns. In particular, the batch-adjusted condition contrast cannot be estimated if the condition column is redundant with the batch columns and intercept.
This gives a short justification for the exact-confounding rule. If \(X\) is the design matrix and two different coefficient vectors produce the same \(X\beta\), then they produce the same linear predictors and fitted means. If \(Xv=0\) for some nonzero vector \(v\), then \(X(\beta+v)=X\beta\). The likelihood cannot distinguish those coefficient vectors, so the corresponding effects are not uniquely identified. When columns are nearly redundant, the estimates may still be unique, but small changes in the data can produce large changes in the coefficients.
Before fitting gene-level models, inspect a condition-by-batch table. Then check the rank of the proposed design matrix and look for sparse combinations or dependence on only one or two samples. These checks answer different questions: a full-rank matrix rules out exact linear dependence, while the allocation table reveals whether the estimate is practically supported across the design. A model fitting without an error is not proof that confounding is harmless.
Use biological sample counts, not the number of cells, to see which combinations are represented.
State the batch reference level and the condition comparison. Include other relevant sample-level predictors only when the design supports their estimation.
Check design-matrix rank and inspect whether the condition comparison is supported within batches, rather than by a single sparse combination.
Report the condition log-fold change or fold change as adjusted for the modeled batch terms, and note design limitations that affect its stability.
Comparing Adjustment Strategies and Common Pitfalls
| Approach | What it estimates or does | Main concern |
|---|---|---|
| Condition-only model | Estimates the condition contrast without modeling batch. | Can mix batch differences into the condition contrast when batch is associated with condition. |
| Condition plus batch in the count model | Estimates the condition contrast conditional on modeled batch effects. | Requires sufficient overlap and a suitable mean structure; it cannot resolve exact confounding. |
| Remove or transform batch effects before testing | Alters the expression values before a separate comparison. | May not preserve the count model, uncertainty, or intended contrast; it is not a substitute for including batch in the inferential model. |
A common mistake is to see that every treated sample is in batch B and every control sample is in batch A, then fit a model containing both factors and interpret the condition coefficient as adjusted. The correction is to recognize that the two factors are inseparable in those data: the fitted difference could be due to condition, batch, or both. The appropriate response is to report the limitation and, where possible, obtain a design with condition represented across batches. Do not interpret a software-selected reference, dropped column, or generalized inverse as scientific separation of the effects.
A second pitfall is to assume that adjustment always improves inference. Adding batch can account for systematic differences when the design supports it, but it also uses degrees of freedom and may increase uncertainty, especially with few samples or weak overlap. A third is to include an interaction without enough observations to estimate it. A condition-by-batch interaction allows the condition effect to differ by batch, but the common-effect interpretation no longer applies; sparse or absent condition groups within batches can make interaction estimates uninformative or unidentifiable.
Batch adjustment is most useful when batch is recorded, the condition groups have meaningful overlap across batches, and a common adjusted contrast matches the scientific question. It does not remove all technical variation, guarantee unbiased inference under omitted or mismeasured factors, or turn a confounded observational comparison into a randomized one. The unit of replication remains the biological sample represented by each pseudobulk, not the individual cells contributing to it.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- In a model with control coded 0 and treated coded 1, what does \(\exp(\beta_{Cg})\) represent when batch is included as a main effect?
- All controls are in batch A and all treated samples are in batch B. Can a model containing both condition and batch separate their effects? Why?
- A condition-by-batch table has at least one sample in every combination, but one combination contains only one sample. What does this suggest about identifiability and precision?
- Why is a condition-only model potentially misleading when condition is associated with batch?
- Does a full-rank design matrix guarantee that the adjusted condition estimate is stable? Explain.
Solutions. (1) It is the treated-to-control expected-count ratio at a fixed batch and library-size offset, under the specified model. (2) No. The condition and batch indicators are identical, so the data identify only a combined contrast, not separate effects. (3) The design may be identifiable if the columns are independent, but the sparse combination provides weak support and can yield an unstable estimate. (4) It omits batch, so batch-associated differences can contribute to the estimated condition contrast. (5) No. Full rank rules out exact linear dependence, but sparse overlap or dependence on a small number of samples can still make the estimate imprecise or sensitive.