Research has traditionally focused on the health impacts of one chemical at a time. In reality, individuals are simultaneously exposed to potentially thousands of environmental chemicals, each of which can add to the cumulative health burden. Based on chemical inventories around the world, over 355,000 chemicals or chemical mixtures have been registered for production and use, including 69,000 in the previous decade alone (2010-2019).¹ For roughly 15% of those registered substances, descriptive chemical names are not provided publicly due to confidential business information, and another ∼15% are only ambiguously described.¹ Considering plastic production alone, more than 16,000 chemicals are used or present in products, of which over 4,200 are potentially hazardous and over 10,000 lack hazard information.^2,3 As a result, humans are exposed to complex mixtures of both known and unidentified chemicals, a large fraction of which have not been comprehensively studied for toxicological safety or health.

A key challenge in identifying new chemical exposures has been the phenomenon of regrettable substitution. Even when one toxic chemical is eventually phased out of production, another chemical—often with a similar structure and same chemical class—may replace it and only later be revealed to have toxicity concerns as well. Such chemical ‘whack-a-mole’ has been observed countless times for phenols in plastic,⁴ phthalates in plastic,⁵ flame retardants in furniture and electronics,⁶ and per- and polyfluoroalkyl substances (PFAS) in consumer products⁷—of which there exist at least 14,000 types of PFAS.⁸ Consequently, traditional targeted research methods that measure only a limited number of pre-selected known chemicals cannot match the pace at which new chemicals enter commerce. For example, large human biomonitoring programs have usually analyzed up to about 300 targeted chemicals, due to limitations from available analytical standards, cost, and time.^9,10

The exposome was first conceptualized in 2005 in response to the evident gulf in the scale at which genomic versus environmental risk factors can be characterized in biospecimens and the need to advance methodologies for measuring biomarkers of environmental exposure at a similar omics level.¹¹ The most recent definition of exposomics is the study of “the comprehensive and cumulative effects of physical, chemical, biological, and psychosocial influences that impact biological systems by integrating data from a variety of interdisciplinary methodologies and streams to enable discovery-based analysis of environmental influences on health.”¹² As a key aspect of the exposome, untargeted chemical profiling using high-resolution mass spectrometry (HRMS) now supports the simultaneous measurement and characterization of over 100,000 chemical signals in human biospecimens, including both the internal exposome (ie, exogenous environmental chemicals and their biotransformation products) and the metabolome (ie, endogenous metabolites).^13-15

The untargeted approach encompasses even unknown chemicals that we cannot yet identify but can still monitor in samples and elucidate structural information for.¹⁶ Each untargeted chemical feature is characterized in the mass spectrometry output by its accurate mass-to-charge ratio (m/z) and retention time.¹⁷ Although the data generated are relative abundances (ie, ion intensities) rather than absolute targeted concentrations,¹⁸ HRMS offers a powerful hypothesis-free discovery approach for detecting previously unknown environmental risk factors of disease as well as the metabolic responses that underlie this risk, without requiring analytical standards for the chemicals.^19,20 With the high dimensionality of untargeted chemical data produced, statistical analysis remains a major challenge in the field of exposomics. The objective of our commentary is to share a statistical workflow that we have developed for evaluating the cumulative mixture effects of untargeted chemicals on disease, prioritizing novel chemicals of concern, and exploring mechanisms of action. This workflow has been applied in our recent studies of exposome drivers of fertility outcomes²¹ and lymphoma risk.

Although statistical mixture methods embrace the complexity of chemical exposures, most approaches cannot handle untargeted data at the omics-scale. The major barrier with untargeted omics data lies in the fact that the number of parameters is far greater than the sample size (p ≫ n), and thus would overfit the model to an impossible degree. In addition, untargeted data can be even more highly correlated because of the presence of some redundant signals, such as adducts, isotopes, or fragments of the same chemical, and multiple biotransformation products arising from the same exogenous chemical.³⁶ For this reason, not all mixture methods can be readily scaled up to the high dimensions of the exposome.

Choice of mixture method largely depends on the research question.²² In our case, we aim to use mixture methods to evaluate cumulative mixture effects from untargeted chemical exposures and identify the ‘bad actor’ chemicals driving the mixture effect the most. There are several common statistical approaches that would not work best for our research questions.

For example, some studies calculate the sum, molar sum, or potency-weighted sum of specific chemical classes to model the total effect on health using fewer variables.³⁷ This approach has been commonly applied to phthalates or phthalate subgroups, such as the summed urinary metabolites of di(2-ethylhexyl) phthalate (DEHP) or high versus low molecular weight phthalates.²⁹ However, chemical sums lose data resolution and can hinder interpretations about individual chemicals for decision-making. In addition, the summation could mask hazards if chemicals in the same class have different molecular weights, abundances, toxicological relevance, interactions, concentrations at which adverse effects occur, or even opposing directions of effects. For instance, the health effect of a low-concentration chemical may not be observed when summed with a high-concentration chemical that does not affect the outcome.³²

Bayesian kernel machine regression (BKMR) is a popular mixture method³⁸ that can model the joint effects of exposures using a flexible kernel function that allows for non-linear, smooth exposure-outcome relationships and interactive effects between chemicals (with some sacrifices to statistical power).^39,40 Although useful for targeted chemical mixtures, especially research questions related to non-linearity and interactions between chemicals, BKMR is currently not practical for higher dimensional chemical data (p ≫ n) because of the large sample size and computational intensity needed for the non-parametric kernel function, which is highly flexible but less statistically powerful.⁴¹ In addition, the visual interpretation of curvilinearity and interaction within high-dimensional mixtures in BKMR would be intractable. For applications of mixture methods to lower-dimensional mixtures, a recent publication offers a helpful workflow for statistical decisions related to the distribution and type of data, variable transformations, missing data, statistical assumptions, specific research questions, and other study design considerations.⁴²

Quantile g-computation (QGcomp) is another common mixture method similar to WQS that builds a weighted summary index based on quantiles of exposures, with a few slight advantages or disadvantages depending on the desired research question and sample size. First, QGcomp relaxes the directional homogeneity assumption by incorporating opposing effects within the same index,⁴³ whereas WQS now sequentially defines separate positive and negative indices to evaluate double associations of the mixture in both directions.^44,45 Second, QGcomp allows for more non-additivity and non-linearity of effects if specified as model terms in advance. Third, QGcomp does not split samples into repeated training/validation sets, which has tradeoffs between statistical power and computational speed versus generalizability and weight stability.^43,46 However, for the purpose of untargeted mixture methods in this paper, QGcomp has not yet been designed for applications to high-dimensional data.⁴⁷

Other statistical approaches seek to systematically reduce the dimensions of data before modeling. For example, principal component analysis (PCA) is an unsupervised approach that converts the chemical data into a smaller set of uncorrelated linear predictors (not taking into account the outcome). PCA may be helpful for assessing patterns in exposures (eg, shared sources of chemicals) without bias from multicollinearity, however, there is a loss of information, generalizability (due to dimensionless units), and interpretability of the individual principal components that makes it difficult to identify specific ‘bad actor’ chemicals or their doses of harm when modeling the effects on health outcomes.⁴⁸ Other dimension reduction approaches that do consider the outcome (ie, are supervised) include shrinkage methods. Ridge regression shrinks the regression coefficients in a manner that decreases variance (with a trade-off of increased bias),⁴⁹ but it keeps all predictors in the model and thus does not reduce the dimensions of the data.^31,50 Lasso regression builds on ridge to allow coefficients to shrink to exactly zero, thus producing a model with fewer predictors.⁵¹ However, lasso saturates at n predictors (no more than the sample size), and it tends to select one exposure arbitrarily from a group of highly correlated exposures,⁵⁰ which may erroneously lead the researcher to conclude that the other unselected exposures are not associated with the health outcome.³¹ Elastic net combines both lasso and ridge regression and can work in a high-dimensional setting; however, it encourages a “grouping effect” that either keeps all or eliminates all the chemicals in a highly correlated set, even when some of the chemicals in the set may only be correlated due to shared exposure routes not health outcomes.³¹ This demonstrated poor specificity in the selection of correlated variables by elastic net compared to WQS is not ideal for studying chemical exposures.^31,52 Adaptive elastic net has been proposed to be more successful than its predecessor in the case of p ≫ n, but it does not provide a cumulative estimate across individual exposures to address research questions focused on combined mixture effects.⁵³

At its base, WQS calculates a mixture effect by combining quantiles of all the exposures (here, chemicals) into one summary weighted index that represents the cumulative mixture effect on the health outcome in a single direction.^31,52 To do so, in a training set of the data (eg, 40% of participants), each chemical is assigned a weight that reflects its individual contribution to the overall mixture effect on the outcome, where all weights sum to one. Every participant then has a score (mixture index) based on this formula for the weighted sum of their personal chemical exposures. In the final regression model using the testing set of the data (eg, the other 60% of participants), WQS simply determines the association of the mixture index with the outcome in a single degree-of-freedom, highly statistically powerful test. Quantiles are used in WQS instead of the fully continuous data because they convert the exposures to all be on the same scale in the mixture regardless of units and avoid extreme values that would cause weights to grow in extremes.³¹

It is best practice to always apply the repeated holdouts (RH) variation of WQS, which randomly partitions the data observations into the training set (used to estimate weights) and the testing/validation set (used to test the mixture effect) many times over.⁴⁶ The entire WQS model is repeated within each of, say, 100 repeated holdouts. This produces a distribution of validated results, which improves data representativeness and generalizability, stabilizes the chemical weights, and characterizes uncertainty in the estimates through confidence intervals.

For the critical extension to untargeted exposome data, the random subsets (RS) implementation of WQS addresses the challenges from having more exposure parameters than sample size (p ≫ n).⁵⁴ In the basic version of WQS, chemical weights are estimated by bootstrapping over observations in the training data (eg, 1,000 randomly selected bootstrap samples with replacement) and then determining the average weights based on the relative signal of the test statistic for the mixture index slope (in association with the outcome) in each bootstrap sample, using the full set of chemicals each time.³¹ By contrast, the RS version of WQS uses “feature bagging” to estimate weights over many different randomly selected subsets of chemicals in the full training data (Figure 1). This allows WQS_RS to aggregate across de-correlated sets of chemicals by repeatedly perturbing the exposure data under different correlation scenarios, thus avoiding multicollinearity and co-confounding, improving generalizability, and preventing overfitting.⁵⁴ The RH and RS variations should be used together, which performs the random subsetting procedure within every repeated holdout. Chemical weights are thus estimated a total of RH x RS times, so the computational intensity of WQS_RS can grow quickly. In summary, WQS_RS represents a major advancement in mixture methods for high-dimensional exposure data with complex correlation patterns.

Although WQS_RS supports high-dimensional untargeted data, inputting a larger and larger number of chemical features as exposures into the model does not necessarily improve the usefulness of the model. A summary of our roadmap is provided in Figure 2. First, depending on the research question, it may be beneficial to focus the mixture index to only environmental risk factors by excluding possible endogenous or pharmaceutical metabolites and early biomarkers of disease that could overinflate the mixture effect. This can be achieved through a variety of approaches, each of which faces a tradeoff between either excluding more possible endogenous metabolites or including more unidentified chemical features of unknown origin. For example, the chemical mixture could focus on only the detected features with annotations as possible environmental chemicals based on database matches, thus minimizing the presence of endogenous features (whether identified or unidentified) but being less inclusive to novel environmental chemicals not yet existing in databases. The coverage offered by a particular annotation database may need to be pruned depending on the research hypothesis; for example, the Norman Substance Database combines multiple sources of environmental chemical lists, some of which are pharmaceutical-focused or cover drinking water contaminants that include pharmaceuticals and hormone-related compounds.⁵⁵ Further, any additional annotated adducts or isotopes of the detected environmental chemicals could be removed to reduce feature redundancy (retaining the primary M + H or M-H adduct), although this should only be done among the annotations meeting the highest confidence level of three by xMSannotator (ie, adducts predictably clustered into the same correlation modules, retention time sub-modules, and mass-defect sub-groups).⁵⁶ An alternative approach for environmental chemical selection could focus on including all the chemical features that do not have identities or annotations as possible endogenous features, which would be inclusive to more chemicals but vary in effectiveness depending on the metabolomic database’s degree of coverage and endogenous/exogenous distinction. The strictest approach would be to focus the mixture on specific known classes or source groups of environmental chemicals. The potential interference by endogenous biomarkers is more pronounced in data from liquid chromatography (LC) HRMS, which widely integrates both the metabolome and the chemical exposome in its measurement of polar molecules with specific functional groups. Depending on sample preparation and extraction methods, gas chromatography (GC) HRMS can allow more focused detection of environmental chemicals, which is why the use of both instruments has been recommended for optimal chemical coverage.^10,15 In summary, the purpose of the WQS index and what mixture it intends to represent should be carefully considered in advance of modeling.

Second, the more chemicals included, the more repetitions needed, so the computational intensity can grow rapidly. Depending on the number of detected environmental chemicals in the research study, a layered pre-filtering strategy may help. The first layer is to restrict chemicals by some detection rate threshold to ensure there is sufficient spread of exposure across the quantiles for each chemical. Stricter thresholds (such as 50% or 75%) may miss out on important chemicals that are mostly found only in those that develop the disease, so this should be tailored to the study; in our 1:1 matched case-control study, for example, we are using 25% to allow for chemicals that hypothetically only present in half of the cases and none of the controls. However, quantizing data with such low detection thresholds requires a specialized approach (see next section).

As the next pre-filtering criterion, the mixture inputted into WQS_RS could be narrowed to only exposures with potential health relevance, while being cautiously inclusive. This is because including higher numbers of irrelevant chemicals may unexpectedly attenuate the mixture effect if the random subset parameter is not equivalently scaled, despite that the model trends weights of unimportant chemicals close to zero. To illustrate this, in the hypothetical case of a model with 18,000 chemicals compared to a model with 1,800 relevant chemicals, if the number of random subsets of chemicals remains the same (eg, 2,000 repetitions, even with slightly larger subset sizes), there would likely be more random subsets with mixtures full of the irrelevant chemicals that exhibit null effects, which then get averaged into the overall weighted index. Then, weight estimates may be based more on noise than signal. Instead, the mixture could first be further filtered to only the chemicals univariately associated with the outcome under a loose significance threshold without correcting for multiple testing (such as unadjusted p < 0.10 or even p < 0.20 in the adverse direction) in an exposome-wide association study (EWAS). We err on the side of inclusivity since univariate analyses are prone to issues from multicollinearity or chemical co-confounding and, under stricter significance thresholds, might miss chemicals that would show a stronger effect in the de-correlated mixture model. Although not without limitations, this pre-filtration approach has advantages over high-dimension shrinkage methods. As we described in the introduction, elastic net regression tends to either keep all or eliminate all chemicals among a highly correlated group,³¹ and thus could over-exclude potentially relevant chemicals from further consideration. However, future work should use simulated data to compare different approaches to selecting high-dimensional mixtures.

Finally, a third possible filtering layer is to conduct the WQS_RS models separately for the data from different instruments (GC-HRMS Versus LC-HRMS) or different instrument columns (LC hydrophilic interaction [HILIC] chromatography versus C18 hydrophobic interaction chromatography), if applicable. This not only reduces potential overlap in chemicals detected through each methodology but can also demonstrate robustness of results across platforms. However, the tradeoff is not identifying the fullest cumulative mixture effect from incorporating all relevant chemicals. Overall, the decisions about the chemical mixture to input into WQS_RS should depend on the research question and the dimensions of the exposome data in the specific study, and conducting sensitivity analyses under different mixture sizes and model parameters is helpful to understand the consistency of results.

For the choice of the number of repeated holdouts of observations to implement, 100 RHs are typical and sufficient to improve generalizability over traditional epidemiologic analyses and to characterize uncertainty of weights in mixture indices.⁴⁶ Although more RHs (such as >1,000) would approximate the normal distribution better, the computational requirements could be prohibitive: 1,000 RHs would take 10 times as long as 100 RHs would. By contrast, the choice of the number and size of random subsets of chemicals is more important to customize to each research question. We recommend first trying various choices of RSs under RH = 1 (for lower computational intensity) to see how the distribution of chemical weights change. If the number of chemicals with weights extremely near zero decreases with more RSs, that might mean that some chemicals did not have sufficient chance of being included under the lower number of RSs and thus that a higher number of RSs is still offering benefits. The size of each RS by default is the square root of the number of chemicals included in the mixture, and its input has a trade-off between giving a particular chemical more chances to be included in the model (under a fixed number of RSs) versus better de-correlating the data. In sensitivity analyses in our fertility study, we found the results to be sensitive to the choice of RS size. When we increased the number of chemicals by loosening filter criteria to include potentially less relevant exposures, and the RS size increased automatically, there were fewer important mixture contributors meeting Busgang criteria than when fixing the RSs at the original smaller size.²¹ We recommend to choose a relatively small size of each RS that allows for sufficient perturbation of chemical correlation patterns to discover important exposures, assuming that the number of RS repetitions is high enough to still give each chemical enough chance of being included in RSs. The previous WQS simulation study of the random subset implementation used a mixture of 472 untargeted metabolites with 1,000 RSs and a default RS size of 22 (ie, &cenveo_unknown_entity_Symbol_F0D6;472),⁵⁴ but simulation analyses have not yet been conducted on implications of higher RS sizes due to larger high-dimensional mixtures.

The average weights of chemicals are estimated based on the relative signal of the test statistic for the index slope in each random subset, and this signal function parameter has multiple possible values. For example, an “expt” signal function would apply the exponential of the t-statistic, which allows the most important chemicals to be much higher weighted than others, if that is desired. In addition, the ß slopes of the indices can be constrained to a single direction, which excludes random subsets with slopes in the opposite direction and ensures that only those in the relevant direction will contribute to the estimation of weights; thus, the mixture effect represents the adverse direction of harm.

There are two key deviations from default parameters that are important to consider. First, we recommend that the exposure data are manually quantized into deciles or quartiles before input into the WQS model (with the q parameter then set to null) in such a way that the non-detect values are put into their own zero-quantile and the detected values are quantized separately (eg, into 9 quantiles for a total of 10 deciles). The quantization is done individually for each chemical and modified from the source code for function gwqs_rank (see example code at github.com/anna-s-young/exposome-statistics). This approach avoids chemicals having, for example, five deciles that all refer to zero values and then large jumps between the latter five deciles. Second, matched case-control studies may need to consider how data are partitioned. The default is to randomly split the observations (eg, participants) into separate training and validation sets for each repeated holdout. However, this may lose some of the benefits of individually matched case-control pairs if a particular set no longer has a similar distribution of confounders between cases and controls. An alternative is to manually partition the pairs, instead of individual observations, randomly within each repeated holdout. This is supported by the “validation_rows” parameter in the model that takes a list of vectors (one for each RH) indicating the rows of observations to include in the validation set (see code available at github.com/anna-s-young/exposome-statistics).

The WQS_RS models with high-dimensional exposome data should be performed on a high-performance computing (HPC) cluster rather than a standard computer due to the computational intensity. The speed can be accelerated through parallel processing by using the “multisession” or “multicore” option for the plan strategy in the WQS function call. Although computational time can vary greatly for a cluster job depending on factors such as network performance, processing speed, cluster utilization, and available memory, in our unpublished lymphoma nested case-control study we completed most WQS sensitivity models in approximately 7-60h each under multisession (using mixtures of up to 2,700 chemicals, up to 4,000 RSs, 100 RHs, 444 samples, and our slower customizations of manual quantization and manual partitions) on the HPC at the Emory Rollins School of Public Health. Without the manual partitioning of case-control pairs into repeated holdouts, the models took an order of magnitude less time. In our fertility study, we also completed most of the sensitivity models in 4-50h each (using mixtures of up to 2,700 chemicals, up to 4,000 RSs, 100 RHs, 82 samples, and manual quantization).²¹ Increasing the number of chemicals and/or numbers of RSs to much higher degrees were our primary limiting factors for cluster load. However, we observed in the sensitivity analyses that the margins of return diminished with increasing numbers of RSs after a certain point, which justified us keeping 2,000 RSs for our main models (although this number will change depending on the study and mixture size). As described earlier, the choice of RSs can first be decided upon while using RH = 1 in faster preliminary models. Importantly, the independence of the repeated holdouts and random subsets makes them easily parallelizable and able to be scaled up in line with available high-performance computing resources.

With repeated holdout validation, the mixture effect (ß estimate for the index’s association with the outcome) is interpreted as an aggregation across repeated holdouts. The function will provide a mean and median of the mixture effect, along with confidence intervals. The median and percentile-based interval are preferred because they do not make assumptions about symmetry. The units of the mixture effect are per quantile. However, a histogram of the distribution of the cumulative mixture index across observations may show that it bunches in the middle quantiles, as it may be rare for a participant to have consistently low (or consistently high) quantiles for most exposures. If so, the per-quantile mixture effect may represent a difference covering a wide range of the exposures. An alternative is to transform the units to a per standard deviation (SD) change in the index. To do so, the ß estimate within each repeated holdout must be manually extracted from the output and multiplied by the standard deviation of that holdout’s index, then the median or mean of the estimates is used as the mixture effect per SD increase. This approach may also aid in the comparison of the cumulative mixture effect to the magnitudes of individual chemical effects, which can be estimated (and transformed to per-SD) in basic univariate regression models among the chemicals deemed important. Such a comparison can show the usefulness of WQS_RS for evaluating chemicals as mixtures instead of as single exposures.

The weights estimated for each chemical component can reveal potential ‘bad actor’ chemicals driving the overall mixture effect. However, because the weights always sum to one, they should only be interpreted in the case when the mixture effect itself is significant or borderline significant, based on its overall p value or percent of repeated holdouts in which it reaches significance. In the simplest interpretation of weights, chemicals are considered important if they have an average weight higher than the equi-weight threshold (1/p), which represents the hypothetical scenario of equal contributions by all chemicals. However, with the repeated holdout validation, we can better characterize the uncertainty in the chemical weights. Following the “Busgang criteria”, chemicals with an average weight above the threshold within at least 90% of repeated holdouts can be defined as “probable contributors”, else within at least 50% of holdouts as “possible contributors”, else within at least 10% of holdouts as “possibly not contributors”, else within less than 10% of holdouts as “probably not contributors.”^57,58 These criteria characterize how replicable the results are across data samples from the same underlying population, especially because certain chemicals may be misclassified as concerning—or not concerning—when looking at only one holdout of data.⁴⁶ We can graph these distributions of weights (eg, among the “possible contributors”) using the individual repeated holdout data from the WQS output. An example for chemicals across different levels of contribution is shown in Figure 3 (code available at github.com/anna-s-young/exposome-statistics).

WQS_RS is a powerful approach to identify cumulative mixture effects and bad actor chemicals while uniquely embracing the complexity and high dimensionality of untargeted exposome data. To our knowledge, it is currently the only mixture method that can be tailored to epidemiologic studies where the number of exposures exceeds the sample size. Compared to traditional univariate regression analyses, this mixture method avoids bias from multi-collinearity and co-confounding of chemicals and thus more accurately prioritizes chemicals of concern. At the same time, WQS_RS is highly statistically powerful, as it only conducts a single degree-of-freedom test, without losing information critical for interpretation of individual chemical risk factors. The summary mixture index not only provides a measure of cumulative mixture effects from the many chemicals that can simultaneously interfere with health, but it can also be used in other analyses as a single variable representing harmful exposure (see next section).

As limitations, WQS assumes that there are no interactions between exposures and that there is a constant change in risk between quantiles. These assumptions can be tested in sensitivity analyses considering other quantiles (eg, quartiles instead of deciles) and in other mixture methods considering small mixtures (such as Bayesian kernel machine regression).³⁹ With the random subsets implementation, WQS also assumes that all component effects are in the same direction, when in reality, the chosen chemical mixture may include substances that operate in the non-adverse direction, such as endogenous metabolites or some chemicals seeming to be protective due to unknown confounding (eg, pesticides related to diet and nutrition). Or, the health endpoint may be adverse in either extreme direction, as opposed to a simple dichotomy. In some cases, especially under hypothesis discovery, determining the direction of interest is not straightforward. Despite the limitation, this unidirectionality assumption actually helps prevent the reversal paradox arising from complex multicollinearity of exposures,⁴⁵ while also focusing the index for better interpretability. There is a recent extension of WQS for double (positive and negative) indices with a penalization term, but it does not yet support random subsets of chemicals for high-dimensional data at the same time.⁴⁵ As another limitation that should be acknowledged in matched case-control studies, WQS does not yet allow for conditional logistic regression in its models, only adjustment for the matching variables. WQS currently supports linear, logistic, Poisson, quasi-Poisson, and negative binomial regression in its current version 3.0.5 of the R package gWQS. Furthermore, WQS only supports continuous or ordinal exposure variables, however, categorical variables can be used if transformable to an ordinal structure. For example, previous work has evaluated quartiles of scales for post-traumatic stress disorder symptoms, depressive symptoms, and life stressors.^59,60 Although quantization of continuous exposures loses the full range of levels, it prevents extreme weights from outliers.³¹ Finally, the mixture method can become quite computationally intensive depending on the number of exposures, random subsets, repeated holdouts, and samples, and it is best used on a high-performance computing cluster.

Because of the simultaneous measurement of both the chemical exposome and metabolome using untargeted HRMS,^13-15 multi-omics analysis offers the opportunity to identify biological mechanisms that may underly the associations between chemical exposure and adverse health outcomes. For example, a meet-in-the-middle (MITM) strategy could: (1) identify the metabolic pathways associated with disease, (2) identify the metabolic pathways associated with exposure, and then (3) determine which significant pathways overlap “in the middle” between exposure and disease (Figure 4); however, researchers have approached MITM in different ways.⁶¹ Because the untargeted metabolomic data are also high dimensional, the supervised WQS_RS index is useful as a single exposure variable representing the cumulative mixture effect of untargeted chemicals on disease, thus reducing the complexity of one of the multi-omics layers and focusing the exposure index to the chemicals that are most relevant to the outcome and its mechanisms. Ideally, prospective, longitudinal samples would be used such that the metabolome data lie chronologically in the middle between the exposome and the health outcome and thus avoid problems with reverse causality. However, achieving this temporality is not always possible due to limited sample availability or budget restrictions, in which case care should be taken to acknowledge the potential for reverse causality or health treatment effects.

To accomplish steps 1 and 2 of MITM (separately), p values are first calculated for the univariate associations of each detected untargeted feature with the dependent variable (either the health outcome for step 1 or the exposure for step 2). Then, functional pathway enrichment analysis can predict functional activity by mapping all possible metabolite annotations for each feature (based on m/z and/or retention time) to a metabolic network and then finding the significant features that are locally enriched on a structure (ie, represent biological activity), whereas the false matches would only be randomly distributed in the network.⁶² This method leverages metabolic interconnections to improve prediction of pathway activity without having to identify each metabolite a priori.⁶² Mummichog (which is also implemented in MetaboAnalyst) has been the most common algorithm for pathway enrichment analysis,^62,63 however, we currently use Metapone’s permutation-based weighted hypergeometric test for several reasons.⁶⁴ Metapone jointly analyzes both positive- and negative-ion mode HRMS data when applicable (while avoiding double counting), accounts for matching uncertainty by using fractional counts of features (thus down-weighting those with higher numbers of matches), has an R package for suitability in our workflow alongside environmental mixture methods, and combines pathway information from multiple databases for higher relevance to xenobiotic pathways.⁶⁴ It leverages three databases: mummichog,⁶² the Kyoto Encyclopedia of Genes and Genomes (KEGG),⁶⁵ and the Small Molecule Pathway Database (SMPDB).⁶⁶ The overall p value of enrichment significance for each pathway is defined as the proportion of permutations in which the total fractional count of significant features in that pathway based on the real feature p values is lower than the total fractional count of significant features in that pathway based on randomly permuted p values for each feature.⁶⁴ As the final step in MITM, the results from pathway enrichment for exposure versus outcome can be compared for potential overlap (example in Figure 5).

When using Metapone, a high number of permutations, such as 1,000, will help stabilize results. The list of adducts ions to consider in annotations can be determined in consultation with the laboratory producing the data. The R package comes with a default data frame of pathway information, however, a data frame with a flag to filter to only the human pathways (flag = 1) is available at github.com/EMERGE-EXPOSOME/Metapone-pathway. The default threshold for significance of features within a pathway is raw p < 0.05, and the overall pathway p value has an optional adjusted value (by the Benjamini-Hochberg procedure) to account for multiple testing and reduce false positives. We recommend that significant pathways are also filtered to those with at least three significant metabolites (fractional counts) and that thus more comprehensively represent biological activity. The interpretations about specific pathways should be appropriately caveated due to the exploratory nature of this type of analysis. Finally, when conducting pathway enrichment using the WQS_RS mixture index, it is possible to determine significantly enriched pathways within each repeated holdout of WQS_RS (instead of using the average index) and investigate how frequently each pathway is significant; however, this would be computationally intensive.

While the MITM pathway enrichment analysis with WQS mixtures can reveal relevant mechanisms overall, the summary mixture indices may mask some significant biological pathways due to the diverse modes of action through which different chemicals operate. For this reason, it may also be useful to investigate the potential metabolic role of specific individual chemicals. For example, a few of the top contributors to the WQS mixture effect could be selected for additional pathway enrichment on their own. As another example, network analyses can explore correlations between the important chemicals in the WQS mixture effect (such as the “possible” and “probable” contributors) and each significant disease-associated metabolic pathway, where a pathway is represented by the first component (PC1) in principal component analysis (PCA) of the metabolites in that pathway (Figure 6). The PCA could be performed on only the significant metabolites within the pathway or, to understand how chemicals impact the pathway as a whole, on all the mapped metabolites regardless of significance. It is also possible to include in the network multiple PCs per pathway (such as the subset of PCs that explain most of the variance), or to use another method for creating a summary index of the pathway. It is important to note that correlation networks are exploratory and not suitable for causal interpretation, especially given the multicollinearity issues between chemicals, so further research would be required to confirm mechanisms of action of chemicals.⁶⁷

High-dimensional mediation analyses offer an alternative approach to investigate effects of exposure on disease through multiple metabolites (or preferably, pathway groups) and may improve causal inference depending on assumptions met.^67-70 Again, the WQS index could be helpful here as a single exposure variable representing the mixture effect, while the metabolome is retained as high-dimensional. The fact that WQS is a supervised approach that takes into account the outcome should be acknowledged when using the index to represent exposure. Other high-dimensional mediation research has employed separate exposure risk scores for different classes of chemicals,⁶⁸ however, this would be challenging with untargeted data where not all chemicals are identifiable or groupable. Pairwise mediation analyses between each chemical and each metabolite or pathway group are also possible to understand specific toxicant-mediator relationships, but they are again prone to bias from chemical multicollinearity unlike mixture effect indices.⁶⁸

Pathway enrichment is a practical exploratory approach to predict functional biological activity by leveraging pathway knowledge and bypassing the bottleneck of metabolite identification.^62,64 It produces interpretable and parsimonious results by mapping individual metabolites into pathway groups, which can generate hypotheses for future experimental research with the ultimate goal of identifying chemical mechanisms, therapeutic targets, or early biomarkers of disease. There are several limitations to note. Pathways are not mutually exclusive, and many metabolites are involved in multiple pathways. Thus, effects seen with a given metabolite do not always mean effects with all its pathways. The reliance only on m/z (and retention time in some cases) also limits the accuracy of metabolite annotations and may lead to false discovery, although the network mapping does help filter out the randomly structured and thus potentially irrelevant matches. In addition, the significance of enriched pathways does not reveal whether the involvement was harmful or beneficial to the health outcome, which is a challenge because of the different directions in which metabolites of the same pathway may act. Results are sensitive to the choice of databases, some of which overlook xenobiotic pathways,⁷¹ and to the selected thresholds for the significance and size of pathways.⁷² Furthermore, any pathway enrichment relies on known pathway definitions, which are subjective in their method of imposing order onto a biochemical network.⁷² Finally, if MITM is conducted on cross-sectional exposome-metabolome data, there is the possibility for reverse causation, which limits interpretation. In general, it is best to treat pathway enrichment results as exploratory.