Empirical Bayes Shrinkage of Gene Dispersions

Why Shrink Gene-Wise Dispersions?

A pseudobulk negative binomial model uses dispersion to represent variability beyond the Poisson mean–variance relationship. Estimating that dispersion separately for every gene is appealing because genes need not have the same variability. But with only a few biological samples, a gene’s counts provide limited information about its own dispersion. The resulting estimate can be very noisy, sometimes extremely large or close to zero.

Empirical Bayes shrinkage addresses this instability by combining gene-specific evidence with information shared across genes. The shared information is summarized by a mean–dispersion trend: genes with similar average abundance are expected to have dispersions in a similar range. Shrinkage does not assert that all genes have the same dispersion. It moderates estimates according to how much information each gene provides and how much gene-to-gene variation is supported by the data.

Model and Shrinkage Assumptions

For gene \(g\) in biological sample \(j\), let \(Y_{gj}\) denote the pseudobulk count and \(\mu_{gj}\) its modelled mean. Under the negative binomial parameterization used here,

The dispersion \(\alpha_g\) is positive. When \(\alpha_g=0\), the variance reduces to the Poisson variance \(\mu_{gj}\); larger values allow greater variation among biological samples at a given mean. The means are specified by the negative binomial generalized linear model, including the design and library-size offsets. Dispersion shrinkage concerns estimation of \(\alpha_g\), not the definition of the mean model.

The rationale relies on several assumptions. First, the negative binomial mean model and its sampling structure must be adequate for the pseudobulk data; a trend cannot repair an incorrect design or unmodelled dependence. Second, after accounting for average abundance, many genes must provide useful information about the distribution of dispersions. This is a pooling assumption, not a claim that genes are identical. Third, the trend must represent the typical abundance–dispersion relationship well enough that it is a reasonable prior centre. Finally, gene-wise uncertainty estimates must reflect how little information the available biological samples provide.

A convenient description uses the log dispersion \(\eta_g=\log(\alpha_g)\), since dispersion is positive and often varies multiplicatively. Let \(m_g=\log\{\alpha_0(\bar{\mu}_g)\}\). A normal approximation to the gene-wise estimate is

Here \(\hat{\eta}_g\) is the log of a gene-wise dispersion estimate, \(s_g\) is its approximate standard error, and \(\tau^2\) describes variation in true log dispersions around the trend. The notation \(\dot{\sim}\) signals an approximation: practical methods may instead use the negative binomial likelihood or a profile-likelihood approximation directly. The normal model makes the shrinkage logic explicit.

Deriving the Shrinkage Estimate

Under the normal approximation, the likelihood for \(\eta_g\) is proportional to \(\exp\{-(\hat{\eta}_g-\eta_g)^2/(2s_g^2)\}\), and the prior density is proportional to \(\exp\{-(\eta_g-m_g)^2/(2\tau^2)\}\). Multiplying them adds their quadratic terms. Completing the square gives a normal posterior whose mean and mode are

The weights are nonnegative and sum to one. When \(s_g^2\) is large relative to \(\tau^2\), the data estimate is imprecise and the trend receives more weight. When \(s_g^2\) is small, the gene-wise estimate receives more weight. The shrunken dispersion is \(\tilde{\alpha}_g=\exp(\tilde{\eta}_g)\); on the original scale this is a weighted geometric, not arithmetic, average of the gene-wise estimate and trend.

The trend and \(\tau^2\) are estimated empirically across genes. In broad terms, a method fits the abundance–dispersion relationship, assesses how much gene-to-gene spread remains around it after allowing for gene-wise sampling error, and uses that residual spread to set the prior strength. If estimated spread is large, shrinkage is weak; if it is small, estimates are pulled more strongly toward the trend. The details vary by implementation. Some methods use a maximum a posteriori estimate from a likelihood with a prior penalty rather than the normal-approximation posterior mean shown here.

Estimating and Checking the Prior

In an analysis, the trend and prior spread are not known in advance. They are estimated from the collection of genes, usually after low-expression genes have been filtered and the mean model has been fitted. The abundance trend can be estimated with a smooth curve or another specified function. Prior spread is then assessed from the dispersion of gene-wise estimates around that trend, accounting for their differing uncertainty. These estimates should be treated as modelling choices to inspect, not as automatic truths.

Useful diagnostics include plotting gene-wise dispersions and shrunken dispersions against average abundance, with the fitted trend visible. A reasonable result usually shows noisy gene-wise values moderated toward the trend, while well-supported departures can remain. Inspect whether the residual spread around the trend changes systematically with abundance, whether many genes pile up at estimation boundaries, and whether a small set of extreme genes unduly influences the trend or prior spread. Compare results under reasonable trend specifications when the relationship is visibly irregular.

A fitted prior with almost no residual spread can force nearly every gene to the trend, even when distinct biological subgroups of genes are plausible. Conversely, an overly broad prior produces little regularization and leaves unstable estimates largely untouched. Very small sample sizes can also make gene-wise standard errors and the prior fit themselves uncertain. Shrinkage reduces variance in estimates; it does not create biological replicates or make the original information limitation disappear.

Alternatives, Uses, and Limitations

A common mistake is to say, “The gene-wise dispersion is \(0.40\), so the gene is definitely overdispersed relative to other genes.” That reasoning ignores uncertainty: with few samples, the estimate may be far from the gene’s underlying dispersion simply because the data are weak. The correction is to compare the estimate with its uncertainty and the appropriate abundance trend. A large raw estimate with a large standard error may shrink substantially; a modest departure supported by precise data may remain.

The main limitation is dependence on the shared structure. If the trend is misspecified, if genes form groups with systematically different dispersion patterns not represented by abundance, or if outlying genes distort hyperparameter estimation, shrinkage can pull estimates in an inappropriate direction. Robust fitting can reduce sensitivity to outliers, but it cannot guarantee that the trend is scientifically suitable. As with any regularization, the resulting estimate reflects both the observed counts and the assumptions used to stabilize them.