Getting the data in order

To get started with an exposure-outcome-only workflow, we need to create our exposomicset which stores this data. This time, we are going to add in the omics data too!

# Load Libraries
library(tidyexposomics)

# 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
  )
)

Here you will see we have two lists: (1) the omics matrices, where the rows are omics features (e.g. genes, proteins, metabolites) and the columns are samples, (2) row data with information about these omics features. We also gave it a codebook with information about the exposures (see more about codebooks here), and the meta data which contain the exposures and covariate information.

For the sake of simplicity we will focus on a short list of exposures:

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

Missingness

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

# 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:

# 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.

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

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.

# 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.

# 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.

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

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.

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

## # A tibble: 6 × 27
##   variable   n_values  n_na   min   max range    sum median  mean    se ci_lower
##   <chr>         <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl>    <dbl>
## 1 h_pm10_ra…       48     0 16.6  25.8   9.25 1.01e3  21.3  21.0   0.31    20.4 
## 2 h_pm25_ra…       48     0 10.6  18.2   7.55 7.12e2  14.7  14.8   0.24    14.3 
## 3 hs_bpa_ma…       48     0  0     2.31  2.31 7.30e1   1.58  1.52  0.05     1.41
## 4 hs_mibp_c…       48     0  0.1   0.22  0.11 7.83e0   0.17  0.16  0        0.16
## 5 hs_pb_c_L…       48     0  1.44  5     3.56 1.59e2   3.32  3.32  0.1      3.11
## 6 hs_pfhxs_…       48     0  0     2.14  2.13 6.72e1   1.49  1.4   0.07     1.25
## # ℹ 16 more variables: ci_upper <dbl>, variance <dbl>, sd <dbl>,
## #   coef_var <dbl>, period <chr>, location <chr>, description <chr>,
## #   var_type <chr>, transformation <chr>, selected_ontology_label <chr>,
## #   selected_ontology_id <chr>, root_id <chr>, root_label <chr>,
## #   category <chr>, category_source <chr>, transformation_applied <chr>

Exposure Visualization

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.

# Plot aerosol exposure distributions by sex
plot_exposures(
  exposomicset = expom,
  group_by = "e3_sex_None",
  exposure_cols = c(
    "h_pm25_ratio_preg_None",
    "h_pm10_ratio_preg_None"
  ),
  plot_type = "boxplot",
  ylab = "Values",
  title = "Aerosol Exposure by Sex"
) +
  ggpubr::stat_compare_means(label = "p.signif")

Here we see that levels of PM2.5 were higher in females compared to males.

Coverting the Exposomicset back to dataframes

If you need to grab your codebook, feel free to do so using the extract_results function:

extract_results(
  exposomicset = expom,
  result = "codebook"
) |>
  head()

## # A tibble: 6 × 13
##   variable               period    location description  var_type transformation
##   <chr>                  <chr>     <chr>    <chr>        <chr>    <chr>         
## 1 h_pm10_ratio_preg_None Pregnancy Home     pm10 value … numeric  None          
## 2 h_pm25_ratio_preg_None Pregnancy Home     pm25 value … numeric  None          
## 3 hs_pm10_yr_hs_h_None   Postnatal Home     pm10 value … numeric  None          
## 4 hs_pm25_yr_hs_h_None   Postnatal Home     pm25 value … numeric  None          
## 5 h_mbmi_None            Pregnancy <NA>     Maternal pr… numeric  None          
## 6 e3_gac_None            Pregnancy <NA>     Gestational… numeric  None          
## # ℹ 7 more variables: selected_ontology_label <chr>,
## #   selected_ontology_id <chr>, root_id <chr>, root_label <chr>,
## #   category <chr>, category_source <chr>, transformation_applied <chr>

If you want to grab the exposure+metadata do so using the pivot_sample function:

expom |>
  pivot_sample() |>
  head()

## # A tibble: 6 × 19
##   .sample e3_sex_None hs_child_age_None h_age_None h_mbmi_None e3_gac_None
##   <chr>   <fct>                   <dbl>      <dbl>       <dbl>       <dbl>
## 1 s812    female                   6.84       35          33.2        41.4
## 2 s296    female                   6.75       32.8        25.1        39.1
## 3 s242    female                   6.97       33          24.6        40.4
## 4 s376    male                     6.35       27.1        28.3        37  
## 5 s1183   female                   6.74       32.7        33.3        38.1
## 6 s98     female                   6.93       23          31.2        39  
## # ℹ 13 more variables: h_native_None <fct>, h_parity_None <fct>,
## #   hs_asthma <dbl>, hs_zbmi_who <dbl>, h_pm10_ratio_preg_None <dbl>,
## #   h_pm25_ratio_preg_None <dbl>, hs_pm10_yr_hs_h_None <dbl>,
## #   hs_pm25_yr_hs_h_None <dbl>, hs_pb_c_Log2 <dbl>, hs_pfhxs_m_Log2 <dbl>,
## #   hs_pfna_c_Log2 <dbl>, hs_bpa_madj_Log2 <dbl>, hs_mibp_cadj_Log2 <dbl>

Exposome-wide association study (ExWAS)

Inspired by the concept of genome-wide association studies, or GWAS, the ExWAS is a high-throughput way of identifying many exposures with an association to an outcome. The run_association allows us to see which exposure variables are associated with our outcome, hs_asthma, which is a 0 or 1, so we set the family to binomial. If the variable was continuous, we’d usually set it to gaussian. Just one more note: since we are running this model on many different exposures, we use .feature as a placeholder:

# Perform ExWAS
expom <- run_association(
  exposomicset = expom,
  formula = hs_asthma ~ .feature,
  source = "exposures",
  feature_set = exp_vars,
  action = "add",
  family = "binomial"
)

We can visualize these associations using the plot_association function and set our p-value filter to color associations below that threshold, here we set it to 0.3 just to demonstrate the coloring, but we suggest using the p_adjust column for filter_col and a cutoff of 0.05:

# Plot the association forest plot
plot_association(
  exposomicset = expom,
  effect = "main",
  source = "exposures",
  terms = exp_vars,
  filter_thresh = 0.3,
  up_col = "#8E0152",
  down_col = "#006666",
  filter_col = "p.value",
  r2_col = "r2"
)

Interaction Modeling

Sometimes, the relationship between an exposure an outcome depends on another variable. For this folks tend to leverage interaction modeling. Let’s look at how we would see if the the relationship between asthma and an exposure depends on age:

run_association(
  exposomicset = expom,
  overwrite = TRUE,
  formula = hs_asthma ~ .feature * h_age_None,
  source = "exposures",
  feature_set = exp_vars,
  action = "add",
  family = "binomial"
) |>
  plot_association(
    effect = "h_age_None",
    source = "exposures",
    terms = exp_vars,
    filter_thresh = 0.3,
    up_col = "#8E0152",
    down_col = "#006666",
    filter_col = "p.value",
    r2_col = "r2"
  )

Notice, how we switched effect = "main" to effect = "h_age_None" so that we could grab the interaction results.

Stratifying Models

We may also want to run a model on a portion of our data, like say just male patients:

expom |>
  run_association(
    formula = hs_asthma ~ .feature + h_age_None,
    source = "exposures",
    overwrite = TRUE,
    feature_set = exp_vars,
    family = "binomial",
    strata = "e3_sex_None",
    name = "by_sex"
  ) |>
  plot_association(
    effect = "main",
    source = "exposures",
    model = "by_sex_male",
    terms = exp_vars,
    filter_thresh = 0.3,
    up_col = "#8E0152",
    down_col = "#006666",
    filter_col = "p.value",
    r2_col = "r2"
  )

We can also filter our object with the subset argument. Let’s examine the association in children with a higher BMI z-score:

expom |>
  run_association(
    formula = hs_asthma ~ .feature + h_age_None,
    source = "exposures",
    overwrite = TRUE,
    feature_set = exp_vars,
    family = "binomial",
    subset = hs_zbmi_who > 0,
    name = "high_bmi"
  ) |>
  plot_association(
    effect = "main",
    source = "exposures",
    model = "high_bmi",
    terms = exp_vars,
    filter_thresh = 0.3,
    up_col = "#8E0152",
    down_col = "#006666",
    filter_col = "p.value",
    r2_col = "r2"
  )

Comparing models

Now sometimes we want to compare the results of different models. Here we demonstrate how to make many different models using two approaches: (1) specifying a list of formulas, and (2) using the strata or subset functionality. This will create results for four different models:

expom <- expom |>
  run_association(
    formula = list(
      crude = hs_asthma ~ .feature,
      adj = hs_asthma ~ .feature + h_age_None
    ),
    source = "exposures",
    overwrite = TRUE,
    feature_set = exp_vars,
    family = "binomial"
  ) |>
  run_association(
    formula = hs_asthma ~ .feature + h_age_None,
    source = "exposures",
    feature_set = exp_vars,
    family = "binomial",
    strata = "e3_sex_None",
    name = "by_sex"
  )

To see what model results are available you can use the following:

association_models(expom)

## # A tibble: 5 × 8
##   source    model         formula     subset stratum family n_samples n_features
##   <chr>     <chr>         <chr>       <chr>  <chr>   <chr>      <int>      <int>
## 1 exposures default       hs_asthma … <NA>   <NA>    binom…        48          9
## 2 exposures crude         hs_asthma … <NA>   <NA>    binom…        48          9
## 3 exposures adj           hs_asthma … <NA>   <NA>    binom…        48          9
## 4 exposures by_sex_female hs_asthma … <NA>   female  binom…        28          9
## 5 exposures by_sex_male   hs_asthma … <NA>   male    binom…        20          9

Now let’s demonstrate how to compare two models, here we will compare effect estimates for males and females:

expom |>
  compare_association(
    "exposures",
    models = c(
      "by_sex_female",
      "by_sex_male"
    )
  ) |>
  plot_association_comparison(type = "forest")

We could also choose to visualize it using a scatter plot:

expom |>
  compare_association(
    "exposures",
    models = c(
      "by_sex_female",
      "by_sex_male"
    )
  ) |>
  plot_association_comparison(type = "scatter")

Here we see that the effect estimates are not correlated with one another, suggesting that the relationship between asthma and exposure differs by sex.

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.

# 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
)

Exposure-Omics Association

Above we saw that there are not too many omics features associated with asthma. Which may be due to the subsampling in this example or because exposures are driving different biology. Let’s examine what omics features exposures are associated with.

# Run association testing between every exposure and omics feature
expom <- run_exposure_omics_association(
  exposomicset = expom,
  exposures = exp_vars,
  action = "add"
)

Now let’s see how many omics features each exposure is associated with using the plot_exposure_omics_association function where we can either plot by the individual exposure or exposure category:

# Plot the number of exposure-omics associations
plot_exposure_omics_association(
  exposomicset = expom,
  plot_type = "category",
  pval_col = "p.value",
  pval_thresh = 0.05
)

Here we can see that there are omics features associated with exposures, while there are fewer that are associated with asthma directly.

Enrichment Analysis

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 omics features associated with exposures. Specifically, we will grab omics features associated with aerosols using the extract_results function. This function allows us to pull any of the results we have been generating so far. We will then filter these association results to significant associations (p-value < 0.05) and those with the category “aerosol”.

# Extract omics features associated with aerosols
var_map <- extract_results(
  exposomicset = expom,
  result = "association"
) |>
  purrr::pluck(
    "exposure_omics",
    "results_df"
  ) |>
  dplyr::filter(
    p.value < 0.05,
    category == "aerosol"
  ) |>
  dplyr::select(
    exp_name,
    variable = feature
  )

Now we will perform enrichment analysis and 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 (Fenr, 2025). Note that we specify a clustering_approach. This will cluster our enrichment terms by the molecular feature overlap.

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

Enrichment Visualizations

To visualize our enrichment results we provide several options:

• dotplot: 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.

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

## NULL

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

## NULL

# 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",
  plot_type = "network",
  label_top_n = 3
)

## NULL

# Build the fill table
# one estimate per experiment x feature
fill_df <- extract_results(
  exposomicset = expom,
  result = "association"
) |>
  purrr::pluck(
    "exposure_omics",
    "results_df"
  ) |>
  dplyr::filter(
    p.value < 0.05,
    category == "aerosol"
  ) |>
  dplyr::reframe(
    estimate = median(estimate),
    .by = c(exp_name, feature)
  )

# Plot a heatmap of genes and corresponding GO terms
plot_enrichment(
  exposomicset = expom,
  feature_type = "omics",
  go_groups = "Group 2",
  plot_type = "heatmap",
  fill_data = fill_df,
  feature_col = "feature_clean"
)

## NULL

# 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",
  go_groups = "Group 2",
  plot_type = "cnet",
  fill_data = fill_df,
  feature_col = "feature_clean"
)

## NULL

sessionInfo()

## R version 4.6.2 (2026-10-29)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] tidyexposomics_1.1.3        MultiAssayExperiment_1.39.1
##  [3] SummarizedExperiment_1.43.0 Biobase_2.73.2             
##  [5] GenomicRanges_1.65.4        Seqinfo_1.3.2              
##  [7] IRanges_2.47.5              S4Vectors_0.51.10          
##  [9] BiocGenerics_0.59.12        generics_0.1.4             
## [11] MatrixGenerics_1.25.0       matrixStats_1.5.0          
## [13] BiocStyle_2.41.0           
## 
## loaded via a namespace (and not attached):
##   [1] naniar_1.1.0         RColorBrewer_1.1-3   sys_3.4.3           
##   [4] rstudioapi_0.19.0    jsonlite_2.0.0       magrittr_2.0.5      
##   [7] farver_2.1.2         rmarkdown_2.32       vctrs_0.7.3         
##  [10] memoise_2.0.1        base64enc_0.1-6      rstatix_1.1.0       
##  [13] BiocBaseUtils_1.15.1 htmltools_0.5.9      S4Arrays_1.13.2     
##  [16] forcats_1.0.1        itertools_0.1-3      missForest_1.6.1    
##  [19] curl_8.1.0           broom_1.0.13         SparseArray_1.13.4  
##  [22] Formula_1.2-6        sass_0.4.10          bslib_0.12.0        
##  [25] htmlwidgets_1.6.4    httr2_1.3.0          cachem_1.1.0        
##  [28] buildtools_1.0.0     igraph_2.3.4         mime_0.13           
##  [31] lifecycle_1.0.5      iterators_1.0.14     pkgconfig_2.0.3     
##  [34] Matrix_1.7-6         R6_2.6.1             fastmap_1.2.0       
##  [37] rbibutils_2.4.1      shiny_1.14.0         digest_0.6.39       
##  [40] colorspace_2.1-4     patchwork_1.3.2      Hmisc_5.3-0         
##  [43] RSQLite_3.53.3       ggpubr_1.0.0         filelock_1.0.3      
##  [46] labeling_0.4.3       randomForest_4.7-1.2 polyclip_1.10-7     
##  [49] httr_1.4.9           abind_1.4-8          compiler_4.6.2      
##  [52] rngtools_1.5.2       bit64_4.8.6          withr_3.0.3         
##  [55] htmlTable_2.5.0      S7_0.2.2             backports_1.5.1     
##  [58] viridis_0.6.5        carData_3.0-6        DBI_1.3.0           
##  [61] ggforce_0.5.0        ggsignif_0.6.4       MASS_7.3-66         
##  [64] DelayedArray_0.39.8  ggsci_5.2.0          tools_4.6.2         
##  [67] foreign_0.8-92       ranger_0.18.0        otel_0.2.0          
##  [70] httpuv_1.6.17        visdat_0.6.1         nnet_7.3-21         
##  [73] glue_1.8.1           promises_1.5.0       grid_4.6.2          
##  [76] checkmate_2.3.4      cluster_2.1.8.3      gtable_0.3.6        
##  [79] tzdb_0.5.0           tidyr_1.3.2          data.table_1.18.6.1 
##  [82] hms_1.1.4            tidygraph_1.3.1      utf8_1.2.6          
##  [85] car_3.1-5            XVector_0.53.0       ggrepel_0.9.8       
##  [88] foreach_1.5.2        pillar_1.11.1        stringr_1.6.0       
##  [91] vroom_1.7.1          limma_3.99.4         later_1.4.8         
##  [94] splines_4.6.2        tweenr_2.0.3         dplyr_1.2.1         
##  [97] BiocFileCache_3.3.1  lattice_0.23-1       bit_4.6.0           
## [100] tidyselect_1.2.1     maketools_1.3.2      knitr_1.52          
## [103] gridExtra_2.3.1      edgeR_4.99.12        xfun_0.62           
## [106] graphlayouts_1.2.5   statmod_1.5.2        factoextra_2.2.0    
## [109] DT_0.34.0            stringi_1.8.9        yaml_2.3.12         
## [112] evaluate_1.0.5       codetools_0.2-20     ggraph_2.2.2        
## [115] fenr_1.11.3          tidybulk_2.3.3       tibble_3.3.1        
## [118] BiocManager_1.30.27  cli_3.6.6            rpart_4.1.27        
## [121] Rdpack_2.6.6         xtable_1.8-8         jquerylib_0.1.4     
## [124] Rcpp_1.1.2           dbplyr_2.6.0         parallel_4.6.2      
## [127] ggh4x_0.3.1          ggplot2_4.0.3        readr_2.2.0         
## [130] assertthat_0.2.1     blob_1.3.0           doRNG_1.8.6.3       
## [133] viridisLite_0.4.3    ggridges_0.5.7       scales_1.4.0        
## [136] purrr_1.2.2          crayon_1.5.3         rlang_1.3.0