Tutorials › Statistical Genetics › Covariate-Adjusted Regression for Gene-Set Enrichment

Gene Prioritization and Pathway Tests · Tutorial 685 of 989

Covariate-Adjusted Regression for Gene-Set Enrichment

Use a gene-level regression to estimate and test the association between pathway membership and gene-level association scores, conditional on prespecified covariates.

Advanced 12 min read

What You'll Learn

  • Specify the gene-level outcome, pathway indicator, and adjustment covariates.
  • Interpret the pathway coefficient as a covariate-adjusted score difference.
  • Derive the coefficient by residualizing membership and scores against covariates.
  • Calculate and interpret a t test and confidence interval for the pathway coefficient.
  • Recognize assumptions and reasons ordinary regression standard errors may be unreliable.
  • Distinguish a competitive regression question from a claim that a pathway causes disease.

Why Adjusted Regression Is Useful

A pathway may contain genes with stronger observed association scores for reasons that have little to do with pathway biology. For example, longer genes or genes with more tested variants may have more opportunities to accumulate strong variant-level evidence. A regression can ask whether pathway membership predicts a gene-level association score after accounting for prespecified gene characteristics.

The analysis uses one row per gene in a defined gene universe. Its outcome is the gene-level score, and its main predictor indicates whether the gene belongs to the pathway being tested. The resulting comparison is competitive: it asks whether pathway genes have higher or lower scores than genes outside the pathway, conditional on the included covariates. It does not, by itself, establish that pathway activity causes the trait or that every gene in the pathway is associated.

Definition: For genes \(g=1,\ldots,n\), let \(Y_g\) be a gene-level association score, \(S_g=1\) if gene \(g\) belongs to the pathway under study and \(S_g=0\) otherwise, and \(C_{g1},\ldots,C_{gq}\) be prespecified gene-level covariates. A covariate-adjusted linear regression is $$ Y_g=\beta_0+\beta_S S_g+\sum_{j=1}^{q}\gamma_j C_{gj}+\varepsilon_g. $$ The coefficient \(\beta_S\) is the adjusted difference in mean score between pathway-member and nonmember genes at the same values of the included covariates, under this linear mean model.

The outcome might be a gene-level score such as a transformed association p-value, where larger values represent stronger evidence. Its scale and direction must be stated: a positive \(\beta_S\) means higher scores for members, not necessarily larger effects in a particular biological direction. Gene-level scores can be constructed in different ways, so the score definition and the gene universe should be fixed before fitting the model.

Covariates may include gene length, variant count, or other prespecified characteristics that could relate to score magnitude or measurement opportunity. Their inclusion defines the conditional comparison; it does not automatically remove every source of bias. Choose them based on the analysis question and design, rather than searching among many covariate sets for the most favorable pathway result.

Model, Hypothesis, and Interpretation

For a single pathway, the usual regression test concerns the pathway coefficient:

$$ H_0:\beta_S=0 \qquad\text{versus}\qquad H_a:\beta_S\ne 0. $$

A directional alternative can be used if a direction was specified in advance, but a two-sided test is appropriate when either higher or lower scores are of interest. In the ordinary linear model, the estimated pathway coefficient \(\hat{\beta}_S\) is tested using \(t=\hat{\beta}_S/SE(\hat{\beta}_S)\). Under the classical model assumptions, the null reference distribution is a t distribution with residual degrees of freedom \(df=n-r\), where \(r\) is the number of fitted regression coefficients, including the intercept.

For the classical t test and its corresponding confidence interval to have their stated exact reference distributions, the model must have a correctly specified linear conditional mean, errors with mean zero that are independent and have a common variance, normally distributed errors, and no perfect multicollinearity. The design must also contain both pathway members and nonmembers with enough covariate overlap to estimate an adjusted contrast. Without overlap, the fitted coefficient may depend heavily on extrapolation.

In genetic data, gene scores are not always independent: nearby genes may share association evidence, and scores may have unequal variances. Thus, the ordinary regression coefficient can still be a useful adjusted summary while the classical standard error and p-value may be poorly calibrated. The test’s inferential procedure must be appropriate to the dependence and variance structure of the data; a small ordinary-model p-value is not a substitute for checking that issue.

Worked Example: Testing an Adjusted Pathway Difference

Hypothetical data. Suppose the analysis includes \(n=100\) genes. The outcome is a gene-level score on a scale where larger values indicate stronger association evidence. The fitted model includes an intercept, pathway membership, standardized gene length, and standardized variant count. The pathway coefficient estimate is \(0.42\) score units, with standard error \(0.15\) score units.

There are four fitted coefficients, so the residual degrees of freedom are \(df=100-4=96\). Under the classical model assumptions, the test statistic for \(H_0:\beta_S=0\) is

$$ t=\frac{\hat{\beta}_S}{SE(\hat{\beta}_S)} =\frac{0.42}{0.15} =2.80. $$

The two-sided p-value is \(2P(T_{96}\ge 2.80)\), which is approximately \(0.006179\). If the classical assumptions and the prespecified testing procedure are appropriate, this is evidence against a zero adjusted difference. It is not the probability that the null hypothesis is true.

Using \(t_{0.975,96}\approx1.985\), a 95% confidence interval is

$$ 0.42\pm 1.985(0.15) =0.42\pm0.29775 \approx[0.1223,\ 0.7178]. $$

In this hypothetical analysis, the interval estimates a positive adjusted difference of about \(0.12\) to \(0.72\) score units for pathway members relative to nonmembers, conditional on gene length and variant count. The interval describes uncertainty under the fitted model; it does not give a 95% probability that the fixed coefficient lies in this particular interval.

Why the Coefficient Is an Adjusted Comparison

The adjustment can be understood by separating the pathway indicator from the covariates. Write \(Z\) for the design matrix containing the intercept and all covariates, and write \(S\) and \(Y\) for the vectors of pathway indicators and scores. Let \(P_Z=Z(Z^\mathsf{T}Z)^{-1}Z^\mathsf{T}\) be the projection onto the columns of \(Z\), assuming those columns are linearly independent. The residualized vectors are \(\tilde S=(I-P_Z)S\) and \(\tilde Y=(I-P_Z)Y\).

The least-squares normal equations require the fitted residual to be orthogonal to every predictor. After accounting for \(Z\), the remaining equation for the pathway coefficient is \(\tilde S^\mathsf{T}(\tilde Y-\hat{\beta}_S\tilde S)=0\). Solving it gives the coefficient formula below. This shows that the estimate compares the part of the scores associated with pathway membership after both membership and scores have been adjusted for the same covariates.

$$ \hat{\beta}_S = \frac{\tilde S^\mathsf{T}\tilde Y} {\tilde S^\mathsf{T}\tilde S}, \qquad \tilde S^\mathsf{T}\tilde S>0. $$

If \(\tilde S^\mathsf{T}\tilde S=0\), pathway membership is completely determined by the covariates in the fitted sample, so the pathway coefficient cannot be estimated separately. More generally, if membership is nearly determined by them, the estimate can be unstable and its standard error large.

Worked Example: Calculating an Adjusted Coefficient

Hypothetical data. Consider eight genes, four pathway members and four nonmembers. Let \(x\) be standardized gene length. The score values are shown below. Within each membership group, the length values are \(-1,-1,1,1\).

MembershipLength \(x\)Scores \(Y\)
Nonmember-1, -11.0, 1.2
Nonmember1, 11.4, 1.6
Member-1, -11.5, 1.7
Member1, 11.9, 2.1

Fit \(Y_g=\beta_0+\beta_S S_g+\gamma x_g+\varepsilon_g\). Because each group has the same length values, membership is orthogonal to centered length: adjusting membership for the intercept and \(x\) leaves \(\tilde S_g=S_g-0.5\). The mean score for members is \((1.5+1.7+1.9+2.1)/4=1.80\), and for nonmembers it is \((1.0+1.2+1.4+1.6)/4=1.30\). The adjusted coefficient is therefore

$$ \hat{\beta}_S = \frac{4(0.5)(1.80-1.30)} {4(0.5)^2+4(0.5)^2} = \frac{1.00}{2.00} =0.50. $$

The fitted length slope is \(0.20\) score units per one-standard-deviation increase in length: the mean score rises from \(1.10\) to \(1.50\) among nonmembers as \(x\) rises from \(-1\) to \(1\). The estimated pathway contrast remains \(0.50\) score units because length is balanced between the two groups. This example demonstrates the coefficient calculation; with only eight constructed observations, it is not a basis for a reliable pathway conclusion.

Applying the Model Carefully

A defensible analysis begins by fixing the gene universe and the outcome construction. Genes excluded because they lack a score should be handled transparently, since exclusion can change which genes the pathway is compared with. Define membership consistently, including how genes mapped to multiple annotations or pathway definitions are treated. Then specify the covariates and their units or transformations before interpreting the coefficient.

For a score on a logarithmic scale, \(\beta_S\) is measured in logarithmic score units, not in odds or disease-risk units. If the score is \(-\log_{10}(p)\), higher values indicate smaller gene-level p-values, but the coefficient is not a change in a gene’s causal effect. The genes are the analysis units, and the response is the chosen summary score.

Worked Example: Interpreting a Fitted Model

Hypothetical fitted equation. Suppose an analysis uses a score measured in score units and fits

$$ \widehat{Y}_g = 1.10+0.35S_g+0.20Z_{\text{length},g} +0.08Z_{\text{count},g}, $$

where \(Z_{\text{length}}\) and \(Z_{\text{count}}\) are standardized gene length and variant count. Consider otherwise comparable genes with \(Z_{\text{length}}=1\) and \(Z_{\text{count}}=0\). For a nonmember, the predicted score is \(1.10+0+0.20(1)+0.08(0)=1.30\). For a member, it is \(1.10+0.35+0.20(1)+0.08(0)=1.65\). The difference is \(1.65-1.30=0.35\) score units, exactly the pathway coefficient.

This is a conditional model comparison at the stated covariate values. It is not evidence that changing a gene’s pathway membership would change its association score, nor does the equation alone provide uncertainty or a p-value; those require the coefficient’s standard error and a suitable inferential model.

Assumptions, Limits, and a Common Mistake

Adjustment is not a guarantee of unbiased inference. The linear mean specification may be inadequate if score relationships with gene size or variant count are strongly nonlinear. A covariate measured poorly may not control the feature it was intended to represent. Covariate adjustment can also change the scientific comparison: \(\beta_S\) is conditional on the selected covariates, not a universal measure of pathway enrichment.

A further distinction matters: the coefficient test asks whether pathway membership is associated with a difference in scores relative to genes outside the pathway after adjustment. It does not test the separate question of whether genes within the set show any association compared with a null of no association. These are different estimands, and the interpretation must match the question actually tested.

  • Mistake: “The pathway coefficient is positive, so the pathway causes the trait.” Correction: The coefficient describes an adjusted association between pathway membership and the chosen gene score. It does not establish causation or identify a biological mechanism.
  • Mistake: “Including gene length and variant count makes the ordinary p-value valid.” Correction: Covariate adjustment addresses the modeled conditional comparison; it does not ensure independent, constant-variance errors or correct calibration when gene scores are dependent.
  • Mistake: “A significant coefficient means all pathway genes are associated.” Correction: It is a set-level contrast in average modeled scores, not a conclusion about every member gene.

Regression is useful because it makes the target comparison and its covariate conditions explicit. Its value depends on the score definition, gene universe, covariate choices, and inferential assumptions being made clear alongside the coefficient. Report the estimate in score units, its uncertainty, the model covariates, and the limitations of the reference distribution used for the test.

Check Your Understanding

Answer each question using the regression estimand and assumptions described above.

  1. In a model of gene scores on pathway membership, gene length, and variant count, what does the pathway coefficient estimate?
  2. Why can a pathway coefficient be unestimable when membership is completely determined by the covariates?
  3. A positive coefficient is estimated for a score where larger values indicate stronger evidence. What does that say, and what does it not say?
  4. What assumptions are needed for the classical t reference distribution to be exact?
  5. Why does adjusting for gene size and variant count not, by itself, guarantee a valid standard error?
Clinical Biostats