Estimating The Dispersion-Mean Trend

Why Estimate a Dispersion-Mean Trend?

In RNA-seq count models, genes with similar average expression can still differ in variability. Yet a gene-wise dispersion estimate can be imprecise, especially when the number of biological replicates is small. Estimating how dispersion generally changes with mean expression across genes gives a reference for judging each gene’s variability. The reference is a cross-gene pattern; it does not establish that every gene follows the pattern exactly.

The earlier material on gene-wise dispersion defines the negative-binomial mean and dispersion for each gene. Here the focus is how to estimate a smooth relationship between those dispersions and the genes’ mean expression. The central distinction is between a gene’s noisy estimate and the overall trend across genes.

Definition: For gene $i$ in sample $j$, let $Y_{ij}$ be the observed count and $s_j$ the sample’s size factor. A negative-binomial model has mean $\mu_{ij}$ and dispersion $\alpha_i$, with $$ \operatorname{Var}(Y_{ij})=\mu_{ij}+\alpha_i\mu_{ij}^{2}. $$ The gene’s normalized mean, denoted $\bar{\mu}_i$, is typically calculated as the average of its fitted means on the normalized-count scale across samples. A dispersion-mean trend is a smooth function $f$ describing the general relationship between $\bar{\mu}_i$ and gene-wise dispersion estimates $\hat{\alpha}_i$.

A standard negative-binomial GLM uses a log link and an offset for the size factor, for example $\log(\mu_{ij})=\log(s_j)+\mathbf{x}_j^\mathsf{T}\boldsymbol{\beta}_i$. The offset places samples on a comparable scale; consequently, the trend should use fitted means on that normalized scale, not unadjusted average raw counts. The model also assumes that, conditional on its mean and dispersion, a gene’s counts follow the specified negative-binomial distribution. The trend is only as useful as these modeling and normalization assumptions.

Estimating the Trend Across Genes

A general estimation procedure has three ingredients: estimate a mean and a gene-wise dispersion for each gene, select a smooth functional form, and fit that form to the gene-level estimates. A common parametric choice is

f(\mu)=a_0+\frac{a_1}{\mu},

where the coefficients are chosen from the cross-gene data, subject to a fitted curve that gives nonnegative dispersions over the relevant range. This form allows dispersion to decline as mean expression increases and to approach a positive level $a_0$ at high means when $a_0>0$. It is a useful model, not a law: the data need not follow this exact curve.

A more flexible alternative estimates a smooth curve, for example by smoothing transformed dispersions against transformed means. A log scale can be helpful when both quantities span wide ranges, but it requires positive values; zero or boundary estimates need suitable handling rather than taking an undefined logarithm. Smoothers also require choices about flexibility and behavior near the ends of the observed range. Regardless of the method, the goal is to capture the broad pattern, not to trace every fluctuation in the gene-wise estimates.

Obtain gene-level summaries.
For each retained gene, fit the count model and record its normalized mean $\bar{\mu}_i$ and gene-wise dispersion estimate $\hat{\alpha}_i$.

Choose a smooth form.
Select a parametric curve or a flexible smoother appropriate to the range and scale of the estimates. State the form because different choices can produce different trends.

Fit and inspect the relationship.
Estimate the curve from genes across the expression range, then check that it follows the broad pattern without being driven by a few unusual estimates or extrapolated far beyond observed means.

One justification for smoothing is that gene-wise estimates contain sampling noise. Suppose, for a fixed mean $\mu$, an estimate $Z$ has conditional expectation $m(\mu)$. For any proposed value $c$, squared-error loss has conditional expectation

E\left[(Z-c)^2\mid \mu\right] =\operatorname{Var}(Z\mid\mu)+\left(m(\mu)-c\right)^2.

The variance term does not depend on $c$, so the expected squared error is minimized at $c=m(\mu)$. A smooth fit approximates this systematic conditional pattern using information across genes. This argument assumes the estimates have a meaningful conditional expectation and that squared error is an appropriate fitting criterion. It does not guarantee that an arbitrary chosen curve is correctly specified, nor that a single gene’s estimate equals the trend.

Worked Example: Calculating a Gene’s Normalized Mean

Consider one hypothetical gene measured in three samples with size factors $0.8$, $1.0$, and $1.2$, and observed counts $8$, $20$, and $36$. For this simple illustration, use count divided by size factor as the normalized value and average those values.

The normalized values are $8/0.8=10$, $20/1.0=20$, and $36/1.2=30$ counts per size-factor unit. Thus,

\bar{\mu}_i=\frac{10+20+30}{3}=20

The gene’s normalized mean is 20 counts per size-factor unit. Suppose its fitted gene-wise dispersion is $\hat{\alpha}_i=0.12$. The point $(20,0.12)$ is one input to a trend fit; it is not itself enough to determine the trend. In a full analysis, fitted means from the count model are generally used, particularly when covariates or condition effects make the sample means differ systematically.

A Small Parametric Trend Fit

The following hypothetical example shows the arithmetic for fitting the parametric form $f(\mu)=a_0+a_1/\mu$. For clarity, suppose three genes have the following normalized means and gene-wise dispersions. Real analyses would generally use many genes and consider the precision and distribution of their estimates.

Worked Example: Fitting a Reciprocal-Mean Curve

The three hypothetical gene summaries are $(\bar{\mu}_i,\hat{\alpha}_i)=(10,0.15)$, $(20,0.10)$, and $(50,0.07)$. Set $x_i=1/\bar{\mu}_i$ and $y_i=\hat{\alpha}_i$; then fit the straight line $y=a_0+a_1x$.

The transformed mean values are $0.10$, $0.05$, and $0.02$. Their averages are $\bar{x}=0.05667$ and $\bar{y}=0.10667$. The ordinary least-squares slope and intercept are

\hat{a}_1= \frac{\sum_i(x_i-\bar{x})(y_i-\bar{y})} {\sum_i(x_i-\bar{x})^2} =1,\qquad \hat{a}_0=\bar{y}-\hat{a}_1\bar{x}=0.05.

For example, the numerator is approximately $0.003267$, and the denominator is also approximately $0.003267$, giving the slope of 1. The fitted curve is therefore

\hat{f}(\mu)=0.05+\frac{1}{\mu}.

At a normalized mean of 25 counts per size-factor unit, the trend gives $\hat{f}(25)=0.05+1/25=0.09$. This is the fitted reference dispersion at that mean. This toy fit uses three points that lie exactly on the curve; real gene-wise estimates are noisy, so actual fitted values need not pass through every point.

The reciprocal-mean form also gives a useful interpretation of the count variance. At mean $\mu$, its variance is $\mu+\left(a_0+a_1/\mu\right)\mu^2$, or $\mu+(a_0\mu^2+a_1\mu)$. The curve therefore allows the amount of extra-Poisson variability to change with mean. The specific coefficients and functional form must still be estimated from the data; the variance expression does not imply that every RNA-seq data set should use this curve.

Worked Example: Interpreting a Gene Relative to the Trend

Suppose a fitted hypothetical trend is $\hat{f}(\mu)=0.05+1/\mu$. A gene has normalized mean $\bar{\mu}_i=40$ counts per size-factor unit and gene-wise dispersion estimate $\hat{\alpha}_i=0.12$.

The trend reference at this mean is

\hat{f}(40)=0.05+\frac{1}{40}=0.075.

The gene-wise estimate, 0.12, is above this reference. To see the scale of the difference, the negative-binomial variance at mean 40 using the trend dispersion is $40+0.075(40^2)=40+120=160$ count-squared units. Using the gene-wise estimate gives $40+0.12(40^2)=40+192=232$ count-squared units. These are model-based variances for counts at that mean, not variances in normalized-count units.

This comparison says that the gene’s estimated dispersion is higher than the fitted cross-gene reference at a similar mean. It does not, by itself, show that the gene is an outlier, establish a formal statistical difference, or identify a biological cause. The gene-wise estimate can be uncertain, and the trend itself is estimated.

Choices, Diagnostics, and Limitations

The trend depends on how the curve is fitted. A parametric form is compact and easy to evaluate, but may miss a genuine bend or plateau. A flexible smoother can follow more complex patterns, but can also follow noise if given too much flexibility. A common-dispersion model instead uses one dispersion reference for all genes, which can be inadequate when variability systematically changes with mean. At the other extreme, gene-wise estimates retain local differences but can be unstable when each gene has little information. These approaches answer different modeling needs; none is automatically best for every data set.

Approach	What it represents	Main limitation
Common dispersion	A single reference value across genes	Cannot represent a changing mean-variance relationship
Parametric trend	A smooth relationship with a specified functional form	Can be misleading if that form poorly matches the data
Flexible smoother	A smooth relationship with fewer fixed shape assumptions	Depends on smoothing choices and can behave poorly at boundaries
Gene-wise estimates alone	Separate estimates for individual genes	Can be noisy, especially with few replicates

Useful checks include plotting or tabulating gene-wise estimates against normalized means, checking whether the fitted curve remains nonnegative over the data range, and looking for systematic patterns left above or below the curve. A trend that changes sharply when a small number of genes are removed may be sensitive to those genes. Sparse data at low or high means also make the curve less certain there. The later material on reading a dispersion-mean plot develops visual interpretation; here the key point is to assess whether the fitted reference represents the broad pattern over the range where it will be used.

A common mistake is to treat the fitted trend as the dispersion of every gene. For example, a student might see a trend value of 0.075 at a gene’s mean and report that the gene’s dispersion is 0.075. That confuses the cross-gene reference with the gene-wise estimate. The correction is to keep both quantities distinct: the trend estimates the general relationship at that mean, while the gene-wise estimate describes that gene based on its own count data. A trend can be useful even when individual genes depart from it, but it cannot repair a badly specified count model or compensate for inadequate normalization.

Key takeaway: Estimate the dispersion-mean trend from gene-wise summaries on the normalized expression scale. Treat it as a smooth cross-gene reference, check its adequacy over the observed range, and do not interpret it as the exact dispersion of each gene.

Check Your Understanding

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

If a gene’s fitted means across three samples are 12, 18, and 30 normalized counts, what is its average fitted mean?
For the trend $f(\mu)=0.04+0.8/\mu$, calculate the reference dispersion at a normalized mean of 20.
At mean 20, compare the negative-binomial variances implied by dispersions 0.08 and 0.12. State the units.
Give one advantage and one limitation of a flexible smoother relative to a parametric curve.
Why is an observed gene-wise dispersion above the trend not, by itself, proof of an unusual gene?

Solutions. (1) The mean is $(12+18+30)/3=20$ normalized counts. (2) $f(20)=0.04+0.8/20=0.08$. (3) With $\operatorname{Var}(Y)=\mu+\alpha\mu^2$, the variances are $20+0.08(400)=52$ and $20+0.12(400)=68$ count-squared units. (4) A flexible smoother can represent a pattern that does not match a chosen parametric shape, but its result depends on smoothing choices and may follow noise. (5) Both the gene-wise estimate and trend are estimated with uncertainty; being above the reference alone is not a formal test and does not establish that the gene is an outlier.

Estimating The Dispersion-Mean Trend

What You'll Learn