Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Contrasts and Log-Fold Changes

Differential Expression With Pseudobulk and Mixed Models · Tutorial 484 of 700

Contrasts and Log-Fold Changes

You will learn to express a specified condition comparison as a linear contrast of fitted coefficients and convert its estimate to a fold change on the desired log scale.

Advanced 10 min read

What You'll Learn

  • Write fitted condition comparisons as linear combinations of model coefficients.
  • Convert a natural-log contrast to a log2 fold change and a fold ratio.
  • Build contrasts for multi-level factors and interaction models.
  • Check whether a contrast is estimable under the fitted design.
  • Diagnose common errors involving offsets, coefficient order, and reference levels.

Why Contrasts Matter

A fitted model produces coefficients, but a scientific question usually asks about a comparison: how does expression in one condition differ from expression in another? A coefficient may answer that question directly under one coding scheme, or it may represent only part of the comparison. A contrast provides a precise way to combine coefficients into the quantity actually being asked about.

For pseudobulk differential expression, the comparison is made for a gene in a specified cell type, using donor-level pseudobulk observations and a fitted model. This tutorial focuses on constructing the comparison and calculating its log-fold change. The value depends on the model parameterization, the conditions compared, and—in models with interactions—the covariate values at which the comparison is defined.

Model and Contrast Definition

For a gene, a common pseudobulk model uses a negative binomial response \(Y_i\) for sample \(i\), with mean \(\mu_i\) and dispersion \(\alpha\). Its linear predictor is

$$ Y_i \sim \operatorname{NB}(\mu_i,\alpha), \qquad \log(\mu_i)=\log(s_i)+\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta}. $$

Here \(s_i\) is the library-size normalization factor used as an offset, \(\mathbf{x}_i\) is the design row, and \(\boldsymbol{\beta}\) is the coefficient vector. The offset accounts for the sample's sequencing exposure; it is not an estimated condition effect. The model's coefficients and linear predictor use the natural logarithm.

To compare conditions while holding other design variables fixed, compare their fitted log means after removing the offset. If \(\mathbf{x}_1\) and \(\mathbf{x}_0\) are the design rows for the two conditions at the same specified covariate values, then the log ratio of normalized expected expression is \((\mathbf{x}_1-\mathbf{x}_0)^\mathsf{T}\boldsymbol{\beta}\). The offset cancels because the comparison uses the same exposure.

Definition: A linear contrast of the coefficient vector is a quantity \(c^\mathsf{T}\boldsymbol{\beta}\), where \(c\) is a fixed vector of weights. For two condition predictions with design rows \(\mathbf{x}_1\) and \(\mathbf{x}_0\), the comparison contrast is \(c=\mathbf{x}_1-\mathbf{x}_0\). Its estimate is \(c^\mathsf{T}\hat{\boldsymbol{\beta}}\). With a natural-log link, this estimate is the natural-log fold change; dividing it by \(\log(2)\) gives the log2 fold change.

This calculation assumes the fitted model's log-mean specification is appropriate for the comparison, the design columns have the stated coding, and the requested contrast is estimable. A contrast is estimable when the data and design identify that linear combination, even if individual coefficients are not uniquely determined. In a full-rank design, the usual coefficient contrasts are estimable. If a design is rank-deficient, the comparison must be checked rather than inferred from an arbitrary coefficient solution.

From a Contrast to a Fold Change

Let \(\Delta\) be the estimated contrast on the natural-log scale. Exponentiating gives the estimated ratio of normalized expression for condition 1 relative to condition 0. The log2 fold change is \(\Delta/\log(2)\). Thus the two quantities answer related but distinct presentation needs: the log2 fold change is additive on a base-2 log scale, while the fold ratio is multiplicative on the original scale.

$$ \widehat{\Delta} =c^\mathsf{T}\hat{\boldsymbol{\beta}}, \qquad \widehat{\log_2 FC} =\frac{\widehat{\Delta}}{\log(2)}, \qquad \widehat{FC} =\exp(\widehat{\Delta}). $$

A positive contrast means the numerator condition has higher fitted normalized expression; a negative contrast means it has lower fitted normalized expression. Reversing the comparison changes the sign of the log-fold change and takes the reciprocal of the fold ratio.

Worked Examples

Worked Example: Two Conditions With Treatment Coding

Suppose a hypothetical pseudobulk model for one gene uses an intercept and a condition indicator: control is coded 0 and stimulated is coded 1. The fitted coefficients are \(\hat{\beta}_0=2.10\) and \(\hat{\beta}_1=0.70\), with a natural-log link. Calculate stimulated relative to control.

The design rows are \(\mathbf{x}_{control}=(1,0)\) and \(\mathbf{x}_{stimulated}=(1,1)\). Therefore \(c=\mathbf{x}_{stimulated}-\mathbf{x}_{control}=(0,1)\), so the contrast is \(0\hat{\beta}_0+1\hat{\beta}_1=0.70\). The intercept cancels.

The log2 fold change is \(0.70/\log(2)=0.70/0.6931=1.010\). The fold ratio is \(\exp(0.70)=2.014\). In this hypothetical comparison, fitted normalized expression under stimulation is about 2.01 times the control level, or about 1.01 log2 units higher. The interpretation is for the modeled gene and cell type, not for raw counts with unequal library sizes.

Worked Example: Comparing Two Nonreference Levels

Suppose condition has three levels: A is the reference, and the design contains indicators for B and C. For a hypothetical gene, the fitted coefficients are \(\hat{\beta}_0=1.40\), \(\hat{\beta}_B=0.30\), and \(\hat{\beta}_C=-0.20\). Find the comparison C relative to B.

The fitted linear predictor for B is \(\beta_0+\beta_B\); for C it is \(\beta_0+\beta_C\). Subtracting gives \((\beta_0+\beta_C)-(\beta_0+\beta_B)=\beta_C-\beta_B\). The contrast weights in coefficient order \((\beta_0,\beta_B,\beta_C)\) are therefore \(c=(0,-1,1)\).

Substitution gives \(\widehat{\Delta}=0(1.40)-1(0.30)+1(-0.20)=-0.50\). Hence \(\widehat{\log_2 FC}=-0.50/0.6931=-0.7213\), and \(\widehat{FC}=\exp(-0.50)=0.6065\). Fitted normalized expression in C is about 0.607 times the level in B, or about 39.3% lower. The reference level A is part of the coefficient coding but does not appear in this particular comparison.

Worked Example: A Condition Effect With an Interaction

Consider a hypothetical model with treatment coded 0 for control and 1 for treated, and a cell-state indicator coded 0 for state X and 1 for state Y. In coefficient order \((\beta_0,\beta_T,\beta_Y,\beta_{TY})\), suppose the estimates are \((1.00,0.40,0.20,0.30)\). Calculate the treated-versus-control log2 fold change separately for each state.

In state X, the state indicator is 0. The control predictor is \(\beta_0\), and the treated predictor is \(\beta_0+\beta_T\). Their difference is \(\beta_T=0.40\). Thus \(\log_2 FC=0.40/0.6931=0.5771\), and the fold ratio is \(\exp(0.40)=1.492\).

In state Y, the state indicator is 1. The control predictor is \(\beta_0+\beta_Y\); the treated predictor is \(\beta_0+\beta_Y+\beta_T+\beta_{TY}\). The difference is \(\beta_T+\beta_{TY}=0.40+0.30=0.70\). Thus \(\log_2 FC=0.70/0.6931=1.010\), and the fold ratio is \(\exp(0.70)=2.014\). The treatment coefficient alone describes the treatment effect in the reference state X; it is not the treatment effect in state Y.

Deriving the Comparison and Checking It

The contrast formula follows directly from subtracting the model's two linear predictors. For conditions 1 and 0 at the same offset and covariate settings, the normalized log means are \(\mathbf{x}_1^\mathsf{T}\boldsymbol{\beta}\) and \(\mathbf{x}_0^\mathsf{T}\boldsymbol{\beta}\). Therefore,

$$ \log\left(\frac{\mu_1/s}{\mu_0/s}\right) =\log(\mu_1)-\log(\mu_0) =(\mathbf{x}_1-\mathbf{x}_0)^\mathsf{T}\boldsymbol{\beta}. $$

Exponentiating both sides yields the ratio of normalized expected expression, \(\exp\{(\mathbf{x}_1-\mathbf{x}_0)^\mathsf{T}\boldsymbol{\beta}\}\). This establishes why subtracting design rows, rather than choosing a coefficient by its label, is the general construction. If two predictions have different offsets, subtracting their raw fitted log means also subtracts \(\log(s_1)-\log(s_0)\). That is a ratio of expected counts at different exposures, not solely the condition-associated normalized expression comparison.

A reliable calculation can be organized as follows:

1
Specify the comparison.
Name the numerator and denominator conditions, gene, cell type, and any covariate values at which the comparison is intended.
2
Construct the design rows.
Apply the fitted model's actual factor coding and interaction terms to each condition. Keep the coefficient order explicit.
3
Subtract and evaluate.
Set \(c=\mathbf{x}_1-\mathbf{x}_0\), then calculate \(c^\mathsf{T}\hat{\boldsymbol{\beta}}\). Confirm that the contrast is estimable.
4
Convert and report.
Divide by \(\log(2)\) for log2 fold change, or exponentiate the natural-log contrast for the fold ratio. State the comparison direction and units of the log scale.

A useful arithmetic check is to calculate both fitted normalized log means and subtract them. This should agree with the contrast calculation. For a multi-level factor, releveling the reference category changes the coefficient coordinates but should not change the fitted comparison when the equivalent contrast is constructed correctly.

Diagnostics, Alternatives, and Limitations

Several checks help identify a mistaken contrast. First, verify the coefficient order against the model matrix; a weight vector applied to a different order can produce a plausible but wrong number. Second, compute the numerator and denominator linear predictors independently and check their difference. Third, reverse the comparison: the contrast should change sign, while the fold ratio should become its reciprocal. In an interaction model, explicitly state the covariate values because a condition effect can differ across groups.

An offset mistake is especially consequential. The offset describes sequencing exposure, so the biological comparison should use fitted normalized expression at a common exposure, or subtract the offset from each fitted log mean. Comparing raw expected counts from samples with different offsets can mix library size with expression. Also check that the fitted model and normalization are appropriate for the intended pseudobulk estimand; the contrast does not repair a misspecified design or an unsuitable aggregation.

The direct-coefficient approach is convenient when a coefficient exactly matches the requested comparison, as in a two-condition treatment-coded model. A general contrast is safer for nonreference comparisons, sum coding, and interactions. Alternatively, one can calculate two fitted linear predictors and subtract them; this is algebraically equivalent and often easier to audit. None of these approaches changes the fitted model—they are different ways to extract a specified comparison from it.

A coefficient contrast is conditional on the model's covariates and design. If an interaction is present, there is no single condition effect independent of the interacting variable unless a particular value or a defined summary is specified. Nor does a log-fold change by itself establish statistical significance, precision, or a causal effect. Those require additional inferential assumptions and evidence. For now, the essential task is to define and calculate the intended estimand correctly.

Common Mistakes

  • Using a coefficient without checking its meaning. A treatment coefficient in a model with an interaction is usually the effect at the reference level, not a universal effect. Construct the design rows for the desired comparison.
  • Reversing numerator and denominator. Control relative to treated is the negative of treated relative to control on the log scale. State the direction before interpreting the sign.
  • Calling a natural-log contrast a log2 fold change. Divide by \(\log(2)\) to report log2 units; exponentiate the natural-log value to obtain a fold ratio.
  • Including unequal offsets in the biological contrast. A ratio of raw expected counts can reflect library-size differences. Compare normalized expected expression at common exposure.
  • Ignoring coding and estimability. A coefficient label alone does not establish that the requested comparison is identified. Inspect the design coding and verify the contrast.
Key takeaway: Define the comparison first, form the difference between its design rows, and apply that contrast to the fitted coefficients. This gives the natural-log fold change; divide by \(\log(2)\) for log2 fold change or exponentiate for a fold ratio.

Check Your Understanding

Work each comparison from the stated coefficient order and coding. Brief solutions are included for checking.

  1. A two-level model has coefficients \((\beta_0,\beta_T)=(1.5,-0.35)\), with control coded 0 and treatment coded 1. What is the treatment-versus-control log2 fold change and fold ratio?
  2. A three-level factor uses A as reference, with coefficient order \((\beta_0,\beta_B,\beta_C)\). Write the contrast for B relative to C.
  3. In an interaction model, the treatment effect at reference group X is \(\beta_T\), and the interaction coefficient for group Y is \(\beta_{TY}\). If their estimates are 0.25 and -0.10, calculate the treatment effect in Y on the natural-log and log2 scales.
  4. Why should a contrast of raw fitted log counts with unequal offsets not automatically be reported as a normalized expression log-fold change?

Solutions. (1) The contrast is \(-0.35\). The log2 fold change is \(-0.35/0.6931=-0.5050\), and the fold ratio is \(\exp(-0.35)=0.7047\). Treatment is estimated at about 0.705 times control normalized expression. (2) For B relative to C, \((\beta_0+\beta_B)-(\beta_0+\beta_C)=\beta_B-\beta_C\), so the contrast weights are \((0,1,-1)\). (3) The effect in Y is \(\beta_T+\beta_{TY}=0.25-0.10=0.15\). Its log2 value is \(0.15/0.6931=0.2164\). (4) The raw fitted log-count difference contains the offset difference \(\log(s_1)-\log(s_0)\); it can therefore reflect unequal sequencing exposure as well as expression. Remove or equalize the offsets for the normalized comparison.

Clinical Biostats