---
title: "Introduction to tidyexposomics"
author: "Jason Laird"
header-includes:
   - \usepackage{amsmath}
   - \usepackage{amsfonts}
output:
  BiocStyle::html_document:
    toc_float: true
    toc: true
vignette: >
  %\VignetteIndexEntry{tidyexposomics: integrated exposure-omics analysis powered by tidy principles}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---



```{r setup,include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


Exposure to environmental factors is a major determinant of health and disease. The exposome is a term that represents the totality of environmental exposures (or in practical terms, the integrated compilation of all physical, chemical, biological, and (psycho)social influences that impact biology) from conception to death, that drive disease and overall health (Miller & Banbury Exposomics Consortium, 2025; Wild, 2005). Characterizing the relationship between a multiplicity of exposures and disease outcomes can be a daunting task. Epidemiological studies have moved towards exposure-wide association studies (ExWAS) where sets of exposures are associated with an outcome (Chung et al., 2024). The advent of high-throughput technologies has enabled profiling of several layers of biological information. These layers can also be integrated to better understand the relationship between exposures and disease outcomes.

The `tidyexposomics` package is designed to facilitate the integration of exposure and omics data to identify exposure-omics associations and their relevance to health outcomes.`tidyexposomics` extends the tidy-Bioconductor ecosystem (e.g., tidybulk, tidySummarizedExperiment) to exposome multi-omics integration using the MultiAssayExperiment container. It provides tidyverse-style accessors and functions for association testing, multi-omics integration, and ontology-driven enrichment, in an effort to complement existing tidy-Bioc tools.


```{r pkg-overview, echo=FALSE,out.width="120%",fig.align='center',fig.cap="tidyexposomics pipeline overview. QC, association testing, integration, and enrichment steps on a MultiAssayExperiment."}
knitr::include_graphics("./overview.png")
```

# Installation

```{r install, eval=FALSE}
# install the package
BiocManager::install("tidyexposomics")

# load the package
library(tidyexposomics)
```

# Command Structure
<hr>

To make the package more user-friendly, we have named our functions to be more intuitive. For example, we use the following naming conventions:

```{r command-str, echo=FALSE,fig.align='center',fig.cap= "Command naming conventions used throughout `tidyexposomics`. More complex pipelines begin with the `run_` prefix, visualizations with `plot_`, and data processing with `filter_`, `transform_`, `pivot_`, or `extract_` prefixes."}

knitr::include_graphics("./command_str.png")
```


We provide functionality to either add results to the existing object storing the omics/exposure data or to return results directly using `action = "get"`. We suggest adding results, given that pipeline steps are tracked and can be output to the R console, plotted as a workflow diagram, or exported to an Excel worksheet.

# Exposure Metadata and Ontology Annotation
<hr>

## Codebook Setup

Before starting an exposomics data analysis we recommend having a codebook, with information on your exposure variables. Some suggestions:

- **Variable Name**: The name of the variable in the data set.

- **Variable Description**: A concise description of what the variable measures, including units (e.g., “urinary bisphenol A (ng/mL)”).

- **Variable Type**: The type of variable, such as continuous, categorical, or binary.

- **Variable Period**: The period of time over which the variable was measured, such as “lifetime”, “year”, “month”, or “day”.

- **Variable Location**: The location where the variable was measured, such as “home”, “work”, “school”, or “geospatial code”.

- **Variable Ontology**: The ontology term associated with the variable.

## Ontology Choices

Variables captured in the codebook should be annotated with ontology terms to provide a standardized vocabulary for the variables. We recommend using the following ontologies for exposure and outcome variables:

- [Environment Exposure Ontology](https://www.ebi.ac.uk/ols4/ontologies/ecto) 
to annotate your exposure variables. 

- [Human Phenotype Ontology](https://www.ebi.ac.uk/ols4/ontologies/hp) to
annotate your outcome variables and phenotypic data.

- [Chemical Entities of Biological Interest](https://www.ebi.ac.uk/ols4/ontologies/chebi) 
to annotate your chemical exposure variables.

**Why annotate with ontologies?**

- **Interpretability**: Ontology labels clarify ambiguous or inconsistently named variables.

- **Harmonization**: You can compare and combine variables across datasets when they map to the same term.

- **Grouping**: Ontologies allow you to collapse fine-grained exposures into broader categories.

- **Integration**: Many public tools, knowledge graphs, and repositories are ontology-aware. This can make your results more interoperable and reusable.

## Ontology Annotation App

To help annotate exposure variables, we provide a lightweight shiny app:

```{r ontology-app,eval=FALSE}
# Launch the shiny app to annotate exposure variables
shiny::runApp(build_ont_annot_app())
```

**To use the app:**

- Click **Browse** to select your exposure metadata file.

- Then you can click the variable you’d like to link to an annotation term and search in the **Choose Ontology Term** dropdown.

- After you select a term, you will see a short description of the term.

- After you are done, click **Apply Annotate** to save the annotation.

- Now you can group exposures into larger categories by selecting each line and then choose your ontology and root depth level (where a lower number means a more general term).

- Then you can click **Apply Categorization** to apply the selected categorization to the selected rows.

- If the ontology has nothing to do with your variable, you may manually enter a category in the **Category** column. This will change the **Category Source** to manual and will not be linked to the ontology.

- Once you have annotated all your variables click **Download Annotated CSV** to save the annotated metadata file.


```{r ontology-app-img, echo=FALSE,out.width="150%",fig.align='center',fig.cap="Screenshot of the ontology annotation app. The sidebar has an upload button to load your exposure metadata file. The main panel displays the uploaded exposure metadata where users can select variables to annotate and categorize. The sidebar also contains buttons to apply annotations/categorizations and download the annotated file."}
knitr::include_graphics("./ont_annot.png")
```


# Loading Data
<hr>

To get started we need to load the data. The `create_exposomicset` function is used to create a `MultiAssayExperiment` object that contains exposure and omics data. As a quick introduction, a `MultiAssayExperiment` object is a container for storing multiple assays (e.g., omics data) and their associated metadata:


```{r mae-overview,echo=FALSE,fig.align='center',fig.cap="Overview of the MultiAssayExperiment structure linking samples, assays, and metadata."}
knitr::include_graphics("./mae_rep.png")
```

We use the MultiAssayExperiment object to store the exposure and omics data. The create_expomicset function has several arguments:

- `codebook`: is a data frame that contains information about the variables in the exposure metadata. The column names must contain variable where the values are the column names of the exposure data frame, and category which contains general categories for the variable names. This is the data frame you created with the ontology annotation app!

- `exposure`: is a data frame that contains the exposure and other metadata.

- `omics`: is a list of data frames that contain the omics data.

- `row_data`: argument is a list of data frames that contain information about the rows of each omics data frame.


We are going to start by loading in example data pulled from the [ISGlobal Exposome Data Challenge 2021](https://doi.org/10.1016/j.envint.2022.107422) (Maitre et al., 2022). Specifically, we will examine how exposures and omics features relate to asthma status in asthma patients with a lower socioeconomic status (SES).


```{r load-libraries,message=FALSE,warning=FALSE}
# Load Libraries
library(tidyverse)
library(tidyexposomics)
```


```{r load-data}
# Load example data
data("tidyexposomics_example")

# Create exposomic set object
expom <- create_exposomicset(
    codebook = tidyexposomics_example$annotated_cb,
    exposure = tidyexposomics_example$meta,
    omics = list(
        "Gene Expression" = tidyexposomics_example$exp_filt,
        "Methylation" = tidyexposomics_example$methyl_filt
    ),
    row_data = list(
        "Gene Expression" = tidyexposomics_example$exp_fdata,
        "Methylation" = tidyexposomics_example$methyl_fdata
    )
)
```

We are interested in how the exposome affects health outcomes, so let’s define which metadata variables represent exposure variables.

```{r exp-vars}
# Grab exposure variables
exp_vars <- tidyexposomics_example$annotated_cb |>
    filter(category %in% c(
        "aerosol",
        "main group molecular entity",
        "polyatomic entity"
    )) |>
    pull(variable) |>
    as.character()
```

# Quality Control
<hr>

## Missingness

Oftentimes when collecting data, there are missing values. Let’s use the `plot_missing` function to determine where our missing values are:

```{r missing-bar,fig.height=3,fig.width=6,fig.cap="Count of features with missing data above a 0% missingness threshold by data layer. Exposure data have variables with missingness."}

# Plot the missingness summary
plot_missing(
    exposomicset = expom,
    plot_type = "summary",
    threshold = 0
)
```

Here we see that there are 4 variables in the exposure data that are missing data. Let’s take a look at them:

```{r missing-bar-lollipop,fig.height=2,fig.width=4,fig.cap= "Percent missingness per exposure variable. Parity, `h_parity_None`, shows the highest missingness."}

# Plot missing variables withing exposure group
plot_missing(
    exposomicset = expom,
    plot_type = "lollipop",
    threshold = 0,
    layers = "Exposure"
)
```

Here we see that one variable, `h_parity_None`, has about 4% missing values. We can apply a missingness filter using the `filter_missing` function. However, given that this level of missingness is quite low, we will not be applying a missingness filter and instead impute the missing data.

## Imputation

The `run_impute_missing` function is used to impute missing values. Here we can specify the imputation method for exposure and omics data separately.

The `exposure_impute_method` argument is used to set the imputation method for exposure data, and the `omics_impute_method` argument is used to set the imputation method for omics data. The `omics_to_impute` argument is used to specify which omics data to impute. Here we will impute the exposure data given using the `missforest` method, but other options for imputation methods include:

- `median`: Imputes missing values with the median of the variable.

- `mean`: Imputes missing values with the mean of the variable.

- `knn`: Uses k-nearest neighbors to impute missing values.

- `mice`: Uses the Multivariate Imputation by Chained Equations (MICE) method to impute missing values.

- `missforest`: Uses the MissForest method to impute missing values.

- `lod_sqrt2`: Imputes missing values using the square root of the lower limit of detection (LOD) for each variable. This is useful for variables that have a lower limit of detection, such as chemical exposures.

```{r impute-missing}

# Impute missing values
expom <- run_impute_missing(
    exposomicset = expom,
    exposure_impute_method = "missforest",
    exposure_cols = exp_vars
)
```

## Filtering Omics Features

We can filter omics features based on variance or expression levels. The `filter_omics` function is used to filter omics features. The method argument is used to set the method for filtering. Here we can use either:

- **Variance**: Filters features based on variance. We recommend this for omics based on continuous measurements, such as log-transformed counts, M-values, protein intensities, or metabolite concentrations.

- **Expression**: Filters features based on expression levels. We recommend this for omics where many values may be near-zero or zero, such as RNA-seq data.


The `assays` argument is used to specify which omics data to filter. The `assay_name` argument is used to specify which assay to filter. The `min_var`, `min_value`, and `min_prop` arguments are used to set the minimum variance, minimum expression value, and minimum proportion of samples exceeding the minimum value, respectively.


```{r filter-omics}

# filter omics layers by variance and expression
# Methylation filtering
expom <- filter_omics(
    exposomicset = expom,
    method = "variance",
    assays = "Methylation",
    assay_name = 1,
    min_var = 0.05
)

# Gene expression filtering
expom <- filter_omics(
    exposomicset = expom,
    method = "expression",
    assays = "Gene Expression",
    assay_name = 1,
    min_value = 1,
    min_prop = 0.3
)
```

## Normality Check

When determining variable associations, it is important to check the normality of the data. The `run_normality_check` function is used to check the normality of the data.

```{r normality-check}

# Check variable normality
expom <- run_normality_check(
    exposomicset = expom,
    action = "add"
)
```

The `transform_exposure` function is used to transform the data to make it more normal. Here the transform_method is set to `boxcox_best` as it will automatically select the best transformation method based on the data. The `transform_method` can be manually set to `log2`, `sqrt`, or `x_1_3` as well. We specify the exposure_cols argument to set the columns to transform.

```{r trasform-vars}

# Transform variables
expom <- transform_exposure(
    exposomicset = expom,
    transform_method = "boxcox_best",
    exposure_cols = exp_vars
)
```

To check the normality of the exposure data, we can use the `plot_normality_summary` function. This function plots the normality of the data before and after transformation. The `transformed` argument is set to `TRUE` to plot the normality status of the transformed data.

```{r norm-plot,fig.width= 4,fig.height= 4, fig.cap= "Normality status of numeric exposure variables after Box-Cox transformation."}

# Examine normality summary
plot_normality_summary(
    exposomicset = expom,
    transformed = TRUE
)
```


## Principal Component Analysis

To identify the variability of the data, we can perform a principal component analysis (PCA). The `run_pca` performs a joint PCA across all numeric exposures and omic assays after standardization, identifying shared axes of variation across layers. The resulting PCs in `colData()` reflect integrated sample-level variance across all data types, and outliers are defined in that joint multi-omics PC space. 

Here we specify that we would like to log-transform the exposure and omics data before performing PCA using the `log_trans_exp` and the `log_trans_omics` arguments, respectively. We automatically identify sample outliers based on the Mahalanobis distance, a measure of the distance between a point and a distribution.


```{r pca}

# Perform principal component analysis
expom <- run_pca(
    exposomicset = expom,
    log_trans_exp = TRUE,
    log_trans_omics = TRUE,
    action = "add"
)
```



```{r pca-plot, fig.align='center', fig.width= 9, fig.height= 7,fig.cap= "PCA of sample and feature space with sample outlier detection."}

# Plot the PCA plot of sample and feature space
plot_pca(exposomicset = expom)
```

Here we see one sample outlier, and that most variation is captured in the first two principal components for both features and samples. We can filter out the outlier using the `filter_sample_outliers` function. 

```{r filt-outliers}

# Filter out sample outliers
expom <- filter_sample_outliers(
    exposomicset = expom,
    outliers = c("s1231")
)
```

To understand the relationship between the principal components and exposures we can correlate them using the `run_correlation` function. Here we specify that the `feature_type` is `pcs` for principal components, specify a set of exposure variables, `exp_vars`, and the number of principal components, `n_pcs.` We set `correlation_cutoff` to `0` and `pval_cutoff` to `1` to initially include all correlations.


```{r corr-pcs}

# Run the correlation analysis
expom <- run_correlation(
    exposomicset = expom,
    feature_type = "pcs",
    exposure_cols = exp_vars,
    n_pcs = 20,
    action = "add",
    correlation_cutoff = 0,
    pval_cutoff = 1
)
```


We can visualize these correlations with the `plot_correlation_tile` function. We specify we are plotting the `feature_type` of `pcs` to grab the principal component correlation results. We then set the significance threshold to 0.05 with the `pval_cutoff` argument.


```{r plot-pc-corr,fig.width= 7,fig.height= 2.2,fig.cap= "Correlation heatmap of exposures versus principal components. Child lead levels (`hs_pb_c_Log2`) and maternal BPA levels (`hs_bpa_madj_Log2`) are associated with the most principal components."}


# Plot the correlation tile plot
plot_correlation_tile(
    exposomicset = expom,
    feature_type = "pcs",
    pval_cutoff = 0.05
)
```


## Exposure Summary

We can summarize the exposure data using the `run_summarize_exposures` function. This function calculates summary statistics for each exposure variable, including the number of values, number of missing values, minimum, maximum, range, sum, median, mean, standard error, and confidence intervals. The `exposure_cols` argument determines which variables to include in the summary.


```{r exp-sum, fig.width=6, fig.height=4}

# Summarize exposure data
run_summarize_exposures(
    exposomicset = expom,
    action = "get",
    exposure_cols = exp_vars
) |>
    head()
```

## Exposure Visualization
<hr>

To visualize our exposure data, we can use the `plot_exposures` function. This function allows us to plot the exposure data in a variety of ways. Here we will plot the exposure data using a boxplot. The `exposure_cat` argument is used to set the exposure category to plot. Additionally, we could specify `exposure_cols` to only plot certain exposures. The `plot_type` argument is used to set the type of plot to create. Here we use a boxplot, but we could also use a ridge plot. 

```{r plot-aerosol, fig.width=4, fig.height=3.5, fig.cap="Distribution of aerosol exposures by sex."}

# Plot aerosol exposure distributions by sex
plot_exposures(
    exposomicset = expom,
    group_by = "e3_sex_None",
    exposure_cat = "aerosol",
    plot_type = "boxplot",
    ylab = "Values",
    title = "Aerosol Exposure by Sex"
)
```

Here we do not see any significant differences in aerosol exposure between males and females.

# Sample-Exposure Association
<hr>

## Sample Clustering

The `run_cluster_samples` function is used to cluster samples based on the exposure data, clustering approaches are available by setting the `clustering_approach` argument. Here we use the `dynamic` approach, which uses a dynamic tree cut method to identify clusters. Other options are:

- `gap`: **Gap statistic method (default)**; estimates optimal `k` by comparing within-cluster dispersion to that of reference data.

- `diana`: **Divisive hierarchical clustering (DIANA)**; chooses `k` based on the largest drop in dendrogram height.

- `elbow`: **Elbow method**; detects the point of maximum curvature in within-cluster sum of squares (WSS) to determine `k`.

- `dynamic`: **Dynamic tree cut**; adaptively detects clusters from a dendrogram structure without needing to predefine `k`.

- `density`: **Density-based clustering** (via `densityClust`); identifies clusters based on local density peaks in distance space.

```{r sample-clusters}

# Sample clustering
expom <- run_cluster_samples(
    exposomicset = expom,
    exposure_cols = exp_vars,
    clustering_approach = "dynamic",
    action = "add"
)
```

We plot the sample clusters using the `plot_sample_clusters` function. This function plots z-scored values of the exposure data for each sample, colored by the cluster assignment. The `exposure_cols` argument is used to set the columns to plot.


```{r plot-sample-clusters, fig.height=4, fig.width=12, fig.cap="Sample clustering heatmap using exposure profiles (z-scored). Clusters appear mostly driven by aerosol exposure during pregnancy."}

# Plot the sample clusters
plot_sample_clusters(
    exposomicset = expom,
    exposure_cols = exp_vars
)
```

Here we see two clusters, largely driven by particulate matter/aerosol exposure during pregnancy (`h_pm25_ratio_preg_None` and `h_pm10_ratio_preg_None`).

## Exposure Correlations

The `run_correlation` function identifies correlations between exposure variables. We set `feature_type` to `exposures` to focus on exposure variables and use a correlation cutoff of `0.3` to filter for meaningful associations. This cutoff can be adjusted based on your data and analysis needs.

```{r correlate-exposures}

# Run correlation analysis
expom <- run_correlation(
    exposomicset = expom,
    feature_type = "exposures",
    action = "add",
    exposure_cols = exp_vars,
    correlation_cutoff = 0.3
)
```


To visualize the exposure correlations, we can use the `plot_circos_correlation` function. Here we will plot the circos plot. This function creates a circular plot of the exposure correlations. The `correlation_cutoff` argument is used to set the minimum correlation score for the association. Here we use a cutoff of `0.3`.

```{r exposure-circos-corr, fig.height=10,fig.width=10,fig.align='center',fig.cap="Circos view of exposure-exposure correlations (threshold 0.3)."}

# Plot exposure correlation circos plot
plot_circos_correlation(
    exposomicset = expom,
    feature_type = "exposures",
    corr_threshold = 0.3,
    exposure_cols = exp_vars
)
```

## Exposure-wide association (ExWAS)

The `run_association` function performs an ExWAS analysis to identify associations between exposures and outcomes. We specify the data source, outcome variable, feature set, and covariates for the analysis. Since we have a binary outcome, we set the model family to binomial.


```{r assoc-exposures}

# Perform ExWAS
expom <- run_association(
    exposomicset = expom,
    source = "exposures",
    outcome = "hs_asthma",
    feature_set = exp_vars,
    action = "add",
    family = "binomial"
)
```

To visualize the results of the ExWAS analysis, we can use the `plot_association` function, which will plot results for the the specified features. The `terms` argument is used to set the features to plot. The `filter_thresh` argument is used to set the threshold for filtering the results. The `filter_col` argument is used to set the column to filter on. Here we use `p.value` to filter on the p-value of the association. We can also include the R^2 or adjusted R^2 (if covariates are included) using the `r2_col` argument.


```{r plot-exposure-assoc,fig.height=3, fig.width=5, fig.align='center',fig.cap="ExWAS associations of exposures with asthma status. No exposures are significantly associated (P < 0.05) with asthma status."}

# Plot the association forest plot
plot_association(
    exposomicset = expom,
    source = "exposures",
    terms = exp_vars,
    filter_thresh = 0.05,
    filter_col = "p.value",
    r2_col = "r2"
)
```

Here we see that no exposure variables are significantly associated with our asthma status. Although we do see that confidence interval for child Mono-iso-butyl phthalate (MiBP) levels (`hs_mibp_cadj_Log2`) does not cross 0, indicating a negative, albeit not significant (P < 0.05) association.

We can also associate our omics features with an outcome of interest using the `run_association` function. Here we specify an additional argument, `top_n`, which is used to set the top number of high variance omics features to include per omics layer.

```{r associate-top-omics,fig.height=3,fig.width=4.5,fig.align='center'}

# Perform ExWAS
expom <- run_association(
    exposomicset = expom,
    outcome = "hs_asthma",
    source = "omics",
    top_n = 500,
    action = "add",
    family = "binomial"
)
```


Now we can visualize these results with a manhattan plot.

```{r plot-manhattan,fig.height=7, fig.width=5, fig.align='center',fig.cap="Manhattan plot of omics-wide associations with asthma status."}

# Plot the manhattan plot
plot_manhattan(
    exposomicset = expom,
    min_per_cat = 0,
    feature_col = "feature_clean",
    vars_to_label = c(
        "TC19001180.hg.1",
        "TC01000565.hg.1",
        "cg01937701",
        "hs_mibp_cadj_Log2"
    ),
    panel_sizes = c(1, 3, 1, 3, 1, 1, 1),
    facet_angle = 0
)
```

## Exposome Scores
<hr>

We can also calculate exposome scores, which are a summary measure of exposure. The `run_exposome_score` function is used to calculate the exposome score. The `exposure_cols` argument is used to set the columns to use for the exposome score. The `score_type` argument is used to set the type of score to calculate. Here we could use:

- `median`: Calculates the median of the exposure variables.

- `mean`: Calculates the mean of the exposure variables.

- `sum`: Calculates the sum of the exposure variables.

- `pca`: Calculates the first principal component of the exposure variables.

- `irt`: Uses Item Response Theory to calculate the exposome score.

- `quantile`: Calculates the quantile of the exposure variables.

- `var`: Calculates the variance of the exposure variables.

The `score_column_name` argument is used to set the name of the column to store the exposome score in. Here we will define a score for aerosols using a variety of different methods and demonstrate their use in association with asthma status.

```{r calc-exposome, results='hide'}

# determine which aerosol variables to use
aerosols <- c("h_pm25_ratio_preg_None", "h_pm10_ratio_preg_None")

# Create exposome scores
expom <- expom |>
    run_exposome_score(
        exposure_cols = aerosols,
        score_type = "median",
        score_column_name = "exposome_median_score"
    ) |>
    run_exposome_score(
        exposure_cols = aerosols,
        score_type = "pca",
        score_column_name = "exposome_pca_score"
    ) |>
    run_exposome_score(
        exposure_cols = aerosols,
        score_type = "irt",
        score_column_name = "exposome_irt_score"
    ) |>
    run_exposome_score(
        exposure_cols = aerosols,
        score_type = "quantile",
        score_column_name = "exposome_quantile_score"
    ) |>
    run_exposome_score(
        exposure_cols = aerosols,
        score_type = "var",
        score_column_name = "exposome_var_score"
    )
```


We can then associate these exposome scores with asthma status using the `run_association` function, just like we did before. However, this time we specify our `feature_set` to be the exposome scores we just calculated.

```{r associate-exposome-score}

# Associate exposome scores with outcome
expom <- run_association(
    exposomicset = expom,
    outcome = "hs_asthma",
    source = "exposures",
    feature_set = c(
        "exposome_median_score",
        "exposome_pca_score",
        "exposome_irt_score",
        "exposome_quantile_score",
        "exposome_var_score"
    ),
    action = "add",
    family = "binomial"
)
```

```{r plot-exposome-scores, fig.height=2.5, fig.width=5,fig.cap="Associations of aerosol exposome scores with asthma status. The variance-based score has the strongest association with asthma status."}

# Plot the association forest plot
plot_association(
    exposomicset = expom,
    source = "exposures",
    terms = c(
        "exposome_median_score",
        "exposome_pca_score",
        "exposome_irt_score",
        "exposome_quantile_score",
        "exposome_var_score"
    ),
    filter_col = "p.value",
    filter_thresh = 0.05,
    r2_col = "r2"
)
```

# Differential Abundance
<hr>

## Differential Abundance

We provide functionality to test for differentially abundant features associated with an outcome across multiple omics layers. This is done using the `run_differential_abundance` function, which fits a model defined by the user (using the `formula` argument) and supports several methods. Here we apply the `limma_trend` method, a widely used approach for analyzing omics data. Users can also specify how features are scaled (e.g. none, quantile, TMM) before fitting. 

```{r diff-abundance,results='hide'}

# Run differential abundance analysis
expom <- run_differential_abundance(
    exposomicset = expom,
    formula = ~hs_asthma,
    method = "limma_trend",
    scaling_method = "none",
    action = "add"
)
```


We can summarize the results of the differential abundance analysis with a volcano plot, which highlights features with a high log fold change and that are statistically significant. The `plot_volcano` function generates this visualization, with options to set thresholds for p-values and log fold changes, and to label a subset of top-ranked features. In this example, we use the `feature_clean` column to display interpretable feature names.

**Note:** we set the `pval_col` to `P.Value` for the purposes of this example, but we recommend keeping the default of `adj.P.Val` to use the adjusted p-values.


```{r volcano-plot,fig.height=3.5, fig.width=6, fig.align='center',fig.cap="Volcano plot of differentially abundant features across omics layers."}

# Plot the volcano plot
plot_volcano(
    exposomicset = expom,
    top_n_label = 3,
    feature_col = "feature_clean",
    logFC_thresh = log2(1),
    pval_thresh = 0.05,
    pval_col = "P.Value",
    logFC_col = "logFC",
    nrow = 1
)
```

## Sensitivity Analysis

Depending on pre-processing steps, the results of the differential abundance analysis may vary. The `sensitivity_analysis` function is used to perform a sensitivity analysis to determine the robustness of the results. Here we determine if a feature is still differentially abundant if different minimum values, proportions, scaling methods are used, the inclusion of covariates, and after bootstrapping. We then define a stability score based on the number of times a feature is found to be differentially abundant under different conditions as well as the consistency of the effect size:

Where:


$$ Stability\ Score = \frac{\sum_i{(p_i < \alpha)}}{N} * \frac{1}{1 + \frac{\sigma_{\beta}}{\mu_{|\beta|}}}$$

Where: 

- $p_i$ is the p-value for the $i^{th}$ test

- $\alpha$ is the significance threshold

- $\beta$ is the log fold change

- $N$ is the number of tests

- $\sigma_{\beta}$ is the standard deviation of the effect size estimates

- $\mu_{|\beta|}$ is the mean of the absolute value of the effect size 
estimates. 

The first term captures the proportion of tests that are significant, while the second term captures the consistency of the effect size estimates. The stability score ranges from 0 to 1.

A stability score of 1 indicates that the feature is always found to be differentially abundant, while a stability score of 0 indicates that the feature is never found to be differentially abundant. Besides these, we provide other score metrics as well:

- `presence_rate`: Proportion of runs in which the feature’s p-value is below the specified threshold (selection frequency).

- `effect_consistency`: Inverse of the coefficient of variation of log fold-changes; measures effect size stability across runs.

- `stability_score`: Hybrid score combining presence_rate and effect_consistency, capturing reproducibility and signal strength.

- `mean_log_p`: Average of negative log-transformed p-values; represents overall statistical signal strength.

- `logp_weighted_score`: Product of `mean_log_p` and effect_consistency; highlights consistently strong features.

- `sd_logFC`: Standard deviation of log fold-change estimates; quantifies variability of effect sizes.

- `iqr_logFC`: Interquartile range of log fold-changes; provides a robust measure of effect size spread.

- `cv_logFC`: Coefficient of variation of log fold-changes; reflects relative variability of effect size.

- `sign_flip_freq`: Proportion of runs where the sign of the effect size differs from the overall average direction.

- `sd_log_p`: Standard deviation of log-transformed p-values; indicates variability in statistical signal.


```{r sensitivity-analysis,results='hide'}

# Perform Sensitivity Analysis
expom <- run_sensitivity_analysis(
    exposomicset = expom,
    base_formula = ~hs_asthma,
    methods = c("limma_trend"),
    scaling_methods = c("none"),
    pval_col = "P.Value",
    logfc_col = "logFC",
    logFC_threshold = log2(1),
    pval_threshold = 0.05,
    stability_metric = "stability_score",
    bootstrap_n = 10,
    action = "add"
)
```


Now we can plot the distribution of stability scores and visualize how many features have scores greater than a threshold of 0.34.

```{r plot-stability-score,fig.height=2.5, fig.width=8, fig.cap="Sensitivity summary highlighting robust features by stability score."}

# Plot distriubtion of stability scores
plot_sensitivity_summary(
    exposomicset = expom,
    stability_score_thresh = 0.34,
    stability_metric = "stability_score"
)
```

# Multi-Omics Integration
<hr>

While differential abundance analysis per omics layer can deliver insights into how each omic is associated with a particular outcome, we may want to leverage methods which integrate multiple omics layers. The `run_multiomics_integration` function is used to integrate multiple omics layers. Here we use either `MCIA`, `RGCCA`, `MOFA`, or the `DIABLO` method to integrate omics layers:

- `MCIA`: Multiple Co-inertia Analysis, a method that uses canonical correlation analysis to integrate multiple omics layers. We use the nipalsMCIA algorithm to compute the co-inertia scores from the [`nipalsMCIA` package](https://www.bioconductor.org/packages/release/bioc/html/nipalsMCIA.html) (nipalsMCIA).

- `RGCCA` : Generalized canonical correlation for flexible multi-block integration implemented using the [`RGCCA` package](https://rgcca-factory.github.io/RGCCA/) (Regularized and Sparse Generalized Ca...).

- `MOFA`: Multi-Omics Factor Analysis, a method that uses factor analysis to integrate multiple omics layers. MOFA is implemented using the [`MOFA2` package](https://www.bioconductor.org/packages/release/bioc/html/MOFA2.html) (MOFA2).

- `DIABLO`: Supervised multi-block PLS for outcome-aligned latent factors implemented using the [`mixOmics` package](https://www.bioconductor.org/packages/devel/bioc/html/mixOmics.html) (mixOmics).

Here we are interested in integrating our omics layers with the end goal of identifying multi-omics features that are associated with asthma status. So, we will use the `DIABLO` method.


```{r multiomics-int}

# Perform multi-omics integration
expom <- run_multiomics_integration(
    exposomicset = expom,
    method = "DIABLO",
    n_factors = 5,
    outcome = "hs_asthma",
    action = "add"
)
```
We can then use `plot_factor_summary` to visualize which omics contribute most to which factors.

```{r factor-summary, fig.height=2, fig.width=5, fig.align='center',fig.cap="Contribution of each omics layer to DIABLO latent factors."}

# Plot factor summary
plot_factor_summary(
    exposomicset = expom,
    midpoint = 2.5
)
```

Here we see that these factors are largely driven by features in the gene expression data.


## Factor Analysis

These methods are designed to identify factors that we can then associate with an outcome variable. Here we will use the `run_association` function to identify factors that are associated with asthma status.


```{r associate-factors}

# Identify factors that are associated with the outcome
expom <- run_association(
    exposomicset = expom,
    source = "factors",
    outcome = "hs_asthma",
    feature_set = exp_vars,
    action = "add",
    family = "binomial"
)
```

Now let’s see if any of our factors are associated with asthma status.

```{r plot-associate-factors,fig.width=5, fig.height=3,fig.cap="Associations of latent factors with asthma after covariate adjustment."}

# Plot the forest plot
plot_association(
    exposomicset = expom,
    source = "factors",
    filter_col = "p_adjust",
    filter_thresh = 0.05,
    r2_col = "r2"
)
```

We see that several factors are associated with our outcome. Factors have loading scores which indicate the strength of the association between the factor and the features. Here we can extract the top features, those in the 70th percentile, per omic, associated with our factors of interest. We set the `pval_col` and `pval_thresh` to filter our association results to grab the factors that pass our thresholds. However, we could specify specific factors with the `factors` argument too.


```{r extract-top-features}

# Extract top features that contribute to a factor
expom <- extract_top_factor_features(
    exposomicset = expom,
    method = "percentile",
    pval_col = "p_adjust",
    pval_thresh = 0.05,
    percentile = 0.7,
    action = "add"
)
```

We can visualize the top features associated with each factor using the `plot_top_factor_features` function. The `top_n` argument is used to set the number of top features to plot. The factors argument is used to set the factors to plot.


```{r plot-top-features,fig.height=8,fig.width=6,fig.align='center',fig.cap="Top-loading features per factor and omics layer."}

# Plot the top factor features
plot_top_factor_features(
    exposomicset = expom,
    top_n = 5,
    feature_col = "feature_clean"
)
```

# Exposure-Omics Association
<hr>

## Exposure-Omics Association

Now we have the option to correlate either the top factor features, differentially abundant features, or user-specified omics features (by using a variable map, a data frame with two columns, `exp_name` for the name of the omics assay, and `variable` for the name of the molecular feature) with exposures.

Here we will correlate features driving multiple latent factors with exposures. To grab these features, we will grab them from the metadata in our `MultiAssayExperiment` object. The `correlation_cutoff` is used to set the minimum correlation score, while the `pval_cutoff` is used to set the maximum p-value for the association.

```{r omics-exp-corr}

# Grab top common factor features and ensure
# feature is renamed to variable for the variable_map
top_factor_features <- expom |>
    extract_results(result = "multiomics_integration") |>
    pluck("top_factor_features") |>
    dplyr::select(variable = feature, exp_name)

# Correlate top factor features with exposures
# Perform correlation analysis between
# factor features and exposures
expom <- run_correlation(
    exposomicset = expom,
    feature_type = "omics",
    variable_map = top_factor_features,
    exposure_cols = exp_vars,
    action = "add",
    correlation_cutoff = 0.3,
    pval_cutoff = 0.05,
    cor_pval_column = "p.value"
)

# Perform correlation analysis between factor features
expom <- run_correlation(
    exposomicset = expom,
    feature_type = "omics",
    variable_map = top_factor_features,
    exposure_cols = exp_vars,
    feature_cors = TRUE,
    action = "add",
    correlation_cutoff = 0.3,
    pval_cutoff = 0.05,
    cor_pval_column = "p.value"
)
```

We can plot the results of the exposure-omics association analysis using the `plot_correlation_summary` function. Here we set the mode to `summary`, which will plot the number of associations per exposure and feature type. The `feature_type` argument is used to set the type of features to plot. Here we use `omics`, which are the top factor features.


```{r exp-omic-corr-summary,fig.height=7,fig.width=12,fig.align='center',fig.cap="Summary of exposure-omics and feature-feature associations."}

# Plot correlation summary
plot_correlation_summary(
    exposomicset = expom,
    mode = "summary",
    feature_type = "omics"
)
```


Here we note that gene expression features have a greater number of associations with exposures. Among exposures, we see that aerosols have the greatest number of associations with features, which may reflect the larger number of aerosol variables.

It may also be useful to identify which exposures are correlated with similar molecular features. We can do this with the `plot_circos_correlation` function. This function will plot a circos plot of the exposures and their shared features. The `feature_type` argument is used to set the correlation results to use for the analysis. Here we use the `omics` feature set, which plots the feature-exposure correlation information. The cutoff argument is used to set the minimum number of shared features to plot. Here we set the `shared_cutoff` to 1, which means that only exposures that share at least one feature will be plotted.

```{r exposure-shared-corr, fig.height=10,fig.width=10,fig.cap="Circos of exposures sharing correlated molecular features."}

# Plot the circos plot of shared correlated features
plot_circos_correlation(
    exposomicset = expom,
    feature_type = "omics",
    shared_cutoff = 1,
    midpoint = 1.5
)
```

Here we see that particulate matter exposures are correlated with many of the same features. 

## Network Analysis

The `run_create_network` function is used to create a network of exposures and omics features. The correlation results are used to create the network. The `feature_type` argument is used to set the type of features to use for the network. We could choose any of the following:

- `degs`: Differentially abundant features correlated with exposures.

- `factors`: Factor features correlated with exposures.

- `omics`: User specified omics features correlated with exposures.

Here we create networks based on two correlation tables generated above:

- `omics`: Correlations between omics and exposure variables.

- `omics_feature_cor`: Correlation just between the omics features.


This allows us to quantify a bipartite graph between exposures and omics features and a graph just between the omics features.

```{r create-network, fig.width=10,fig.height=7}

# Create omics within feature correlation network
expom <- run_create_network(
    exposomicset = expom,
    feature_type = "omics_feature_cor",
    action = "add"
)

# Create omics-exposure correlation network
expom <- run_create_network(
    exposomicset = expom,
    feature_type = "omics",
    action = "add"
)
```

To plot the network we can use the `plot_network` function. The network argument is used to set the type of network to plot. The `top_n_nodes` argument is used to set the number of nodes to plot. The `node_color_var` argument is used to set the variable to color the nodes by. The `label` argument is used to set whether or not to label the nodes. We can label the top n nodes, where the default is 5 nodes based on centrality. Additionally, we can choose to label certain nodes using the `nodes_to_label` argument. We can also include the network statistics using the `include_stats` argument.

```{r plot-network,fig.width=7, fig.height=7,fig.cap="Bipartite network of exposures and correlated omics features where node labels reflect most central nodes."}

# Plot the exposure-omics network
# Setting the seed so that the graph layout is consistent 
# given that layouts can have random starting points
set.seed(42)
plot_network(
    exposomicset = expom,
    network = "omics",
    top_n_nodes = 50,
    include_stats = TRUE,
    cor_thresh = 0.2,
    node_color_var = "group",
    label = TRUE,
    label_top_n = 5
)
```


Here we note that the network is highly connected, with particulate matter and methyl paraben exposures having the most connections to omics features. To briefly explain the metrics mentioned:

- **Nodes**: Number of nodes included in the graph.

- **Edges**: Number of connections between nodes in the graph.

- **Components**: Number of disconnected subgraphs.

- **Diameter**: The longest shortest path between any two nodes.

- **Mean Distance**: The average shortest path between any two nodes.

- **Modularity**: A measure of the division in the graph, where higher modularity indicates that nodes within the same community are more densely connected compared to nodes in different communities.

## Exposure-Omics Impact Network Analysis

The bipartite graph above describes how exposures are associated with omics features. However, this does not describe if exposures are associated with features that are more central in the feature-only graph (i.e. are exposures associated with features that are more central or peripheral to the graph?). To answer this question, we can use the `run_exposure_impact` function to calculate centrality metrics for the omics features associated with a given exposure. Centrality is computed on the feature-only graph. Centrality metrics are used to identify features that are more central to the network and may be more important in the context of the exposure-omics relationships. The centrality metrics included are:

- **Mean Betweenness Centrality**: Measures how often a node lies on the shortest path between other nodes where high values indicate potential intermediates.

- **Mean Closeness Centrality**: Measures how close a node is to all others in the network where high values indicate faster access to all nodes.

- **Mean Degree**: The average number of direct connections or edges a node has.

- **Mean Eigenvector Centrality**: Measures how well-connected a node is and how well-connected its neighbors are, where high values suggest influence in a well-connected cluster.


```{r exposure-impact}

# Run exposure-omics impact analysis
expom <- run_exposure_impact(
    exposomicset = expom,
    feature_type = "omics"
)
```

To plot the results of the exposure-omics impact analysis, we can use the `plot_exposure_impact` function. We set `feature_type` to `omics`, which means that we are plotting the network metrics of the omics features correlated with exposures. The `min_per_group` argument is used to set the minimum number of features per group to plot. This is useful for filtering out exposures that do not have enough associated features.

```{r plot-exposure-impact,fig.width=5.5,fig.height=5, fig.cap= "Network centrality metrics of omics features associated with each exposure."}

# Plot the exposure impact summary
plot_exposure_impact(
    exposomicset = expom,
    feature_type = "omics",
    min_per_group = 10
)
```


Here we see that child PFNA (`hs_pfna_c_Log2`) exposure is associated with omics features that are central to our network, despite not having the greatest number of features associated with it. This suggests that childhood PFNA exposure may have a significant impact on the omics features, even if it is not associated with as many features as some of the other exposures.


# Enrichment Analysis
<hr>

Enrichment analysis tests whether a set of molecular features (e.g. differentially abundant genes, metabolites, etc.) is over-represented in a predefined biological process. The benefit of grouping our exposures into categories is that we can now determine how broad categories of exposures are tied to biological processes. The `run_enrichment` function can perform enrichment analysis on the following feature types:

- `degs`: Differentially abundant features.

- `degs_robust`: Robust differentially abundant features from the sensitivity analysis.

- `omics`: User chosen features.

- `factor_features`: Multi-omics factor features either from `factor_type = “common_top_factor_features”` or `“top_factor_features”`.

- `degs_cor`: Differentially abundant features correlated with a set of exposures.

- `omics_cor`: User chosen features correlated with a set of exposures.

- `factor_features_cor`: Multi-omics factor features correlated with a set of exposures.

Here we will run enrichment analysis on factor features correlated with exposures. We specify `feature_col` to represent the column in our feature metadata with IDs that can be mapped (i.e. gene names). We will be performing Gene ontology enrichment powered by the [`fenr` package](https://www.bioconductor.org/packages/release/bioc/html/fenr.html) (Fenr, 2025). Note that we specify a `clustering_approach.` This will cluster our enrichment terms by the molecular feature overlap. 


```{r enrichment}

# Run enrichment analysis on factor features correlated with exposures
expom <- run_enrichment(
    exposomicset = expom,
    feature_type = c("omics_cor"),
    feature_col = "feature_clean",
    db = c("GO"),
    species = "goa_human",
    fenr_col = "gene_symbol",
    padj_method = "none",
    pval_thresh = 0.1,
    min_set = 1,
    max_set = 800,
    clustering_approach = "diana",
    action = "add"
)
```

## Enrichment Visualizations

To visualize our enrichment results we provide several options:

- dot`plot: A dot plot showing the top enriched terms. The size of the dots represents the number of features associated with the term, while the color represents the significance of the term.

- `cnet`: A network plot showing the relationship between features and enriched terms.

- `network`: A network plot showing the relationship between enriched terms.

- `heatmap`: A heatmap showing the relationship between features and enriched terms.

- `summary`: A summary figure of the enrichment results.

### Enrichment Summary

To summarize the enrichment results, we can use the `plot_enrichment` function with the `plot_type` argument set to `summary.` This will plot a summary of the enrichment results, showing:

- The number of exposure categories per enrichment term group.

- The number of features driving the enrichment term group.

- A p-value distribution of the enrichment term group.

- The number of terms in the enrichment term group.

- The total number of terms per experiment name.

- The overlap in enrichment terms between experiments (i.e. between gene expression and methylation).

```{r enrichment-summary,fig.width=14,fig.height=5, fig.cap="Summary of enriched GO terms grouped by overlap and exposure category."}

# Plot the summary diagram
plot_enrichment(
    exposomicset = expom,
    feature_type = "omics_cor",
    plot_type = "summary"
)
```

Here we see that it is just the features associated with “polyatomic entity” exposures that seem to be enriched. Additionally, there appears to be no overlap in terms between methylation and gene expression results.

### DotPlot

By setting the `plot_type` to dotplot we can create a dotplot to show which omics and which exposure categories are associated with which terms. By specifying the `top_n_genes` we can add the most frequent features in that particular enrichment term group.


```{r enrichment-dotplot,fig.height=11, fig.width=9, fig.cap="Dotplot of top enriched GO terms by omics layer and exposure category."}

# Plot a dotplot of terms
plot_enrichment(
    exposomicset = expom,
    feature_type = "omics_cor",
    plot_type = "dotplot",
    top_n = 15,
    add_top_genes = TRUE,
    top_n_genes = 5
)
```

### Term Network Plot

We can set the `plot_type` to `network` to understand how our enrichment terms are individually connected.

```{r enrichment-term-network, fig.width=6, fig.height=5, fig.cap="Network of enriched GO terms connected by shared genes."}

# Plot the term network plot
# Setting a seed so that the plot layout is consistent
set.seed(42)
plot_enrichment(
    exposomicset = expom,
    feature_type = "omics_cor",
    plot_type = "network",
    label_top_n = 1
)
```

At the individual term level, we see that they differ by omics layer, with the gene expression driving terms related to protein binding and the methylation data driving terms related to DNA binding and transcription.

### Heatmap

Setting the `plot_type` to `heatmap` can help us understand which genes are driving the enrichment terms. We have the additional benefit of being able to color our tiles by the Log_2_Fold Change from our differential abundance testing. Here we will examine group 1, given it seems to be driven by the most terms and multiple omics layers.


```{r enrichment-heatmap, fig.height=3, fig.width=8, fig.cap="Heatmap of genes driving enriched GO terms (Group 1) with log2 fold-change overlay."}

# Plot a heatmap of genes and corresponding GO terms
plot_enrichment(
    exposomicset = expom,
    feature_type = "omics_cor",
    go_groups = "Group 1",
    plot_type = "heatmap",
    heatmap_fill = TRUE,
    feature_col = "feature_clean"
)
```

Here we see that methylation features have the strongest fold change with respect to asthma status. Interestingly, we see that gene expression of *GZMH*, a potent immune effector molecule, is upregulated in asthma and correlated with exposures in the polyatomic entity category.

### Cnet Plot

Another way to visualize this information is with the `cnet` plot, where the enrichment terms are connected to the genes driving them.


```{r enrichment-cnetplot, fig.width=6, fig.height=6, fig.cap="Cnet plot linking enriched terms to contributing genes (Group 1)."}

# Plot the gene-term network
# Setting a seed so that the plot layout is consistent
set.seed(42)
plot_enrichment(
    exposomicset = expom,
    feature_type = "omics_cor",
    go_groups = "Group 1",
    plot_type = "cnet",
    feature_col = "feature_clean"
)
```

# Custom Analysis

We provide functionality to access the underlying data in the MultiAssayExperiment object and to construct tibbles for your own analysis:

- `pivot_sample`: Pivot the sample data to a tibble with samples as rows and exposures as columns.

- `pivot_feature`: Pivot the feature metadata to a tibble with features as rows and feature metadata as columns.

- `pivot_exp`: Pivot the sample and experiment assay data to a tibble with samples as rows and sample metadata as columns. Additionally, there will be a column for values for specified features in specified assays.

## Pivot Sample

Let’s check out the `pivot_sample` function. This function pivots the sample data to a tibble with samples as rows and exposures as columns.

```{r pivot-sample}

# Pivot sample data to a tibble
expom |>
    pivot_sample() |>
    head()
```

We could use this functionality to count the number of asthmatics per sex:

```{r pivot-sample-example}

# Count the number of asthma cases by sex status
expom |>
    pivot_sample() |>
    group_by(hs_asthma, e3_sex_None) |>
    reframe(n = n())
```

## Pivot Feature

The `pivot_feature` function pivots the feature metadata to a tibble with features as rows and feature metadata as columns. This can be useful for exploring the feature metadata in a more flexible way.

```{r pivot-feature}

# Pivot feature data to a tibble
expom |>
    pivot_feature() |>
    head()
```

We could use this functionality to count the number of features per omics layer, or to filter features based on their metadata. For example, we can count the number of features per omics layer:


```{r pivot-feature-example}

# Count the number of features per omic layer
expom |>
    pivot_feature() |>
    group_by(.exp_name) |>
    summarise(n = n())
```

## Pivot Experiment

Now if we want to grab assay data from a particular experiment, we can do that with the `pivot_exp.` Let’s try grabbing assay values for the `TC01004453.hg.1` probe (probe for *IL23R*) in the `Gene Expression` experiment:

```{r pivot-exp}

# Pivot experiment data to a tibble
expom |>
    pivot_exp(
        exp_name = "Gene Expression",
        features = "TC01004453.hg.1"
    ) |>
    head()
```

We can use this functionality to create custom plots or analyses based on the exposure and feature data. For example, we can plot the expression levels of *IL23R* by asthma status:

```{r pivot-exp-example, fig.height=4, fig.width=4, fig.align='center',fig.cap="IL23R gene expression levels (probe TC01004453.hg.1) by asthma status."}

# Plot expression of IL23R
expom |>
    pivot_exp(
        exp_name = "Gene Expression",
        features = "TC01004453.hg.1"
    ) |>
    ggplot(aes(
        x = as.character(hs_asthma),
        y = log2(counts + 1),
        color = as.character(hs_asthma),
        fill = as.character(hs_asthma)
    )) +
    geom_boxplot(alpha = 0.5) +
    geom_jitter(alpha = 0.1) +
    ggpubr::geom_pwc(label = "{p.adj.format}{p.adj.signif}") +
    theme_minimal() +
    ggpubr::rotate_x_text(angle = 45) +
    ggsci::scale_color_cosmic() +
    ggsci::scale_fill_cosmic() +
    labs(
        x = "",
        y = expression(Log[2] * "Abd."),
        fill = "Asthma Status",
        color = "Asthma Status"
    )
```

# Pipeline Summary

To summarize the steps we have taken in this analysis, we can use the `run_pipeline_summary` function. This function will provide a summary of the steps taken in the analysis. We can set `console_print` to `TRUE` to print the summary to the console. Setting `include_notes` to `TRUE` will include notes on the steps taken in the analysis.


```{r pipeline-summary}

# Run the pipeline summary
expom |>
    run_pipeline_summary(console_print = TRUE, include_notes = TRUE)
```


## Saving Results

We can export the results in our `MultiAssayExperiment` to an Excel spreadsheet using the `extract_results_excel` function. Here we add all of our results to the Excel file, but we can choose certain results by changing the `result_types` argument.

```{r save-results}

# Save results
extract_results_excel(
    exposomicset = expom,
    file = tempfile(),
    result_types = "all"
)
```


## Session Info

```{r session_info}

sessionInfo()
```

## References

Chung, M. K., House, J. S., Akhtari, F. S., Makris, K. C., Langston, M. A., Islam, K. T., Holmes, P., Chadeau-Hyam, M., Smirnov, A. I., Du, X., Thessen, A. E., Cui, Y., Zhang, K., Manrai, A. K., Motsinger-Reif, A., Patel, C. J., & Members of the Exposomics Consortium. (2024). Decoding the exposome: data science methodologies and implications in exposome-wide association studies (ExWASs). *Exposome, 4*(1), osae001. https://doi.org/10.1093/exposome/osae001  

fenr. (n.d.). Bioconductor. Retrieved August 18, 2025, from https://www.bioconductor.org/packages/release/bioc/html/fenr.html  

Maitre, L., Guimbaud, J.-B., Warembourg, C., Güil-Oumrait, N., Petrone, P. M., Chadeau-Hyam, M., Vrijheid, M., Basagaña, X., Gonzalez, J. R., & Exposome Data Challenge Participant Consortium. (2022). State-of-the-art methods for exposure-health studies: Results from the exposome data challenge event. *Environment International, 168*(107422), 107422. https://doi.org/10.1016/j.envint.2022.107422  

Miller, G. W., & Banbury Exposomics Consortium. (2025). Integrating exposomics into biomedicine. *Science, 388*(6745), 356–358. https://doi.org/10.1126/science.adr0544  

mixOmics. (n.d.). Bioconductor. Retrieved August 18, 2025, from https://www.bioconductor.org/packages/devel/bioc/html/mixOmics.html  

MOFA2. (n.d.). Bioconductor. Retrieved August 18, 2025, from https://www.bioconductor.org/packages/release/bioc/html/MOFA2.html  

MultiAssay Special Interest Group. (2025, April 15). *MultiAssayExperiment: The Integrative Bioconductor Container.* https://www.bioconductor.org/packages/release/bioc/vignettes/MultiAssayExperiment/inst/doc/MultiAssayExperiment.html  

nipalsMCIA. (n.d.). Bioconductor. Retrieved August 18, 2025, from https://www.bioconductor.org/packages/release/bioc/html/nipalsMCIA.html  

Regularized and Sparse Generalized Canonical Correlation Analysis for Multiblock Data. (n.d.). Retrieved August 18, 2025, from https://rgcca-factory.github.io/RGCCA/  

Wild, C. P. (2005). Complementing the genome with an “exposome”: the outstanding challenge of environmental exposure measurement in molecular epidemiology. *Cancer Epidemiology, Biomarkers & Prevention, 14*(8), 1847–1850. https://doi.org/10.1158/1055-9965.EPI-05-0456

