Worked Integration of Expression, Protein, and Metabolite Data

From Separate Layers to One Integrated Analysis

A genetic variant can be associated with expression of a gene, abundance of its protein product, and a downstream metabolite. Integrating these measurements means more than placing three association results side by side. The analysis must keep track of which participants contribute to each result, how uncertainty is shared across layers, and which claims are supported by association versus by a mediation model.

The worked analyses below use a hypothetical cohort and hypothetical estimates. They illustrate a coherent workflow: first define the measurements and samples, then test whether the variant is associated with the molecular outcomes jointly, and then estimate a prespecified pathway model. A joint association test asks whether the outcome vector is associated with genotype; it does not by itself establish that expression causes protein abundance or that protein abundance causes a metabolite change.

Definition: A multi-omic genetic integration analysis combines genetic data with measurements from two or more molecular layers, such as gene expression, protein abundance, and metabolites. A joint association test tests a null hypothesis about the genetic effects across a vector of outcomes. A mediation analysis estimates effects along specified pathways through intermediate variables, under additional modeling and identification assumptions.

Set Up the Hypothetical Dataset

Suppose a study measures one variant and three molecular outcomes in 200 unrelated participants. Genotype \(G\) is coded as the number of copies of a specified allele, from 0 to 2. Expression \(E\) is normalized log2 expression, protein abundance \(P\) is normalized log2 abundance, and metabolite \(M\) is the natural log of concentration. Let \(C\) denote a fixed set of adjustment variables, including age, sex, and prespecified ancestry covariates. The units matter: an effect on log2 expression is not an effect in copies of RNA, and an effect on log metabolite concentration is not an effect in the original concentration scale.

For a layer-specific association analysis, a basic linear model for each outcome \(Y_k\), where \(k\) indexes expression, protein, or metabolite, is

Y_{ik}=\beta_{0k}+\beta_{Gk}G_i+\boldsymbol{\gamma}_k^\mathsf{T}C_i+\varepsilon_{ik}.

Here \(\beta_{Gk}\) is the adjusted mean difference in outcome \(k\) per additional allele, conditional on the included covariates. The usual linear-model assumptions include an appropriate conditional mean specification, independent participants, and a valid estimate of the coefficient covariance. Normal residuals give exact small-sample inference under the model; with a sufficiently large sample, Wald inference can instead rely on an appropriate large-sample covariance estimate.

Worked Example: Mapping Measurements to Participants

In the hypothetical cohort, all 200 participants have genotype data. Expression is measured for 192, protein for 178, and metabolite for 165. A total of 160 participants have all three molecular outcomes observed. A data table has one row per participant, with genotype, covariates, and separate columns for \(E\), \(P\), and \(M\).

The 160 complete records can be used for an analysis requiring all three outcomes observed for the same people. The larger layer-specific samples can be used for separate association models, provided the sample differences are reported. For example, the expression association has 192 observations, not 200, because expression is unavailable for eight participants. Treating a missing outcome as zero would change the measurement definition and would not be a valid way to retain those participants.

Before interpreting estimates across layers, check that the outcomes have the stated scales, genotype is coded consistently, and covariates are harmonized. A coefficient of 0.20 for expression and 0.20 for log metabolite concentration have different units and should not be compared as though they were the same effect size.

Test Association Across the Three Layers

On the 160 complete participants, suppose the estimated genetic effects from three covariate-adjusted marginal models are \(0.28\) log2 expression units per allele, \(0.18\) log2 protein units per allele, and \(0.12\) log metabolite units per allele. Their standard errors are \(0.08\), \(0.07\), and \(0.06\), respectively. The estimates are correlated because they are measured in the same participants and the outcomes themselves may co-vary.

A covariance-aware Wald test evaluates the joint null hypothesis that all three genetic effects are zero. If \(\hat{\boldsymbol{\beta}}\) is the vector of estimates and \(\widehat{\boldsymbol{\Sigma}}\) is its estimated covariance matrix, the statistic is

W=\hat{\boldsymbol{\beta}}^\mathsf{T} \widehat{\boldsymbol{\Sigma}}^{-1} \hat{\boldsymbol{\beta}}.

Under the joint null and the model’s regularity conditions, \(W\) is approximately chi-square distributed with three degrees of freedom. The covariance matrix must describe the estimates together; using only three separate standard errors discards information about their dependence.

Worked Example: A Covariance-Aware Joint Association Test

For the estimates above, suppose the estimated correlation between each pair of coefficient estimates is 0.25. This is a hypothetical covariance structure. The standard errors are \(0.08\), \(0.07\), and \(0.06\), so the covariance matrix is the product of the diagonal matrix of standard errors, the correlation matrix, and the diagonal matrix of standard errors. The standardized estimates are

\mathbf{z}= \left( \frac{0.28}{0.08}, \frac{0.18}{0.07}, \frac{0.12}{0.06} \right)^\mathsf{T} =(3.50,\ 2.571,\ 2.00)^\mathsf{T}.

With an equicorrelation matrix having off-diagonal entries 0.25, the Wald statistic can be calculated from the standardized estimates as

W=\mathbf{z}^\mathsf{T}R^{-1}\mathbf{z} =\frac{1}{0.75} \left[ \sum_{j=1}^{3}z_j^2 -\frac{0.25}{1+2(0.25)} \left(\sum_{j=1}^{3}z_j\right)^2 \right] \approx 16.01.

Under \(H_0:\boldsymbol{\beta}=\mathbf{0}\), this is compared with a chi-square distribution with \(df=3\). The upper-tail p-value is approximately 0.0011. Thus, these data provide evidence against the joint null of no genetic association with any of the three outcomes. The result does not show that every outcome has a nonzero effect, nor does it identify a causal order among the outcomes.

For context, the separate 95% Wald intervals are \(0.28\pm1.96(0.08)=(0.123,0.437)\) log2 expression units per allele, \(0.18\pm1.96(0.07)=(0.043,0.317)\) log2 protein units per allele, and \(0.12\pm1.96(0.06)=(0.002,0.238)\) log metabolite units per allele. These intervals summarize marginal estimates; they are not substitutes for the joint test.

Estimate a Prespecified Molecular Pathway

A mediation model answers a different question. Consider the ordered structure \(G\rightarrow E\rightarrow P\rightarrow M\), while allowing direct genetic effects and shorter pathways. One linear specification is

\begin{aligned} E_i &= \alpha_E+aG_i+\boldsymbol{\gamma}_E^\mathsf{T}C_i+\varepsilon_{Ei},\\ P_i &= \alpha_P+cG_i+dE_i+\boldsymbol{\gamma}_P^\mathsf{T}C_i+\varepsilon_{Pi},\\ M_i &= \alpha_M+tG_i+eE_i+bP_i+\boldsymbol{\gamma}_M^\mathsf{T}C_i+\varepsilon_{Mi}. \end{aligned}

In this model, \(a\) is the genotype-to-expression path, \(d\) is the expression-to-protein coefficient conditional on genotype and covariates, and \(b\) is the protein-to-metabolite coefficient conditional on genotype, expression, and covariates. The serial indirect effect through expression and then protein is \(adb\). Other indirect paths are \(ae\) through expression directly to the metabolite, and \(cb\) through the genotype-to-protein path and then protein to metabolite. The direct effect is \(t\).

These products have a causal interpretation only under suitable identification conditions. In particular, the model needs correct temporal and structural specification, adequate control of relevant confounding, and no unaddressed mediator–outcome confounding induced by genotype. The genotype’s presence alone does not guarantee these conditions. Without them, the products are model-based pathway estimates, not established causal effects.

Worked Example: Decomposing the Genetic Effect on a Metabolite

Suppose the hypothetical fitted pathway model gives \(a=0.30\) log2 expression units per allele, \(d=0.40\) log2 protein units per log2 expression unit, \(c=0.10\) log2 protein units per allele, \(e=0.05\) log metabolite units per log2 expression unit, \(b=0.25\) log metabolite units per log2 protein unit, and \(t=0.02\) log metabolite units per allele.

The serial indirect estimate is the product along \(G\rightarrow E\rightarrow P\rightarrow M\):

\widehat{IE}_{G\to E\to P\to M} =\hat{a}\hat{d}\hat{b} =(0.30)(0.40)(0.25) =0.030

This is an estimated 0.030 log metabolite units per allele through the specified serial path. The two shorter indirect paths are \(ae=(0.30)(0.05)=0.015\) and \(cb=(0.10)(0.25)=0.025\), each in log metabolite units per allele. Substituting the linear models into one another shows why the total genetic effect is the sum of these paths and the direct path:

\begin{aligned} P &: \quad cG+d(aG+\cdots)+\cdots=(c+da)G+\cdots,\\ M &: \quad tG+e(aG+\cdots)+b\{(c+da)G+\cdots\}+\cdots,\\ \beta_{G,\mathrm{total}}&=t+ae+cb+adb\\ &=0.020+0.015+0.025+0.030=0.090. \end{aligned}

The total estimate is therefore 0.090 log metabolite units per allele under this model. It is not a percentage mediated: the pathways are expressed on a log scale, and a ratio of an indirect estimate to the total can be unstable when the total is near zero or when pathways have opposing signs.

Missing Measurements and Uncertainty

A complete-case analysis is transparent, but it discards participants who lack any one layer. An observed-data likelihood can use incomplete records by contributing the probability density for the components that are observed for each participant. For example, under a specified joint model for \((E,P,M)\) conditional on genotype and covariates, a participant with \(E\) and \(M\) observed contributes their bivariate marginal density, rather than a fabricated value for \(P\).

This approach requires a correctly specified model and assumptions that make the missingness mechanism ignorable for the likelihood analysis. A common condition is missing at random: after conditioning on observed measurements and covariates in the model, the probability a value is missing does not depend on that unobserved value. If missingness depends on an unobserved measurement even after conditioning, the mechanism is nonrandom and requires further modeling or sensitivity analysis. The fact that a likelihood can accommodate incomplete records does not itself establish that missingness is ignorable.

Worked Example: Using Incomplete Omics Records

In the hypothetical cohort, eight expression values, 22 protein values, and 35 metabolite values are missing; 160 participants have all three outcomes observed. Suppose an observed-data likelihood fit uses each participant’s available outcomes under a joint multivariate model conditional on genotype and covariates. This fit can use records such as a participant with expression and metabolite observed but protein missing, without assigning that participant an invented protein value.

Suppose the model estimates the serial indirect effect as 0.030 log metabolite units per allele, with a participant-level bootstrap 95% percentile interval of 0.008 to 0.058. In each of 1,000 hypothetical bootstrap resamples, participants are resampled as whole units, the model is refitted, and the product \(\hat{a}\hat{d}\hat{b}\) is recalculated. The interval is formed from the 2.5th and 97.5th percentiles of those 1,000 product estimates. It describes sampling uncertainty under the fitted model and resampling scheme; it does not account automatically for model misspecification or violation of the missing-at-random assumption.

The interval excludes zero, so the fitted analysis provides evidence that the serial product differs from zero under its assumptions. A 95% confidence interval is not a 95% probability statement about the fixed parameter. Also, because participants rather than individual measurements are resampled, the resampling preserves the within-participant structure among the observed layers.

Practical Checks Before Interpretation

An integrated analysis is most useful when its test, pathway, and uncertainty summaries answer clearly distinguished questions. The joint Wald test concerns the vector of genetic associations. The mediation model concerns coefficients along a prespecified structure. The missing-data model concerns how incomplete records contribute information. Neither a small joint-test p-value nor a nonzero pathway product, by itself, verifies the assumptions needed for a causal interpretation.

Define the analysis population.
Report the sample size for each layer and for any analysis requiring participants with multiple outcomes observed.

Specify scales and models.
State genotype coding, outcome units, adjustment variables, and the regression used for each association or mediation path.

Match the test to the question.
Use the covariance of effect estimates for a joint test; do not infer a pathway order from joint association alone.

Summarize uncertainty and assumptions.
Give intervals with units, state the missingness assumption, and separate statistical evidence from causal interpretation.

Common mistake. It is tempting to say, “The variant is jointly associated with all three layers, so the expression-to-protein-to-metabolite chain is confirmed.” That conclusion does not follow. The joint test rejects only the null that all three genetic effects are zero; it neither tests every link in the chain nor establishes mediation. The correction is to report the joint association result as such, estimate the specified path products separately, and qualify causal claims by the mediation and missing-data assumptions.

Key takeaway: Integration requires aligned samples, compatible model definitions, and covariance-aware uncertainty. A joint association test can establish evidence across an outcome vector, while a pathway model estimates specified indirect effects; neither should be interpreted beyond its assumptions.

Check Your Understanding

Use the hypothetical models and estimates above to answer these questions.

What null hypothesis is tested by the three-outcome Wald statistic, and what does rejecting it not establish?
Why is the covariance among the three genetic effect estimates relevant to the joint test?
Using the stated coefficients, calculate the expression-to-metabolite indirect path \(ae\).
What missingness assumption is commonly needed for an observed-data likelihood to ignore the missingness mechanism?
Why should a bootstrap resample participants rather than independently resample expression, protein, and metabolite measurements?