Worked Example With Two Conditions

A Two-Condition Comparison From Counts to Interpretation

A two-condition RNA-seq analysis asks whether a gene's expression differs between two groups while accounting for the fact that samples can have different sequencing depths. The central steps are connected: normalize counts for sample-specific library scale, represent the counts with a model that allows overdispersion, test a clearly defined condition contrast, and report the effect in a direction that matches the chosen reference group.

The small data below are hypothetical and chosen to make the arithmetic transparent. They illustrate the analysis rather than provide evidence about any real biological system. In practice, dispersion estimation, quality checks, and genome-wide multiplicity adjustment are essential parts of interpreting results.

Data, Size Factors, and Model

Let $K_{ij}$ be the observed count for gene $i$ in sample $j$, and let $s_j$ be that sample's size factor. A size factor represents relative sequencing depth; it is not a gene-specific adjustment. One common median-ratio estimate first calculates each gene's geometric mean across samples, then takes the median of the count-to-geometric-mean ratios within each sample. Genes with zero counts in a way that makes their geometric mean unusable are omitted from this calculation. The resulting factors are commonly scaled to have geometric mean one.

Definition: For a two-condition comparison with condition A as the reference and $x_j=1$ for condition B (otherwise $x_j=0$), the negative binomial generalized linear model for gene $i$ is $$ K_{ij}\sim \operatorname{NB}(\mu_{ij},\alpha_i),\qquad \operatorname{Var}(K_{ij})=\mu_{ij}+\alpha_i\mu_{ij}^{2}, $$ $$ \log(\mu_{ij})=\log(s_j)+\beta_{0i}+\beta_{1i}x_j. $$ Here $\log(s_j)$ is an offset, $\alpha_i$ is the gene's dispersion, $\beta_{0i}$ is the log mean for condition A at size factor one, and $\beta_{1i}$ is the log ratio of condition B to condition A. Thus the B-versus-A fold change is $\exp(\beta_{1i})$, and the log base 2 fold change is $\beta_{1i}/\log(2)$.

The model assumes condition labels and offsets are specified correctly, observations are independent conditional on the model, and the negative binomial mean-variance relationship is an adequate description for the gene. The offset adjusts for library scale; it does not correct batch effects or other omitted design variables. The dispersion must also be estimated or supplied appropriately. In a real analysis, dispersion estimation usually borrows information across genes, while allowing gene-specific differences.

For the condition effect, the usual two-sided Wald test uses $H_0:\beta_{1i}=0$ against $H_a:\beta_{1i}\ne 0$. If $\widehat{\beta}_{1i}$ has estimated standard error $SE$, the Wald statistic is $z=\widehat{\beta}_{1i}/SE$; its large-sample reference distribution under the null is standard normal. This tests a zero log fold change, equivalently a fold change of one.

Worked Example: Calculating Size Factors

Suppose four hypothetical samples, ordered A1, A2, B1, B2, have the following counts for three genes.

Gene	A1	A2	B1	B2
G1	50	100	100	200
G2	20	40	40	80
G3	10	20	40	80

For G1, the geometric mean is $(50\cdot100\cdot100\cdot200)^{1/4}=100$, giving ratios $0.5,1,1,2$. G2 has geometric mean $40$ and gives the same ratios. For G3, the geometric mean is $(10\cdot20\cdot40\cdot80)^{1/4}\approx28.28$, giving ratios approximately $0.354,0.707,1.414,2.828$.

Take the median ratio separately within each sample. The four medians are $0.5,1,1,2$: for example, the A1 ratios are $0.5,0.5,0.354$, whose median is $0.5$. Their geometric mean is $(0.5\cdot1\cdot1\cdot2)^{1/4}=1$, so the factors are already scaled to geometric mean one. The size factors are therefore $s=(0.5,1,1,2)$.

Dividing counts by their sample factors gives normalized counts. G1 becomes $100,100,100,100$; G2 becomes $40,40,40,40$; and G3 becomes $20,20,40,40$. The last pattern suggests a twofold B-versus-A difference, which the model will quantify while accounting for count variability.

Fitting and Testing the Condition Contrast

The condition coefficient is a log ratio, not a difference in raw counts. With condition A coded $x=0$ and condition B coded $x=1$, a positive $\beta_{1i}$ means higher modeled expression in B; a negative value means lower modeled expression in B. Reversing the coding reverses the sign, so a reported fold change should always name its numerator and denominator.

A useful justification for the standard error comes from the negative binomial likelihood. For a log-linked mean, the expected information about a sample's log mean is $\mu/(1+\alpha\mu)$. With a separate mean in each condition and independent samples, the approximate variance of the difference between estimated log means is the sum of the inverse total information in each group:

\operatorname{Var}(\widehat{\beta}_1)\ \approx\ \frac{1}{\sum_{j\in A}\mu_j/(1+\alpha\mu_j)} + \frac{1}{\sum_{j\in B}\mu_j/(1+\alpha\mu_j)}.

This expression shows why dispersion and replication matter. Larger means add information, but when $\alpha>0$, information does not grow without bound as the mean increases. More independent samples can increase the group information directly. Software commonly calculates the Wald standard error from the fitted model, using the fitted means and the dispersion estimate.

Worked Example: Fitting G3 and Testing B Versus A

Use the four G3 counts from the table, the size factors $0.5,1,1,2$, and the model with A as reference. For a transparent hand calculation, suppose the dispersion is treated as known to be $\alpha=0.1$. The normalized counts are $10/0.5=20$, $20/1=20$, $40/1=40$, and $80/2=40$. Thus the fitted condition means are $\widehat{\mu}_A=20$ and $\widehat{\mu}_B=40$, matching the equal normalized observations within each group.

The fitted log contrast is $$ \widehat{\beta}_1=\log(40)-\log(20)=\log(2)\approx0.6931. $$ The estimated B-versus-A fold change is $\exp(0.6931)=2.000$, or a log base 2 fold change of $0.6931/\log(2)=1.000$.

For each A sample, the information contribution is $20/(1+0.1\cdot20)=20/3$. The two A samples contribute $40/3$, so the A mean's approximate variance on the log scale is $1/(40/3)=0.07500$. For each B sample, the contribution is $40/(1+0.1\cdot40)=8$. The two B samples contribute $16$, giving variance $1/16=0.06250$. Therefore, $$ SE=\sqrt{0.07500+0.06250}=\sqrt{0.13750}\approx0.3708,\qquad z=0.6931/0.3708\approx1.869. $$

The two-sided standard-normal p-value is approximately $0.0617$. At $\alpha=0.05$, this illustrative gene-level Wald test does not reject $H_0$. The estimated expression is twofold higher in B, but the uncertainty is substantial for only two samples per condition under the assumed dispersion. The p-value is not the probability that the null hypothesis is true. In a real analysis, dispersion would be estimated, and the gene-level p-value would need to be considered alongside multiple-testing adjustment if many genes were tested.

Worked Example: A Contrast With the Opposite Direction

Using the same G3 fit, suppose the scientific question is expressed as A versus B rather than B versus A. The contrast is now $\beta_A-\beta_B=-\widehat{\beta}_1=-0.6931$. Its standard error remains $0.3708$, because reversing a contrast changes its sign but not its variance. The Wald statistic is $-0.6931/0.3708=-1.869$, with the same two-sided p-value, approximately $0.0617$.

Exponentiating the reversed contrast gives $\exp(-0.6931)=0.5000$. This says the modeled mean in A is one-half the modeled mean in B, equivalently that B is twofold higher than A. Both descriptions are consistent; the reference and numerator must be stated to avoid an apparent contradiction. A two-sided test has the same p-value in either direction, while a one-sided alternative would need to be specified before looking at the result.

Worked Example: Why the Offset Changes the Conclusion

For G1, the raw mean in A is $(50+100)/2=75$ counts, while the raw mean in B is $(100+200)/2=150$ counts. Comparing raw means alone gives $150/75=2$, which might be misread as a twofold condition effect.

After dividing by the size factors, G1 has normalized counts $100,100$ in A and $100,100$ in B. Its fitted condition contrast is $\log(100/100)=0$, and its fold change is $\exp(0)=1$. The raw counts doubled because the B samples had larger size factors, not because G1's normalized expression differed. This example shows why a condition comparison must include the offset rather than use raw group means as the effect estimate.

Assumptions, Diagnostics, and Alternatives

A model-based result is only as credible as its inputs and fit. Before interpreting a coefficient, check that the samples are assigned to the intended conditions, the size factors are plausible, and the design can separate condition from any known batch or pairing structure. The simple model above is appropriate only when a two-group comparison is the relevant design and observations are independent. If samples are paired or a batch effect must be modeled, that structure belongs in the design rather than being ignored.

Residual diagnostics for negative binomial GLMs, dispersion-mean plots, and checks for extreme observations help reveal poor fit. Warning signs include residual patterns by condition, a gene with a count far beyond its fitted expectation, or dispersions inconsistent with the assumed mean-variance behavior. A few genes with unusual counts can also affect normalization if they dominate the distribution of ratios. Diagnostics do not prove that a model is correct, but they can identify cases where the stated assumptions are implausible and the result needs investigation.

The Wald test is convenient when the coefficient estimate and standard error are reliable. A likelihood-ratio test compares a full model with a nested model lacking the condition term; it can be useful when the question is framed as a comparison of model fit. Exact two-group tests provide another approach for particular two-group count settings. Quasi-likelihood procedures account for uncertainty in dispersion estimation in their testing framework. These procedures are not interchangeable in every design, and the choice should reflect the model, replication, and inferential target rather than whichever p-value is smallest.

The estimates also have limitations. Size factors correct relative depth under the normalization method's assumptions; they do not establish absolute RNA abundance per cell. A small sample can yield a large fold-change estimate with considerable uncertainty. The p-value describes compatibility with a null contrast under the model, not biological importance. Report the direction, fold change, uncertainty or test result, and units of interpretation (a ratio of normalized expression, not a change in raw read counts). For genome-wide conclusions, account for the many tested genes using an appropriate multiple-testing procedure.

Common Mistakes

Using raw group means as the condition effect. In the G1 example, that produces a spurious twofold ratio. Use a model with the size-factor offset.
Forgetting the contrast direction. A B-versus-A fold change of 2 is an A-versus-B fold change of 0.5. Name both conditions in the interpretation.
Interpreting a large fold change as proof of differential expression. Effect size and uncertainty are distinct; the G3 example has a twofold estimate but a Wald p-value above 0.05.
Treating a gene-level p-value as a genome-wide conclusion. When many genes are tested, multiplicity must be addressed.
Assuming normalization fixes every design problem. A size factor cannot account for omitted batch structure, dependence, or a poorly fitting count model.

Key takeaway: In a two-condition negative binomial GLM, the condition coefficient is a log ratio of normalized means, and its exponent is the fold change in the stated direction. The offset, dispersion, design, and uncertainty are all part of that interpretation.

Check Your Understanding

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

If the fitted B-versus-A coefficient is $-0.4055$, what is the approximate B-versus-A fold change?
Why can a gene have higher raw counts in B but no condition effect after fitting the offset?
What happens to the sign, standard error, and two-sided Wald p-value when the contrast direction is reversed?
A gene has a two-sided p-value of $0.03$ in a genome-wide analysis. Does this alone establish that it is significant after adjustment for multiple testing?

Solutions. 1. $\exp(-0.4055)\approx0.6667$, so the modeled mean in B is about two-thirds the mean in A. 2. B may have larger size factors, so its raw counts can be higher even when counts adjusted for relative library scale have equal means. 3. The sign reverses, the standard error is unchanged, and the two-sided p-value is unchanged. 4. No. The adjusted result depends on the set of tested genes and the multiple-testing procedure; an unadjusted p-value alone is insufficient.