Negative Binomial GLMs for Pseudobulk Counts

Why Use a Negative Binomial Model for Pseudobulk Counts?

Pseudobulk analysis sums cell-level counts within a biological sample and cell type, producing a count for each gene in each sample. These counts are nonnegative integers, but their variability across biological samples is often greater than a Poisson model allows. A negative binomial generalized linear model (GLM) accommodates that extra variability while relating each gene’s expected count to the experimental design.

The observational units for the model are the pseudobulk samples, not the individual cells that contributed to each sum. Thus, the rows of the model correspond to biological samples (or to the relevant sample-level observations), and the response for a given gene is its pseudobulk count across those rows. The model does not make the cells within a sample into independent biological replicates.

The mean and dispersion describe different features. The mean sets the expected count for a sample under the model; dispersion determines how much variability there is around that mean beyond ordinary Poisson sampling variation. A gene can have a large expected count and low relative extra variation, or a modest expected count and substantial extra variation. Dispersion is not the same as the standard deviation: the model-implied standard deviation is \(\sqrt{\mu_{ig}+\alpha_g\mu_{ig}^{2}}\).

Specifying the Gene-Wise GLM

Let \(X\) be the design matrix, with one row per pseudobulk sample and columns encoding the predictors in the experimental design. For gene \(g\), let \(\boldsymbol{\beta}_g\) be its coefficient vector, and let \(\mathbf{x}_i^{T}\) be row \(i\) of \(X\). A basic negative binomial GLM uses a log link:

The linear predictor \(\eta_{ig}\) is the design-row weighted sum of coefficients. Exponentiating it gives the fitted mean count, \(\mu_{ig}=\exp(\eta_{ig})\), which is necessarily positive. The coefficients and design columns are gene-specific in their effects: a condition or covariate can have a different relationship with expected count for each gene.

For example, with an intercept and a binary condition indicator \(T_i\), the design row is \(\mathbf{x}_i^T=(1,T_i)\), and the predictor is \(\eta_{ig}=\beta_{0g}+\beta_{1g}T_i\). When \(T_i=0\), the mean is \(\exp(\beta_{0g})\); when \(T_i=1\), it is \(\exp(\beta_{0g}+\beta_{1g})\). This illustrates how the design matrix maps the experimental groups to mean counts. The choice and interpretation of design columns follow the study design; this model specification does not by itself determine a testing procedure.

The equation above is the basic specification without an exposure or library-size offset. It is appropriate as written only when the mean structure is suitable for the count scale being modeled—for example, when relevant exposure differences are absent or have otherwise been addressed. In practical pseudobulk analyses, library sizes can differ substantially; incorporating an offset is a separate part of the model specification. A model without a needed adjustment can mistake differences in sequencing depth for differences in gene expression.

What the Dispersion Does

One useful representation of the negative binomial distribution is a Poisson model with a sample-specific latent rate that varies across a gamma distribution. Conditional on a latent rate \(\Lambda_{ig}\), let \(Y_{ig}\mid\Lambda_{ig}\sim\operatorname{Poisson}(\Lambda_{ig})\), and let \(\Lambda_{ig}\sim\operatorname{Gamma}(1/\alpha_g,\alpha_g\mu_{ig})\), using shape and scale parameters. This mixture gives the negative binomial variance.

This mixture is a useful explanation of the distribution, not a claim that a particular biological mechanism must generate the variation. In the GLM, the dispersion describes the chosen count distribution’s mean–variance relationship. It does not, by itself, identify the source of excess variability or remove dependence among repeated measurements from the same donor.

Distributional Assumptions and a Probability Calculation

For clarity, a common parameterization sets \(r=1/\alpha_g\) when \(\alpha_g>0\). Under that parameterization, the probability of observing a count \(y\) is

This form has mean \(\mu_{ig}\) and variance \(\mu_{ig}+\alpha_g\mu_{ig}^{2}\). In particular, setting \(y=0\) gives \(P(Y_{ig}=0)=(1+\alpha_g\mu_{ig})^{-1/\alpha_g}\). This can help assess whether a specified mean and dispersion allow plausible zero counts, although a zero frequency alone is not a complete model check.

The basic gene-wise model relies on several conditions. The response values must be counts; the log-mean relationship must represent the systematic effects encoded in the design; and, conditional on the modeled predictors, sample-level observations must satisfy the assumed independence structure. The negative binomial variance function must also be an adequate description of variability. These are modeling assumptions, not consequences of aggregating cells.

If multiple pseudobulk observations come from the same donor, they may be correlated. Treating them as independent rows in an ordinary GLM can understate uncertainty. Donor-aware fixed- or random-effect structures may be needed for repeated or paired observations. Similarly, an omitted relevant sample-level predictor, severe design confounding, or an unmodeled library-size difference can make the fitted mean structure inadequate even if the negative binomial variance formula is flexible.

Diagnostics, Alternatives, and Limits

Diagnostics should evaluate both the mean and the distributional variation. Useful checks include comparing observed and fitted counts across samples, examining residual patterns against fitted means and design variables, and checking whether the frequency of zeros and large counts is broadly compatible with fitted distributions. Patterns of residuals by condition or batch can signal a missing mean effect; residual spread that systematically increases or decreases with the fitted mean can signal a poor variance specification. A small number of extreme counts can also have substantial influence and should be checked against the underlying sample and count-processing records.

The Poisson GLM is a simpler alternative, but it assumes \(\operatorname{Var}(Y_{ig})=\mu_{ig}\). If biological samples vary more than this assumption permits, Poisson-based uncertainty can be too small. A negative binomial GLM relaxes that restriction using \(\alpha_g\). Conversely, the negative binomial model is not automatically correct merely because it permits extra-Poisson variation: it can still fail if the mean is misspecified, observations are dependent, or the assumed distribution does not describe the data adequately.

A further practical limitation is that gene-specific dispersion can be difficult to estimate reliably when there are few biological replicates. The model’s definition specifies the role of dispersion but does not solve the estimation problem. Nor does a gene-wise GLM alone address multiplicity across genes, select a scientifically meaningful contrast, or establish causal effects. Those questions require additional design and inferential choices.

The key modeling sequence is to define the sample-level response, encode the design in \(X\), specify the expected count through the log link, and state the dispersion-dependent variance. Every fitted mean and uncertainty calculation depends on these choices being compatible with the biological replication and sample structure.