Worked Example With Three Conditions

From Three Conditions to a Planned Comparison

With three conditions, a gene can differ in several ways: one condition may differ from both others, or only one pair may differ. A useful analysis begins by stating the scientific comparison of interest, then encoding all three conditions in a model that can estimate it. The model uses the same negative binomial generalized linear model framework as other count-based differential-expression analyses; the emphasis here is on how the design matrix and a planned pairwise contrast work together.

For gene \(g\) in sample \(j\), let \(Y_{gj}\) be the observed count, \(s_j\) the sample’s normalization factor, and \(x_j\) the row of the design matrix. A standard model is

Here \(\alpha_g\) is the gene-specific dispersion and \(\log(s_j)\) is an offset, not an estimated condition effect. The model assumes that counts are conditionally independent across biological samples given the model, that the negative binomial mean-variance relationship is adequate, and that the design includes the relevant systematic effects. The log link makes coefficients differences on the log-mean scale; exponentiating such a difference gives a mean ratio, or fold change, conditional on the model.

Design Matrix for Three Conditions

Let A be the reference condition, with B and C represented by indicator columns. In a model containing only condition, the design matrix has an intercept and two condition indicators. Its rows depend on each sample’s condition, not on its observed count.

The phrase “at offset zero” describes the model’s normalized mean scale: for a sample with size factor \(s_j\), the expected count is \(s_j\) times the modeled normalized mean. The condition comparisons are unaffected by this common offset because the same factor is accounted for in each sample’s model.

Planned Contrasts and Their Uncertainty

A contrast is a specified linear combination of model coefficients. For the B-to-C comparison under reference coding, define \(c=(0,1,-1)^\mathsf{T}\). The estimated log mean ratio is \(\hat{\theta}_g=c^\mathsf{T}\hat{\beta}_g=\hat{\beta}_{Bg}-\hat{\beta}_{Cg}\). Its standard error depends on the estimated covariance matrix of the coefficients, not just on their separate standard errors.

For a planned two-sided Wald test of \(H_0:\theta_g=0\) against \(H_a:\theta_g\ne 0\), the statistic is \(z=\hat{\theta}_g/SE(\hat{\theta}_g)\). A large-sample normal reference distribution gives the approximate p-value. An approximate 95% confidence interval on the log scale is \(\hat{\theta}_g\pm1.96SE(\hat{\theta}_g)\). Exponentiating the estimate and both interval endpoints gives a fold-change estimate and interval. These approximations rely on the fitted coefficient estimates being sufficiently well behaved for the Wald approximation; small samples, weakly informed coefficients, or poor model fit can undermine that approximation.

Justifying the Contrast Calculation

The variance formula follows from the variance of a linear combination. Write the estimation error in the contrast as \(c^\mathsf{T}(\hat{\beta}_g-\beta_g)\). If the coefficient estimation error has covariance matrix \(V_g\), then for any fixed vector \(c\),

For \(c=(0,1,-1)^\mathsf{T}\), expanding this product gives \(\operatorname{Var}(\hat{\beta}_{Bg})+\operatorname{Var}(\hat{\beta}_{Cg})-2\operatorname{Cov}(\hat{\beta}_{Bg},\hat{\beta}_{Cg})\). Omitting the covariance would therefore generally give the wrong standard error. The sign matters: the covariance is subtracted twice for a difference and added twice for a sum.

Adding Covariates, Diagnostics, and Alternatives

The same contrast logic applies when the design includes covariates. For example, if each condition is represented in both of two batches, a model can include an additional batch indicator. In sample order A-batch1, A-batch2, B-batch1, B-batch2, C-batch1, C-batch2, one suitable design matrix is

where the fourth column marks batch 2. In this additive model the adjusted B-to-C contrast is still \((0,1,-1,0)^\mathsf{T}\): the batch coefficient cancels. This requires the design to distinguish condition from batch. If, for example, all A samples occur in batch 1 and all B and C samples in batch 2, some condition and batch effects may be inseparable. A rank-deficient design cannot uniquely estimate all those coefficients; additional software output cannot repair a confounded experiment.

Check that the matrix encodes the actual sample labels and that each condition has biological replication. Inspect residual diagnostics and the dispersion-mean relationship for indications that the count model is inadequate, and investigate extreme observations rather than assuming they are harmless. A single influential sample can change a gene’s estimate or uncertainty substantially. Also verify that the contrast direction matches the question: reversing B-to-C to C-to-B reverses the log estimate and reciprocates the fold change, while leaving a two-sided test unchanged.

A planned pairwise contrast is not the same as an omnibus test of whether any of the three conditions differ. The pairwise test addresses one specified comparison; an omnibus test assesses a broader condition effect, typically using a joint test of the two condition coefficients. Separate pairwise tests can address several comparisons, but the number of tests and the multiplicity strategy must be considered. In RNA-seq, thousands of genes are typically tested, so gene-level p-values should not be treated as if only one hypothesis had been examined; adjusted p-values address that across-gene multiplicity. A contrast also cannot compensate for an inappropriate model, too few biological replicates, or poor experimental balance.

Common Pitfalls

• Reading a coefficient as the desired comparison. With A as reference, \(\beta_{Bg}\) compares B with A, not B with C. Subtract \(\beta_{Cg}\) to obtain B versus C.
• Ignoring coefficient covariance. The variance of a difference contains \(-2\operatorname{Cov}(\hat{\beta}_{Bg},\hat{\beta}_{Cg})\). Adding the two variances alone is generally incorrect.
• Calling a log coefficient a fold change. A coefficient of \(-0.2232\) is a log mean ratio. Exponentiating it gives the ratio, about \(0.8000\).
• Interpreting a non-significant result as equality. A p-value above the chosen significance level does not establish equal expression; assess the estimate and its uncertainty in context.
• Using an unidentifiable design. If condition and batch are completely confounded, the model cannot separate their effects. Balance and replication are design requirements, not post-analysis options.