Common Errors In Count-Based Differential Expression

Why Count-Based Differential Expression Can Go Wrong

A differential-expression result is only as reliable as the model that connects observed counts to the biological comparison. In RNA-seq, a large count can reflect higher expression, greater sequencing exposure, or both. Counts also tend to be more variable than a simple Poisson model predicts, and testing thousands of genes creates a multiplicity problem. A procedure that overlooks any one of these features can produce compelling-looking but misleading results.

This tutorial focuses on diagnosing those errors. The central model for gene $g$ in sample $i$ is a negative-binomial generalized linear model (NB GLM), which represents counts directly, includes an exposure offset, and allows gene-specific extra variation. It is not the only valid approach, but it makes the assumptions behind common mistakes especially clear.

Definition: Let $Y_{gi}$ be the observed count for gene $g$ in sample $i$, and let $\mu_{gi}=E(Y_{gi})$. In the NB GLM convention used here, $Y_{gi}$ has variance $\mu_{gi}+\alpha_g\mu_{gi}^2$, where $\alpha_g\geq 0$ is the dispersion. A typical mean model is $$ \log(\mu_{gi})=\log(e_i)+\mathbf{x}_i^{T}\boldsymbol{\beta}_g. $$ Here $e_i>0$ is an effective exposure or library-size factor, $\mathbf{x}_i$ is the design vector for sample $i$, and $\boldsymbol{\beta}_g$ contains gene-specific coefficients. The term $\log(e_i)$ is an offset: its coefficient is fixed at one rather than estimated.

This model assumes that the design describes the relevant experimental structure, the samples are independent conditional on the model (unless dependence is explicitly modeled), and the mean is correctly linked to the design and exposure. Dispersion must also be estimated or otherwise accounted for. An offset is not a substitute for a valid design, and a negative-binomial distribution does not repair confounding or incorrect sample labels.

Error 1: Treating Raw Counts as Ordinary Continuous Outcomes

An ordinary linear model for a gene typically treats the outcome as a continuous response with an additive mean and errors that are independent, approximately normal, and—in the usual equal-variance version—have the same variance across groups. Raw counts are nonnegative integers, often skewed at low means, and their variance generally changes with their mean. More importantly, counts depend on sequencing exposure. Comparing raw group means without accounting for exposure can mistake a sequencing-depth difference for a biological difference.

A log link with an offset encodes a multiplicative relationship: expected count equals exposure times a modeled expression rate. For a two-group model, if the treatment coefficient is $\beta_g$, then the fitted treatment-to-control ratio of expected counts per unit exposure is $\exp(\beta_g)$, assuming the other modeled covariates are held fixed. This is a relative expression comparison, not automatically a comparison of absolute molecules per cell.

Worked Example: A Raw-Count Difference Caused by Exposure

Suppose four hypothetical control samples have counts $90,100,110,100$, and four treatment samples have counts $180,200,220,200$ for one gene. Assume, for this illustration, that control exposure factors are $e_C=1$ and treatment exposure factors are $e_T=2$. The values are deliberately constructed so that treatment counts are exactly twice control counts at matched positions.

The group means are $\bar{x}_C=(90+100+110+100)/4=100$ counts and $\bar{x}_T=(180+200+220+200)/4=200$ counts. A raw-count linear model therefore estimates a treatment-minus-control difference of $200-100=100$ counts. The sample variances are $s_C^2=[(-10)^2+0^2+10^2+0^2]/3=66.67$ and $s_T^2=[(-20)^2+0^2+20^2+0^2]/3=266.7$. A Welch-style standard error for this difference is $$ SE=\sqrt{\frac{66.67}{4}+\frac{266.7}{4}}=\sqrt{83.34}=9.129. $$ The raw difference is large relative to this standard error, but that comparison does not account for exposure. Dividing each count by its stated exposure gives control rates $90,100,110,100$ and treatment rates $90,100,110,100$. Both mean rates equal $100$ counts per exposure unit. Under the hypothetical exposure information, there is no difference in the rate represented by these data. The raw-count contrast answers a different question and is not evidence of a treatment effect on exposure-adjusted expression.

In real analyses, exposure factors are estimated from the data or supplied by an appropriate measurement design; they are not chosen to force group means to agree. This example checks the arithmetic and shows why raw counts alone do not identify a biological contrast.

A transformed-count analysis can sometimes support a linear-model workflow, but it requires its own variance treatment and assumptions. The point is not that every use of a linear model is invalid; it is that an untransformed raw-count linear model with ordinary constant-variance errors is not automatically appropriate. The earlier material on choosing a test for differential expression distinguishes count-model approaches, including likelihood-ratio, quasi-likelihood, and exact tests, according to the design and assumptions.

Error 2: Treating Total Counts as an Unquestionable Offset

An offset must represent the exposure relevant to the expression comparison. Total library counts are a possible proxy, not a guaranteed truth. RNA-seq counts are compositional: the total reads are distributed among the features sequenced, so a change in a small number of highly abundant genes can alter total counts and the apparent fractions of many other genes. Without external calibration, a global shift in RNA abundance can also be difficult to distinguish from a change in sequencing exposure.

One common alternative estimates sample factors from the pattern across genes, relying on an assumption that most genes are not changing in a strongly one-sided way. Such factors can be less influenced by a few extreme genes than total-count scaling, but they are not assumption-free. The relevant diagnostic question is whether the chosen factors sensibly represent exposure for the biological comparison, not whether one normalization rule is always correct.

Worked Example: One Abundant Gene Changes the Total

Consider two hypothetical samples and three genes. Sample A has counts $100,100,800$, for a total of $1000$; sample B has counts $100,100,1800$, for a total of $2000$. Suppose the third gene is unusually abundant in sample B, while the first two are the majority of the listed genes and have unchanged counts.

If total counts are used directly as exposure factors, the first gene's count per total is $100/1000=0.10$ in A and $100/2000=0.05$ in B. The log rate contrast for B versus A is $$ \log(0.05)-\log(0.10)=\log(0.5)=-0.6931. $$ That calculation implies a twofold decrease in the first gene's normalized rate. It is a consequence of total-count scaling in this example, even though its observed count is $100$ in both samples.

For comparison, consider a median-ratio size-factor calculation for these three genes. Their geometric means across samples are $100$, $100$, and $\sqrt{800\cdot1800}=1200$. Sample A's count-to-geometric-mean ratios are $1,1,800/1200=0.6667$, whose median is $1$. Sample B's ratios are $1,1,1800/1200=1.5$, also with median $1$. These factors leave the first gene's normalized count unchanged between samples. This comparison does not prove the median-ratio factors are correct: they rely on a majority-stable pattern and the toy data contain only three genes. It shows why checking how an abundant feature affects the proposed offsets matters.

The observed counts alone cannot determine whether the third gene's increase reflects more molecules, a change in composition, or another exposure-related feature. If absolute RNA abundance is the target, an appropriate external calibration or experimental measurement may be needed.

Error 3: Ignoring Dispersion

The Poisson model sets the count variance equal to its mean. Biological replication commonly produces extra variability, so the NB variance $\mu_{gi}+\alpha_g\mu_{gi}^{2}$ is often more suitable. When $\alpha_g>0$, the variance exceeds the mean, especially at higher means. Ignoring this extra variation typically makes standard errors too small and tests too optimistic. The size of the effect depends on the mean, the dispersion, the design, and the number of replicates.

A dispersion-mean plot and residual diagnostics help reveal problems, but no single plot proves that a model is correct. The earlier material on reading a dispersion-mean plot and residual diagnostics for negative-binomial GLMs gives tools for checking whether the fitted variability and residual patterns are plausible. Extreme observations and influential samples also deserve investigation rather than automatic deletion.

Worked Example: How Dispersion Changes a Standard Error

Suppose two hypothetical groups each have $n=4$ independent samples, with observed mean counts $20$ and $30$ for a gene. For a transparent approximation, take a common NB dispersion $\alpha=0.20$, and use a Wald comparison of the difference in group means. This simplified calculation assumes the group means are estimated independently and uses the NB variance evaluated at each mean; it is illustrative rather than a replacement for fitting the full design.

Under an NB model, the variances at the two means are $20+0.20(20^2)=100$ and $30+0.20(30^2)=210$. The estimated difference is $30-20=10$, and its approximate standard error is $$ SE_{\mathrm{NB}}=\sqrt{\frac{100}{4}+\frac{210}{4}}=\sqrt{77.5}=8.803. $$ Thus the Wald statistic is $10/8.803=1.136$, corresponding to an approximate two-sided standard-normal p-value of $0.256$.

If one instead imposes the Poisson variance-equals-mean assumption, the standard error would be $$ SE_{\mathrm{Pois}}=\sqrt{\frac{20}{4}+\frac{30}{4}}=\sqrt{12.5}=3.536, $$ giving $10/3.536=2.828$ and an approximate two-sided p-value of $0.0047$. The data summary is the same; the uncertainty changes because the variance assumption changes. The two p-values are approximations under their respective models, not competing facts about the gene. In actual analysis, dispersion is estimated using the experimental design and available replication, and its uncertainty may also be relevant.

With very few replicates, dispersion estimates for individual genes can be unstable. Methods that share information across genes can improve stability, but they do not eliminate the need to inspect fit and design. Quasi-likelihood methods offer a way to account for uncertainty in dispersion estimates; the earlier material on quasi-likelihood tests covers that approach.

Error 4: Interpreting Unadjusted Gene-Wise Tests as a Discovery List

A gene-wise p-value is calculated for a particular null hypothesis, such as no condition effect for that gene. If thousands of such hypotheses are tested, some small p-values can arise by chance even when all null hypotheses are true. A nominal threshold such as $0.05$ does not control the expected proportion of false discoveries among a large list of selected genes. For large-scale testing, report an appropriate multiplicity-adjusted result, commonly a false-discovery-rate adjusted p-value, and state the threshold used.

Worked Example: Nominal Significance Does Not Ensure an FDR Discovery

Suppose $100$ hypothetical genes are tested. Ten have nominal p-values of $0.01$, and the other $90$ have p-values of $0.50$. At a nominal $0.05$ threshold, the first ten would all be called significant. Apply the Benjamini-Hochberg adjustment at rank $r$ among $m=100$ tests by first calculating $p_{(r)}m/r$, then enforcing monotonicity from the largest rank toward the smallest.

At rank 10, the raw adjustment is $0.01(100)/10=0.10$. For ranks 1 through 9, the raw adjustment values are $0.01(100)/r$, all at least $0.10$; after the monotonicity step, the adjusted p-values for those ten tests are also $0.10$. The remaining tests have adjusted p-values $0.50$. Thus none passes an adjusted threshold of $0.05$, even though ten have nominal p-values of $0.01$. The adjusted p-value is not the probability that a particular gene is null; it is a multiplicity-adjusted measure used for controlling false discoveries under the method's assumptions.

Multiplicity correction does not repair a bad count model, poor normalization, confounding, or an unsuitable gene-wise hypothesis. Conversely, an adjusted p-value above a chosen cutoff does not establish that there is no biological effect; it means the analysis did not identify that test as a discovery at that threshold. Effect sizes and their uncertainty should be considered alongside adjusted p-values. The earlier material on interpreting adjusted p-values in RNA-seq develops this distinction.

Checks, Comparisons, and Failure Modes

Before interpreting a differential-expression table, verify that the analysis matches the experimental question and that the data support its assumptions. A practical review asks whether the design includes relevant conditions and covariates, whether exposure factors are defensible, whether dispersion and residual behavior are plausible, and whether the reported discovery threshold accounts for the number of tests. The following comparison summarizes the central distinctions.

Approach or shortcut	What it assumes or does	Diagnostic question
Raw-count ordinary linear model	Uses counts as continuous outcomes, often with normal, constant-variance errors; may omit exposure.	Are scale, variance, and exposure assumptions credible for these counts?
Count GLM with offset	Models counts with a specified distribution and includes log exposure in the mean model.	Do the exposure factors represent the relevant library scaling, and does the fitted model describe the data?
Poisson count model	Sets variance equal to the mean.	Is there evidence of extra variation beyond the Poisson variance?
Negative-binomial count model	Allows variance $\mu+\alpha\mu^2$.	Are dispersion estimates and residual diagnostics plausible for the design?
Nominal gene-wise p-value threshold	Assesses each test without controlling multiplicity across the gene list.	How many hypotheses were tested, and what adjusted criterion supports the discovery list?

Failure signals include systematic residual patterns, extreme observations driving a contrast, dispersion estimates that do not fit the mean-variance pattern, strong sensitivity to plausible normalization choices, and a large gap between nominal and adjusted findings. These signals call for investigation: check sample identities and design coding, inspect the counts and diagnostics, and compare a method suited to the design. Do not select whichever analysis yields the smallest p-value. Methods such as exact tests for two groups, likelihood-ratio tests for nested designs, and quasi-likelihood tests answer related questions under different structures; their assumptions should guide the choice.

These methods also have limitations. Normalization factors cannot, by themselves, identify absolute RNA abundance from purely compositional sequencing data. Low counts and small replication limit information, while model-based uncertainty depends on the mean-variance assumptions and dispersion estimation. Multiplicity adjustment controls a specified error criterion; it does not guarantee that every selected gene is biologically important or that every real effect will be detected. A defensible report gives the design, normalization approach, test, effect estimate, uncertainty or adjusted evidence, and the relevant limitations.

Common Mistake

Common mistake: “The treatment group has larger raw counts, so the gene is upregulated; total library counts are the correct offset; the Poisson test is significant; therefore the gene is a discovery.” This reasoning treats counts as exposure-free, assumes away overdispersion, and ignores the family of tests. The correction is to check the experimental design and effective exposure, use a count model whose variance is appropriate, examine diagnostics, and base a discovery claim on an adjusted criterion while interpreting the effect in context.

A useful discipline is to separate the four questions: what is the count-generating model, what exposure is being adjusted for, how much variation is expected across biological replicates, and what error rate applies to the set of tests? A persuasive answer to only one of these does not settle the others.

Check Your Understanding

Use the calculations and distinctions above to answer these questions. Brief solutions follow each problem.

A gene has mean count $50$ in a group and estimated dispersion $\alpha=0.10$. What variance does the NB model assign, and how does it compare with the Poisson variance?
In the exposure example, what was the raw treatment-minus-control mean difference, and what was the difference in mean count per exposure unit?
Why can a total-count offset make an unchanged observed count appear to have a lower normalized rate?
In the 100-test example, how many genes pass an adjusted threshold of $0.05$, and why does this differ from the nominal result?

Solutions. (1) The NB variance is $50+0.10(50^2)=50+250=300$, compared with Poisson variance $50$; the NB model allows substantial extra variation. (2) The raw difference was $200-100=100$ counts, while the mean rates were $200/2=100$ and $100/1=100$, for a difference of $0$ counts per exposure unit. (3) The total count can increase because other genes contribute more reads, so dividing the same target count by a larger total yields a smaller fraction; that does not by itself establish a decrease in absolute abundance. (4) None pass: the ten nominal p-values of $0.01$ have BH-adjusted p-values $0.10$, which exceed $0.05$; nominal testing did not account for the 100 simultaneous tests.