Microbiome and Metabolome Integration in Inflammatory Disease

Study Question and Analysis Plan

In an inflammatory-disease study, stool samples and blood samples may be collected from the same patients. The stool assay records microbial sequencing counts, while the blood assay measures metabolite concentrations. The scientific question is whether a microbial abundance pattern is associated with a metabolite profile after accounting for measured clinical characteristics such as age, disease severity, medication use, and assay batch.

The patient is the unit of analysis: each row in the analysis corresponds to one patient with paired measurements. The analysis estimates conditional associations, not causal effects. For example, an association that remains after adjusting for age and severity does not establish that a microbe changes a metabolite; unmeasured factors, reverse direction, or selection into the study could still explain the association.

A useful workflow has three linked decisions: represent microbiome data on a scale suitable for relative abundances, transform metabolite measurements if needed, and model their relationship while including prespecified clinical covariates. In a hypothetical study of 60 patients with inflammatory bowel disease, the target might be the association between a log-ratio microbial feature and each of 20 measured metabolites, adjusted for age, severity, medication, and batch.

Representing the Two Omics Layers

Definition: A microbiome composition is a vector of nonnegative component abundances whose relative values are observed, often as sequencing counts. A centered log-ratio (CLR) coordinate for component $j$ in a sample with positive components $x_1,\ldots,x_D$ is $$ \operatorname{clr}(x_j)=\log\left(\frac{x_j}{g(x)}\right), \qquad g(x)=\left(\prod_{k=1}^{D}x_k\right)^{1/D}. $$ It expresses the component relative to the geometric mean of the components in that sample.

Sequencing counts are not direct measurements of absolute microbial load. Because their total is constrained by sequencing depth, an increase in one component's relative share can make other shares decrease even if their absolute abundances do not change. A log-ratio feature is therefore interpreted as a relative abundance pattern, not as an absolute change in a taxon's cell count. CLR coordinates sum to zero within each sample, so they are not independent features; including every CLR coordinate alongside an intercept can create exact collinearity.

Zeros need explicit treatment before taking logarithms. A small pseudocount is one simple option for a teaching example, but its choice can affect results, especially for rare taxa or many zeros. In applied work, report the zero-handling rule and assess whether conclusions depend on it. Filtering rare features may also be appropriate, provided the rule is set without selecting features based on the association results.

Metabolite concentrations are often right-skewed, so a logarithm can make their distribution more suitable for a linear model. The transformation should match the measurement scale and the intended interpretation. If a log-base-2 concentration is the outcome, a coefficient of 1 corresponds to a two-fold concentration difference, conditional on the covariates.

Worked Example: Converting Counts to CLR Coordinates

Suppose one hypothetical stool sample has counts of 40, 20, and 10 for three retained microbial components. Add a pseudocount of 1 to each count, giving $41,21,11$. The geometric mean is $$ g=(41\cdot21\cdot11)^{1/3}=9471^{1/3}\approx21.16. $$

The CLR coordinates are $$ \log(41/21.16)\approx0.6616,\quad \log(21/21.16)\approx-0.0075,\quad \log(11/21.16)\approx-0.6541. $$ Their sum is approximately zero, as expected apart from rounding. The first component is relatively high compared with the sample's geometric mean; this calculation does not show that its absolute abundance is high. The pseudocount is part of the example's specified preprocessing rule, not a universal constant.

Covariate-Adjusted Association Model

Definition: For patient $i$, let $Y_i$ be a transformed metabolite measurement, $X_i$ a prespecified microbial log-ratio feature, and $C_{i1},\ldots,C_{iq}$ measured clinical or technical covariates. The additive linear model is $$ Y_i=\beta_0+\beta_X X_i+\sum_{r=1}^{q}\gamma_r C_{ir}+\varepsilon_i. $$ Here $\beta_X$ is the conditional mean difference in $Y$ associated with a one-unit increase in $X$, holding the included covariates fixed.

The standard small-sample regression test and confidence interval assume independent patient-level errors, a correctly specified conditional mean, no exact collinearity among model terms, and constant error variance. Exact $t$-based inference additionally assumes normally distributed errors; with larger samples, inference may rely on an approximation, while heteroscedasticity-robust standard errors can address unequal variance under their own large-sample conditions. Repeated samples from one patient are not independent observations and require a model that accounts for that dependence.

Covariates should be chosen using subject-matter knowledge and the study design. Include plausible confounders measured before or alongside the exposures, and consider technical variables such as assay batch when they may influence measurements. Adjustment is not automatically beneficial: adjusting for a variable on a causal pathway can change the question being answered, and conditioning on a collider can introduce bias. A model also cannot repair perfect confounding, such as a medication group occurring only in one assay batch.

Worked Example: Estimating an Association Adjusted for Severity

Consider an illustrative subset of eight hypothetical patients. $X$ is a microbial log-ratio feature, $Y$ is a metabolite concentration on a log-base-2 scale, and severity is coded as group 0 or group 1. The model is $Y_i=\beta_0+\beta_X X_i+\gamma\,\text{severity}_i+\varepsilon_i$.

Severity group	$X$ values	$Y$ values
0	1, 2, 2, 3	8, 9, 11, 12
1	2, 3, 3, 4	12, 13, 15, 16

Within each severity group, the mean $X$ values are 2 and 3, and the mean $Y$ values are 10 and 14. After subtracting each group's means, the $X$ deviations in either group are $-1,0,0,1$, and the $Y$ deviations are $-2,-1,1,2$. The adjusted slope is the within-group cross-product divided by the within-group sum of squares: $$ \hat{\beta}_X= \frac{\sum_i (X_i-\bar X_{g(i)})(Y_i-\bar Y_{g(i)})} {\sum_i (X_i-\bar X_{g(i)})^2} =\frac{2(4)}{2(2)}=2. $$ The group means are removed here because the model includes a severity indicator. The slope means that, within the same severity group, a one-unit higher microbial log-ratio is associated with a 2-unit higher log-base-2 metabolite concentration. On the original concentration scale, that corresponds to a four-fold ratio, since $2^2=4$; it is an association, not evidence of a causal four-fold increase.

The fitted residuals within each group are $0,-1,1,0$, so the residual sum of squares is $2+2=4$. There are $8-3=5$ residual degrees of freedom (intercept, microbial feature, and severity indicator), giving $MSE=4/5=0.8$. The within-group $X$ sum of squares is 4, so $$ SE(\hat{\beta}_X)=\sqrt{0.8/4}=0.4472,\qquad t=\frac{2}{0.4472}=4.472,\quad df=5. $$ Under the model assumptions, the two-sided $p$-value is approximately 0.0066. A 95% $t$-interval is $2\pm2.571(0.4472)=[0.850,3.150]$ log-base-2 units per log-ratio unit. This interval describes uncertainty in the population slope under the model; it is not a 95% probability statement about this fixed parameter.

Why the Adjustment Works

The calculation above is a special case of the Frisch–Waugh–Lovell result: in a linear model, the coefficient of $X$ after adjusting for covariates equals the slope obtained by removing the covariate-predicted parts of both $X$ and $Y$, then regressing the remaining part of $Y$ on the remaining part of $X$. For a group indicator, removing the group-predicted part means subtracting each group's mean.

To see why, write $X_i=\bar X_{g(i)}+X_i^*$ and $Y_i=\bar Y_{g(i)}+Y_i^*$, where starred values are deviations from their group means. Within each group, the deviations sum to zero, so group-specific intercept terms contribute nothing to the cross-product of the deviations. Minimizing the sum of squared residuals over the slope therefore gives $$ \hat{\beta}_X=\frac{\sum_i X_i^*Y_i^*}{\sum_i (X_i^*)^2}, $$ provided the denominator is positive. If it is zero, $X$ has no within-group variation, and its association cannot be separated from the group effect in this model. For continuous covariates, the same logic applies after residualizing on all covariates, rather than simply subtracting group means.

Testing Many Metabolites and Checking Diagnostics

If the same microbial feature is tested against many metabolites, each regression answers a separate null hypothesis $H_0:\beta_X=0$. Testing many hypotheses increases the chance of false discoveries. The Benjamini–Hochberg procedure controls the expected false discovery rate under its standard dependence conditions. Sort the $m$ raw $p$-values from smallest to largest, and compare the $r$-th value with $r\alpha/m$. Adjusted $p$-values summarize the corresponding step-up procedure and are more informative than presenting only a thresholded list.

Worked Example: Applying Benjamini–Hochberg to Five Metabolites

Suppose five hypothetical metabolite tests yield sorted $p$-values $0.004,0.012,0.030,0.080,0.200$. At false discovery rate level $\alpha=0.05$, the ordered thresholds are $0.010,0.020,0.030,0.040,0.050$. The first three values meet their thresholds, including equality for the third; the last two do not. The step-up rule therefore identifies the first three tests as discoveries at this level.

The corresponding Benjamini–Hochberg adjusted values are obtained from $mp_{(r)}/r$, then made nondecreasing from the largest rank back to the smallest. Here they are $0.020,0.030,0.050,0.100,0.200$. A value of 0.050 is not the probability that a particular finding is false; the procedure controls the expected proportion of false discoveries among the rejected hypotheses under its assumptions.

Diagnostics should be checked for each fitted association, especially for findings that drive the scientific conclusion. Plot residuals against fitted values to look for curvature or changing variance; inspect residual distributions for severe departures that may undermine small-sample $t$-inference; and examine leverage and influence to determine whether one patient dominates the slope. Check whether the microbial feature varies enough within covariate groups and whether it is strongly collinear with severity, medication, or batch. Inspect missingness and sample pairing as well: a row mismatch between assays can create a technically valid but scientifically meaningless regression.

If a linear relationship is implausible, consider a scientifically justified transformation or nonlinear term and report the change in interpretation. If variance is unequal, robust standard errors may be appropriate. If the number of taxa-metabolite tests is large, control multiplicity over the family of tests that supports the stated claims. These steps do not eliminate the need to report preprocessing, covariate choices, effect estimates, uncertainty, and the full testing scope.

Alternatives and Limits of the Case-Study Approach

An unadjusted correlation is simpler but does not account for clinical or technical covariates; it can reflect differences in severity or batch rather than a microbial–metabolite relationship within comparable patients. Adjusted regression makes the covariates and the target coefficient explicit, and gives an effect estimate with units. A joint multivariate test can instead ask whether a microbial feature is associated with a set of metabolites considered together, potentially accounting for their correlation. That omnibus question is different from identifying which individual metabolites are associated, and it does not replace appropriate follow-up and multiplicity control.

Pairwise regression is also limited when there are many correlated taxa and metabolites relative to the number of patients. Estimates may be unstable, and a result for one feature can depend on which other features are included. CLR transformation handles the relative scale but does not make components independent or recover absolute abundance. Clinical covariates are measured imperfectly, unmeasured confounding may remain, and a cross-sectional association cannot establish temporal order. These limitations should shape the strength of the conclusion rather than be hidden by a small $p$-value.

Common Mistake

Common pitfall: “The taxon's CLR value increased, so its absolute abundance increased and caused the metabolite concentration to rise.” This conclusion makes two unsupported leaps. A CLR value is relative to the sample's geometric mean, and an adjusted observational association is not a causal effect. State the result as a conditional association between a specified log-ratio feature and a transformed metabolite, report the units and uncertainty, and describe causal interpretation as unresolved.

Check Your Understanding

Answer the questions, then compare your reasoning with the solutions.

Why is a CLR coefficient not directly an effect per unit increase in absolute microbial cell count?
In a model adjusted for severity, what does the microbial-feature slope represent?
What problem arises if all samples in one medication group were measured in batch A and all others in batch B?
For sorted $p$-values $0.01,0.04,0.20$, which values meet the Benjamini–Hochberg step-up threshold at $\alpha=0.05$?

Solutions. (1) CLR is a log ratio to the sample geometric mean and depends on the composition, not an absolute count scale. (2) It is the conditional mean difference in the modeled metabolite outcome for a one-unit increase in the microbial feature, holding severity and the other included covariates fixed. (3) Medication and batch are perfectly confounded, so their separate effects cannot be identified from these data. (4) With $m=3$, the thresholds are approximately $0.0167,0.0333,0.0500$. The first value meets its threshold, while the second and third do not; the step-up rule therefore rejects only the first hypothesis.