The environment-first approach starts from collecting and analyzing the extensive data arising from an individual’s environmental exposures before examining the complex internal biological responses.¹² According to Zhang et al. the environmental exposures can be organized into three main domains—the general external domain, the specific external domain, and the internal domain.¹⁸ The general external domain includes natural environments, socioeconomical, psychological factors, etc. The specific external domain is concerned with immediate local environment such as air, water, food, lifestyle, and occupation. The internal domain includes processes intrinsic to the body such as metabolic profiles, microbiome, inflammation status, etc, Comprehensively capturing exposure from specific external domain in a convenient and cost-effective manner is a fast-evolving research area, with many innovations and progresses made in recent years. For instance, environmental sensors provide objective, real-time measurements of factors like air and water pollutants, humidity, temperature, and allergens.¹⁹ With improved sensor technologies, granular data on environmental conditions can now be continuously collected. Recent studies integrating both targeted and non-targeted strategies are producing ever larger sets of data on associated chemical profiles in real-life matrices such as air, food, water or personal care products.^20-22 Wearable devices like smartwatches, fitness trackers, and specialized monitors allow constant physiological and contextual data gathering about activity, sleep, ultraviolet exposure, and other metrics, providing invaluable insights on exposures at personal level.¹⁰ Some devices can be combined with geospatial location data and health parameters, providing some direct bridge across datasets.²³ For example, silicon wristbands used for air pollutant monitoring can also be used to measure sweat metabolites.^24,25

Food and gut microbiome are strong influencers of host development and health trajectories through extensive metabolic interactions with the host.²⁶ The Earth Microbiome Project^27,28 and Human Microbiome Project^29,30 have significantly improved our knowledge on the microbiome. Many initiatives are producing and organizing data on food constituents such as FooDB, the Periodic Table of Food Initiative, Database on Migrating and Extractable Food Contact Chemicals.^31-33 Comprehensive chemical fingerprinting of foods is an active field of research. The vast diversity of natural products, macromolecules, and food processing-induced chemicals suggest that significant efforts are needed to map out all the chemical compounds in the foods. Chemical mixture exposome assessment remains an active research area. Many innovative concepts and approaches have been proposed, such as proactive use of the likelihood of co-exposure patterns depending on lifestyle and dietary habits, and effect-directed analysis to prioritize and assess the specific risk of identified mixtures.^34,35

Overall, the environment-first approach begins with an individual’s external or internal environmental exposure data to build the initial framework for unraveling the complex exposome. The current focus is on better characterizing food and microbiome, better capturing real-life chemical mixtures, individual exposure, and temporal variations in the exposome notably during critical windows of susceptibility to health hazards, such as childhood or pregnancy.³⁶

The phenotype-first approach starts by comprehensively profiling health outcomes before linking them to individuals’ environmental exposures, omics profiles and their interactions. Key sources for health outcomes data include EHRs, health surveys and various observable pathophysiological measures.^12,37 EHRs contain extensive documentation of an individual’s medical history, including diagnoses, medications, medical procedures, and other clinical information compiled over time across healthcare encounters. These records provide a rich context of an individual’s health trajectory over time. EHRs are mainly used in exposomics studies conducted in clinical settings. For studies focusing on environmental and public health, common phenotype data are health surveys administered by public health agencies and self-reported population-level data on health conditions, behaviors, lifestyle factors and environmental exposures collected by researchers.^38,39

Despite the wealth of data, effectively utilizing health outcome information presents significant challenges.⁴⁰ Patient privacy is a major concern given the ethical sensitivities around personal medical data. Heterogeneity also hinders synthesis of data from diverse EHR systems and surveys. The variations in terminology, format, structure, and coding make harmonizing disparate data a challenging task. Specialized platforms have been developed to address these challenges by enabling collaborative analysis of heterogeneous health data. For instance, the Observational Health Data Sciences and Informatics (OHDSI) provides an open-source framework to standardize and integrate diverse health datasets for evidence generation.⁴¹

Overall, the phenotype-first approach relies on accessing and synthesizing comprehensive health data sources through proper privacy protections and data integration techniques. Leveraging the recent developments in ontologies and artificial intelligence (AI) have the potential to transform this area.^42,43 For instance, the latest large language models (LLMs) such as GPT-4 and Gemini can read and interpret both texts and images within the same context. Their potential implementations in healthcare settings, such as to process and interpret EHR, have attracted considerable attention.⁴⁴

The omics-first approach starts with applying various omics technologies to obtain comprehensive profiles on individuals’ genotypes and biological responses associated with environmental exposures. Genomics encodes an individual’s genetic susceptibilities to certain exposures, environmentally induced epigenetic changes (eg, through DNA methylation) can be retained through lifetime or even trans-generationally as a proxy for exposures,^45-47 while miRNAomics, transcriptomics, proteomics, and metabolomics reveals corresponding expression changes at miRNA, mRNA, protein or metabolite levels.^48-50 Genomics are considered static and only need to be measured once. The downstream omics, often known as functional genomics, vary in response to environmental factors and operate at different time scales. They are usually sampled at multiple time points to depict the trajectory of the changes.

A key opportunity of the current multi-modality exposomics studies is integrated analysis to enable comprehensive understanding of exposure-response relationships.^51-53 For example, coupling exposomic profiling with epigenomics provides insights into how exposures may modulate epigenetic patterns underlying development and health outcome.⁵⁴ Integrating transcriptomic, proteomic and metabolomic data can reveal specific pathways modulated by exposures. At the very downstream of omics cascade, metabolomics occupies a unique place as the interface between individuals and environment exposures. Metabolomics platforms such as LC-HRMS can simultaneously measure both chemical exposures as well as metabolic responses from the same samples. In exposomics, metabolomics is often paired with microbiome profiling and food intake questionnaires to study the interactions between gut microbiome, dietary patterns, and the metabolome.^6,55-58

In summary, omics technologies have become well established across life sciences over the past two decades. They can efficiently measure host genotypes and capture downstream biological responses, with current focus on multi-omics, single-cell omics and spatial omics. Effective integration of these high-resolution omics profiles with environmental exposure data for translational discovery remains an active research area.⁵⁹

Exposomics data are primarily observational, therefore the typical dataset contains many covariates that must be accounted for during statistical analysis. The field of transcriptomics has a long history of using generalized linear models in gene expression analysis, which can flexibly handle many different experimental designs.^60-62 The methods have been extended to other omics types such as metabolomics,⁶³ proteomics,⁶⁴ microbiome data,⁶⁵ and in theory could be used to analyze any high-dimensional matrices. In this analysis framework, linear models can be configured to include both continuous and categorical variables, and categorical covariates can be modeled as either fixed or random effects.⁶³ For cases where there is only a chemical exposure matrix (E) and a sample metadata or phenotype matrix (P), a linear model can be configured to identify chemicals associated with a metadata of interest, while accounting for any number of covariates. For example, if the objective in a biomonitoring study is to identify chemicals associated with occupation type while accounting for relevant covariates, a potential model could be chemical_i ∼ occupation + age + sex + years of employment, where the latter three terms are covariates, and p-values and coefficients from the occupation term are used to identify the direction and statistical significance of chemical-occupation relationships. This method is considered univariate because each feature is analyzed independently.

In scenarios where omics data (O) are collected in addition to E and P, the objective may be to identify gene-chemical relationships. In this case, either the genes or the chemicals can be treated as sample metadata, and the same modeling approach described above can be used. Here, a potential model could be gene_j ∼ chemical_i + age + sex + BMI, which will generate results for i*j gene-chemical associations. Since a typical transcriptomics experiment generates measurements for ∼20 000 genes, this approach can produce an overwhelming number of results, and therefore it is recommended to apply stringent filters (including both statistical filters such as false discovery rate, and biological filters such as fold change or tissue specificity) to both the O and E matrices and focus only on the relationships of greatest interest. Walker et al. used this approach to identify compounds associated with occupational trichloroethylene (TCE) exposure, while adjusting for age, sex, and BMI.⁶⁶ Mass features associated with TCE exposure were carefully examined to separate TCE and TCE-derived compounds from metabolites that represent biological responses. This is a good example of how untargeted metabolomics data can simultaneously survey both the external chemical compounds and internal biological pathways. In the next section, we describe multivariate approaches that may be better suited to analyzing multiple high-dimensional matrices.

Dimensionality reduction (DR) methods summarize important sources of variation within high-dimensional datasets into a few component scores that are combinations of the original input feature values. There are many different DR methods, each of which produces sample-level scores and feature-level loadings. The most used method is principal component analysis (PCA), where sets of linear combination coefficients, also called ‘loadings’, are computed such that the resulting scores capture the maximum variance within the dataset. Different DR methods vary in the types of input features that they accept, the pre-processing steps that they employ, and the statistical criteria that they aim to maximize.⁶⁷ Most researchers use DR methods to visualize high-dimensional datasets in low-dimensional space such as 2D or 3D scatter plots. Data points are often annotated with different metadata using colors or text labels, to see if there are patterns correspond to specific experimental factors or covariates. DR can also be used as a feature selection tool. The magnitude of the loading coefficients indicates which features contribute most to each component and highly correlated features will have similar loadings. Therefore, extracting features with the highest magnitude for the top components will retrieve a set of features that are both highly correlated with each other, and that contribute a major source of variability within the dataset. If the scores plots showed that the component is related to a specific metadata variable, then the features in the top loadings can be interpreted as having a relationship with that metadata.

The same concept can be extended to support multi-omics data integration by considering two or more input matrices through joint dimensionality reductions (JDR). In this case, the component scores capture the shared variations across multiple datasets, and the component loadings capture sets of features that have correlated levels both within and across datasets. Multiple co-inertia analysis (MCIA) and multi-omics factor analysis (MOFA) are two commonly used JDR methods without considering sample metadata.^68,69 In comparison, the data integration analysis for biomarker discovery using latent components (DIABLO) is a supervised JDR method.⁷⁰ It takes a single metadata variable as input in addition to the two or more feature tables. Components are computed to maximize association with the metadata variable in addition to correlations between features both within and across the feature tables. In general, JDR methods require strictly matching samples and relatively large sample size (at least 20 per group) to ensure the identified patterns are robust.⁷¹ The top components from JDR are meaningful proxies of the original high-dimensional data by capturing their shared variance. They can be included in linear modeling to help identify important features underlay the observed associations with exposure or other metadata of interest. Effectively integrating JDR methods with linear modeling to reveal global patterns and the key underlying features represents a promising direction for developing computational platforms to support exploratory data analysis for exposomics.

Understanding the interplays of the exposome, metabolome and gut microbiome within the context of host genetics has been an active research area.^72,73 It is of great interest to be able to identify potential interactions between exposures and biological responses directly from their corresponding matrices. For example, detecting microbe-metabolite interactions from paired microbiome and metabolomics data have attracted significant research efforts in many recent studies.^74-77 The linear modeling and multivariate DR methods we discussed above can partially address this issue, but they have significant limitations. Linear models consider each feature individually, missing the opportunity to use the shared information across multiple features to boost performance. DR methods are notoriously difficult to interpret and are especially ill suited when the main goal is to obtain individual exposure-biomarker interactions where non-linear and higher-order interactions are expected to be common.¹² Advanced machine learning methods especially neural networks have shown great promise in this direction.

The mmvec is probably the first neural network algorithm applied to estimate microbe-metabolite interactions by learning the embeddings of microbial and metabolite features to estimates their co-occurrence probabilities.⁷⁴ Another neural network model, MiMeNet (Microbiome-Metabolome Network), uses a multi-layer perceptron (MLP) architecture to model metabolomic profiles based on microbial compositions.⁷⁶ The input layer takes the microbial composition profiles and maps them onto the hidden layers, where the neural network learns to extract complex patterns and representations from the input data through a series of mathematical functions, the final output layer produces the prediction of metabolomic profiles. This architecture enables the modeling of inputs and outputs of different sizes, which is the case for microbial and metabolomic features. A more recent algorithm, mNODE (metabolomic profile predictor using neural ordinary differential equations) has shown great promise in efficiently learning potential interactions from multi-omics data.⁷⁷ Similar to MiMeNet, mNODE adopts a MLP architecture, with an additional neural ODE module in between the hidden layers and employs different mathematical operations. The algorithm can incorporate dietary information to further enhance the performance in predicting metabolomic profiles from microbial compositions.

Although the “black box” characteristic of neural networks may be intimidating to some, their robustness in capturing complex patterns has been widely recognized especially in recent years. Compared to traditional modeling methods, they offer a distinct advantage in automatically learning and extracting relevant features from the data without having to necessarily define the exact number or nature of the features. With the growing acceptance of AI, the neural network-based approaches are expected to find broader applications in exposomics data analysis. In general, neural networks require a reasonably large experimental dataset, with sample size on the scale of hundreds and beyond, to train the model and tune the parameters. This usually will not be a concern for large-scale exposomics studies. When the sample size is small, linear modelling should be preferred due to its robustness and “white box” nature.

Networks are natural representations of complex relationships among different entities. Networks can be created based on our prior knowledge (ie, knowledge-based networks), based on the correlations computed from the current datasets (ie, data-driven networks), or using a mixture of both approaches. After creation, they can then be analyzed by applying graph theory to identify important connections and modules to gain insights into the organizing principles of the systems.⁷⁸ A key appealing feature for network-based approach is its intuitive visualization, which engages and empowers researchers to identify meaningful patterns. Networks can be viewed in hierarchical structures consisting of multiple layers corresponding to different components of exposomics datasets (Figure 2). Each omics layer can be explored separately and then jointly to gain mechanistic insights. We can intuitively find out the key points of connection between different omics layers as well as their interaction partners within the same layer to help understand pathways, biological processes, as well as to identify potential targets of intervention. This type of visualization can greatly facilitate data understanding and hypothesis generation, making networks a common choice for exposomics data integration and result presentation.⁷⁹

The knowledge-driven networks are created by projecting molecules of interest into a comprehensive knowledge network describing the known relationships (based on experimental evidence or computational prediction) among SNPs, genes, proteins, metabolites, diseases, and chemicals. The data-driven integrative networks can be created de novo based on the associations detected from current data using algorithms described in the previous sections. For instance, after covariate adjustments, the significant associations detected based on linear modeling or neural networks can be used to form networks between molecules, exposures and phenotypes. A key feature of the DIABLO-based multi-omics integration is its ability to generate highly interconnected networks with enriched biological themes.⁷⁰ In addition to visual exploration, different algorithms such as topology analysis, causal mediation analysis, integer linear programming can be applied to the networks to reveal key patterns, causal links, and mechanistic hypotheses.^80-82

Mendelian randomization (MR) has emerged as a valuable strategy in exposomics research to assess causal relationships between exposures and health outcomes by leveraging genetic data.^83-85 The core premise of MR is using genetic variants like SNPs as instrumental variables to infer the causal effect of a modifiable exposure on an outcome⁸⁶ This is possible because the random assortment of alleles during meiosis results in independent allocations of genetic variants. Therefore, associations between genetic variants and outcomes are less susceptible to confounding biases that often challenge conventional observational studies.⁸⁷

Two-sample MR (2SMR) approach has been proven very useful in exposomics by leveraging published data in MR analysis.⁸⁸ It utilizes summary statistics from comparable, separate genome-wide association studies (GWAS) for both exposure and outcome. These GWAS are often generated from large cohort studies, with sample sizes ranging from 1000 s to 100 000 s of individuals. For instance, a recent study examined the impact of 4587 environmental exposures on longevity using data from 361,194 individuals in the UK Biobank, and identified factors such as sugar intake, body fat, and various age-related diseases as being causally linked to longevity.⁸⁴ The study effectively demonstrated the capability of 2SMR to unravel relationships between exposures and health outcomes. The reliance on publicly available GWAS summary statistics not only enhances the cost-effectiveness of the 2SMR approach but also maximizes the utility of existing datasets, allowing for a broader examination of exposure-outcome relationships in an exposome-wide fashion.^89-91

The fundamental assumptions underpinning MR is that the genetic variants are strongly associated with the exposure, independent of confounders, and affect the outcome only through the exposure and not via alternative pathways.⁹⁰ A critical aspect of its design is the rigorous selection of appropriate instrumental variables.^92,93 Incorrect selection can lead to biased results, such as those arising from pleiotropy. Therefore, conducting sensitivity analysis is crucial to assess the robustness of the MR estimates against potential biases.⁹⁴ Moreover, careful interpretation of MR findings is necessary to avoid overstating causal conclusions based solely on observed associations.⁹⁵ To ensure the robustness and validity of findings, MR studies must control for Type I and Type II errors. Techniques such as the Bonferroni correction and False Discovery Rate (FDR) control are essential for minimizing the risk of false positives, while careful power analysis is important for reducing the likelihood of false negatives.⁹⁶

Results from causal inference in an observational setting require further experimental validations to establish causality. For complex mixtures, effect-directed analysis (EDA) can help uncover causal relationships between adverse health effects and chemical contaminants. This approach starts with toxicity screening of complex real-life mixtures. The complexity of the mixture is then reduced through fractionation steps and characterized using non-targeted approaches to narrow down on substances of concern.³⁵ Finally, dose-response studies are applied to the shortlist of candidate compounds to pinpoint the compounds that are cause bioactivity upon exposure.^97,98

Dose-response relationship is a foundational concept in toxicology and chemical risk assessment. Recent years have seen the growing applications of dose-response analysis to omics data, particularly transcriptomics data.⁹⁹ In this approach, exposures are conducted at multiple concentrations of a chemical, typically include a control group and at least three different dose groups with the same number of replicates in each group. The sample sizes vary from a few dozens to hundreds depending on the design, cost, and effect sizes. After data collection, a suite of linear and non-linear curve models are fitted to the levels of each omics feature. Each curve is analyzed to compute a feature-level benchmark dose (BMD), which is the minimum concentration of a substance that produces a clear, low level health risk relative to the control group.¹⁰⁰ The collection of benchmark doses can be summarized at the pathway level or across the entire experiment (commonly called the transcriptomic point-of-departure, or tPOD). When in vitro assays are coupled with high-throughput omics platforms, it becomes feasible for a typical research lab to use transcriptomics dose-response analysis to assess 10∼100 s of chemicals.^97,101,102 The omics-based dose-response study not only allows researchers to establish causal links between specific chemical exposures and general bioactivity (ie, which chemicals in the exposome ‘matter’), but also to disentangle specific biological processes that are perturbed by different chemicals.

The previous sections have introduced a handful of methods chosen based on their perceived utilities in exposomics data analysis. Beyond these, many other advanced statistical methods are also used for similar purposes, especially in the field of epidemiology.¹² General machine learning approaches such as random forest, support vector machines, and gradient boosting can be generally applied for biomarker selection and predictive modeling. Bayesian additive regression trees is a more recent method that has the added benefit of quantifying the uncertainty of predictions.¹⁰³ Other methods have been developed for assessing causal effects of complex mixture exposures on clinical outcomes and phenotypes, a blend of the objectives that we presented for 2SMR and EDA. Some examples include G-computation associated with super learners,¹⁰⁴ grouped weighted quantile sum (GWQS) regression,¹⁰⁵ Bayesian kernel machine regression,¹⁰⁶ etc,