Negative Binomial Mean Models for Depth Normalization

Why Model Counts Against Depth?

Cells with greater sequencing depth often have larger observed counts for a gene, even when the gene’s underlying expression is similar across cells. A model can describe this relationship as a change in the expected count associated with depth, while allowing observed counts to vary around that expectation. This is different from dividing every count by a size factor: the model expresses the depth relationship on the count scale and makes its assumptions explicit.

Here $Y_{ig}$ denotes the observed count for gene $g$ in cell $i$, and $L_i$ is a positive relative depth factor for cell $i$. The factor may be derived from a normalization procedure; it represents relative sequencing depth, not an absolute number of reads. The model below describes one gene at a time. It does not by itself establish that any association with depth is purely technical.

The Gene-Wise Negative-Binomial Mean Model

Definition: For a given gene $g$, assume the count $Y_{ig}$ in cell $i$ has a negative-binomial distribution with mean $\mu_{ig}>0$ and dispersion $\alpha_g\geq 0$. Its variance is $$ \operatorname{Var}(Y_{ig}\mid L_i)=\mu_{ig}+\alpha_g\mu_{ig}^{2}. $$ The mean is modeled as $$ \log(\mu_{ig})=\beta_{0g}+\beta_{1g}\log(L_i), \qquad \mu_{ig}=\exp(\beta_{0g})L_i^{\beta_{1g}}. $$ Here $\beta_{0g}$ is the gene-specific intercept, $\beta_{1g}$ is the gene-specific depth coefficient, and $\alpha_g$ controls extra variation beyond the Poisson variance.

The logarithm requires $L_i>0$. The intercept $\beta_{0g}$ is the log expected count when $L_i=1$, so $\exp(\beta_{0g})$ is the expected count at that reference depth. If depth is multiplied by a factor $c$, the model’s expected count is multiplied by $c^{\beta_{1g}}$. Thus, $\beta_{1g}=1$ describes proportional scaling with depth, $\beta_{1g}=0$ describes a mean that does not change with depth, and other values describe a different power-law relationship.

The dispersion parameterization matters. When $\alpha_g>0$, one probability mass function is

P(Y_{ig}=y\mid L_i)= \frac{\Gamma(y+1/\alpha_g)}{\Gamma(1/\alpha_g)y!} \left(\frac{1}{1+\alpha_g\mu_{ig}}\right)^{1/\alpha_g} \left(\frac{\alpha_g\mu_{ig}}{1+\alpha_g\mu_{ig}}\right)^y, \qquad y=0,1,2,\ldots.

For $\alpha_g=0$, this mass function is not evaluated at zero; the model is instead defined as the Poisson distribution with mean $\mu_{ig}$, which is also the limiting distribution as $\alpha_g$ approaches zero. Positive dispersion gives variance greater than the mean. In single-cell data, that flexibility can be useful because counts often vary more than a Poisson model permits.

What the Fitted Mean Separates

Given fitted parameters $\hat{\beta}_{0g}$ and $\hat{\beta}_{1g}$, the fitted mean for cell $i$ is

\hat{\mu}_{ig}=\exp(\hat{\beta}_{0g})L_i^{\hat{\beta}_{1g}}.

This is the model’s expected count at that cell’s depth. The difference $Y_{ig}-\hat{\mu}_{ig}$ describes how far the observed count is above or below that expected value, in count units. It is a raw observed-minus-fitted difference, not a variance-standardized residual and not automatically a corrected expression value. The model attributes the systematic mean trend represented by $L_i$ to the fitted depth relationship; variation around that trend can reflect biological differences, sampling variation, omitted factors, or model failure.

A useful justification for the negative-binomial variance comes from its gamma-Poisson representation. For $\alpha_g>0$, suppose a latent rate $\Lambda_{ig}$ has mean $\mu_{ig}$ and variance $\alpha_g\mu_{ig}^2$, and, conditional on that rate, $Y_{ig}$ is Poisson with mean $\Lambda_{ig}$. The law of total variance gives

\begin{aligned} \operatorname{Var}(Y_{ig}\mid L_i) &=E\{\operatorname{Var}(Y_{ig}\mid \Lambda_{ig},L_i)\mid L_i\} +\operatorname{Var}\{E(Y_{ig}\mid \Lambda_{ig},L_i)\mid L_i\}\\ &=E(\Lambda_{ig}\mid L_i)+\operatorname{Var}(\Lambda_{ig}\mid L_i)\\ &=\mu_{ig}+\alpha_g\mu_{ig}^{2}. \end{aligned}

The first term is Poisson sampling variation; the second allows additional variability in the latent rate. This interpretation does not mean the extra variation has a single known biological source. It is a distributional allowance, not a causal decomposition of every observed deviation.

Worked Example: Calculating a Fitted Mean and Variance

Suppose a hypothetical gene has fitted parameters $\hat{\beta}_{0g}=\log(5)$, $\hat{\beta}_{1g}=0.6$, and $\hat{\alpha}_g=0.2$. A hypothetical cell has depth factor $L_i=3$. The fitted expected count is

\hat{\mu}_{ig}=\exp\{\log(5)\}\,3^{0.6} =5\exp\{0.6\log(3)\} \approx 5(1.933182)=9.666.

The fitted variance is therefore

\widehat{\operatorname{Var}}(Y_{ig}\mid L_i) =\hat{\mu}_{ig}+\hat{\alpha}_g\hat{\mu}_{ig}^{2} \approx 9.665910+0.2(9.665910)^2 \approx 28.35.

The intermediate unrounded mean is about $9.665910$; it rounds to $9.666$, not $9.665$. The variance exceeds the mean because the fitted dispersion is positive. These are model-based quantities, not claims that this particular cell must have a count near 9.666.

Fitting and Interpreting the Depth Relationship

For a gene, the parameters are typically estimated from its counts across cells, with the depth factors treated as known inputs. Maximum likelihood chooses parameters that make the observed counts plausible under the specified negative-binomial distributions. In a simple gene-wise formulation, counts are assumed conditionally independent across cells given their depth factors and the model parameters. The mean is log-linear in log depth, and the dispersion is shared across cells for that gene. A fitted coefficient is an association in the modeled mean; it is not proof that sequencing depth is the only cause of count differences.

The model is most interpretable when the depth factors are positive and reasonably reflect relative depth, the log-mean relationship is adequate over the observed range, and major systematic covariates are not omitted. In particular, cell type or state may be associated with both depth and expression. If so, a gene’s fitted depth slope can partly reflect biological composition rather than a purely technical response.

Worked Example: Separating the Fitted Trend from Count Differences

Consider hypothetical counts for one gene in three cells. Let the fitted model be $\hat{\mu}_{ig}=4L_i$, corresponding to $\hat{\beta}_{0g}=\log(4)$ and $\hat{\beta}_{1g}=1$. Suppose the data are:

Cell	Depth factor $L_i$	Observed count $Y_{ig}$	Fitted mean $\hat{\mu}_{ig}=4L_i$	Observed minus fitted
A	1	3	4	-1
B	2	11	8	3
C	4	13	16	-3

For cell B, for example, $\hat{\mu}_{Bg}=4(2)=8$, so $Y_{Bg}-\hat{\mu}_{Bg}=11-8=3$. The fitted means describe the model’s expected depth-related pattern: doubling depth doubles the expectation under the fitted slope of one. The differences show that the individual counts do not fall exactly on that pattern. They do not, on their own, identify why they differ.

Worked Example: Interpreting a Nonunit Depth Coefficient

Suppose a hypothetical gene has fitted depth coefficient $\hat{\beta}_{1g}=0.8$. If depth changes from $L_i=1$ to $L_i=2$, the expected-count multiplier is

\frac{\hat{\mu}_{g}(2)}{\hat{\mu}_{g}(1)} =\frac{\exp(\hat{\beta}_{0g})2^{0.8}}{\exp(\hat{\beta}_{0g})1^{0.8}} =2^{0.8} =\exp\{0.8\log(2)\} \approx 1.741.

Thus the fitted mean is about 1.741 times as large at twice the depth, an increase of about 74.1%, rather than a doubling. This describes the fitted mean relationship over these depths. It should not be interpreted as a causal effect of changing sequencing depth, particularly if depth is associated with biological characteristics of the cells.

Assumptions, Diagnostics, and Alternatives

Before relying on this model, check that every $L_i$ is positive and that the depth factors correspond to the intended comparison. Examine observed counts against depth for each gene or for representative genes: a pronounced curved pattern on the log-mean scale, a sharp change in slope, or systematic grouping by cell type can indicate that the mean model is inadequate or missing relevant structure. Very sparse genes may provide little information about a gene-specific slope. A fitted coefficient based on few nonzero counts can be unstable; this is a reason for caution, not evidence that the model has correctly removed depth effects.

The principal alternative in the earlier material on offset-based normalization is to use log depth as an offset, which fixes its coefficient at one. That specification assumes expected counts scale proportionally with depth. Estimating $\beta_{1g}$ instead allows the fitted mean to rise more or less than proportionally, but adds a gene-specific parameter and can mistake biological depth-associated patterns for technical effects. Simple size-factor scaling also imposes a proportional adjustment, but it does not model a negative-binomial count distribution. Neither flexibility nor a count distribution repairs inappropriate depth factors or confounding.

Common pitfall: “The fitted mean is the technical count, so everything above or below it is biological.” This does not follow. A fitted mean is an expectation under a chosen model; observed-minus-fitted differences combine biological variation, sampling noise, omitted influences, and possible misspecification. The model separates a specified depth-related mean trend from remaining variation, not technical causes from biological causes with certainty.

The model also has limits. It summarizes a gene’s mean-depth relationship and variance with a compact form; it cannot guarantee that all genes share an appropriate relationship, nor can it distinguish confounded technical and biological effects from counts alone. The scale of $L_i$ is relative: multiplying every depth factor by a constant changes the intercept but not the depth coefficient or the model’s fitted comparisons when the intercept is adjusted accordingly. Interpret fitted values only in the context of the depth factors and covariates used to estimate them.

Check Your Understanding

Try each problem, then compare your reasoning with the solutions.

If $\hat{\beta}_{1g}=1$, what happens to the fitted mean when depth is multiplied by three?
For $\hat{\mu}_{ig}=6$ and $\hat{\alpha}_g=0.1$, calculate the negative-binomial variance.
A fitted model gives an expected count of 12, while the observed count is 8. What is the raw observed-minus-fitted difference, and what does it not establish?
State one assumption that distinguishes a model with an estimated depth coefficient from an offset model.

Solutions. (1) The mean is multiplied by $3^1=3$, because the fitted model is proportional to depth. (2) $\operatorname{Var}(Y)=6+0.1(6^2)=6+3.6=9.6$ counts squared. (3) The difference is $8-12=-4$ counts; it does not identify the cause of the deviation or show that it is biological. (4) The estimated-coefficient model assumes a log-linear mean in log depth but estimates its slope from the gene’s counts; the offset model fixes that slope at one, imposing proportional scaling.