--- title: > The `EMMA` User's Guide - Enrichment Method MAtters author: - name: Najla Abassi affiliation: - Institute of Medical Biostatistics, Epidemiology and Informatics (IMBEI), Mainz email: abassina@uni-mainz.de - name: Annekathrin Nedwed affiliation: - Institute of Medical Biostatistics, Epidemiology and Informatics (IMBEI), Mainz email: anneludt@uni-mainz.de - name: Federico Marini affiliation: - Institute of Medical Biostatistics, Epidemiology and Informatics (IMBEI), Mainz - Research Center for Immunotherapy (FZI), Mainz email: marinif@uni-mainz.de date: "`r BiocStyle::doc_date()`" package: "`r BiocStyle::pkg_ver('EMMA')`" output: BiocStyle::html_document: toc: true toc_float: true code_folding: show code_download: yes vignette: > %\VignetteIndexEntry{EMMA} %\VignetteEncoding{UTF-8} %\VignettePackage{EMMA} %\VignetteKeywords{Pathways, GeneSetEnrichment, Annotation, ReproducibleResearch, Transcriptomics} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console bibliography: EMMA_bibliography.bib --- ```{r knitr, include = FALSE} knitr::opts_chunk$set( comment = "#>", error = FALSE, warning = FALSE, eval = TRUE, message = FALSE ) ``` # Introduction {#introduction} Functional Enrichment Analysis (FEA) is a key downstream step in omics workflows, commonly applied after differential expression analysis to support biological interpretation and generate pathway-level hypotheses. A wide range of tools and methods are available, mainly Over-Representation Analysis (ORA) [@Khatri2012] and Gene Set Enrichment Analysis (GSEA) [@Subramanian2005], leading to substantial heterogeneity in analytical choices and reported results. Despite its widespread use, FEA is often insufficiently documented [@Wijesooriya2022]. Critical parameters such as background gene sets or multiple testing correction methods are frequently missing or inconsistently reported in scientific papers, limiting reproducibility and interpretability. Currently, no standardized framework exists to ensure transparent and reproducible documentation of FEA workflows, comparable to the MIAME guidelines [@Brazma2001]. Similar approaches have been developed outside the Bioconductor ecosystem to improve transparency and reproducibility in other analytical contexts. For example, `tidylog` records operations performed with `dplyr` and `tidyr`, while [`omicslog`](https://github.com/tidyomics/omicslog) tracks transformations applied to omics-oriented objects. However, no tool specifically addresses the metadata and provenance requirements of FEA workflows. To address this gap, we introduce `r BiocStyle::Biocpkg("EMMA")` (standing for **E**nrichment **M**ethods **MA**tters), a framework that automatically captures key analytical parameters and provenance information during the execution of FEA methods. This vignette demonstrates how `EMMA` integrates with existing tools (e.g. `r BiocStyle::Biocpkg("clusterProfiler")`, `r BiocStyle::Biocpkg("topGO")`, `r BiocStyle::Biocpkg("Enrichr")`, `r BiocStyle::Biocpkg("gprofiler2")`) to execute enrichment analyses while systematically capturing analysis parameters and provenance information during runtime, returning enrichment results **without altering their original format**, alongside structured and reusable metadata. ## What do you get with `EMMA`? Using `EMMA` allows you to: * Record the exact function call and parameters used for FEA * Automatically track annotation sources (e.g. organism, gene set database) * Retain provenance directly within the results object * Generate reproducible summaries of the analysis (e.g. Methods sections) * Facilitate sharing of results together with their analysis context ## How does `EMMA` work? `EMMA` works by wrapping an enrichment call, executing it, and capturing relevant provenance information and parameters during runtime. The recorded metadata is then attached directly to the result object using R’s **attribute system**. This approach enables `EMMA` to preserve provenance information **without modifying the original result structure or introducing new classes**, allowing users to continue working seamlessly with standard outputs from existing tools. While Bioconductor provides dedicated metadata slots for certain S4 classes (e.g. via `metadata()`), these are not consistently available across all enrichment result types. By relying on attributes, provenance information can be attached to any result object regardless of its underlying class. ![](EMMA_workflow.png) # Installation {#installation} To install this package, start R and enter: ```{r install, eval = FALSE} if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager")} BiocManager::install("EMMA") ``` Once installed, the package can be loaded and attached to the current workspace as follows: ```{r setup} library("EMMA") ``` # Usage example of `EMMA`: the `macrophage` dataset ## About the data In the remainder of this vignette, we will illustrate the main features of `r BiocStyle::Biocpkg("EMMA")` on a publicly available dataset from Alasoo, et al. "Shared genetic effects on chromatin and gene expression indicate a role for enhancer priming in immune response", published in Nature Genetics, January 2018 [@Alasoo2018]. The data is made available via the `r BiocStyle::Biocpkg("macrophage")` Bioconductor package, which contains the files output from the Salmon quantification (version 0.12.0, with GENCODE v29 reference), as well as the values summarized at the gene level, which we will use to exemplify. In the `macrophage` experimental setting, the samples are available from 6 different donors, in 4 different conditions (naive, treated with Interferon gamma, with SL1344, or with a combination of Interferon gamma and SL1344). Let's start by loading all the necessary packages: ```{r loadlibraries} library("EMMA") library("macrophage") library("DESeq2") library("org.Hs.eg.db") library("clusterProfiler") library("mosdef") library("topGO") library("GO.db") ``` We will show an example of how `r BiocStyle::Biocpkg("EMMA")` fits into a regular bulk RNA-seq data analysis workflow. ## Getting a list of Differentially Expressed Genes For this, we will load the `macrophage` data and perform Differential Expression Analysis with `r BiocStyle::Biocpkg("DESeq2")` ```{r differential_expression} # load data data(gse, package = "macrophage") # set up design dds_macrophage <- DESeqDataSet(gse, design = ~ line + condition) # preprocess rownames(dds_macrophage) <- substr(rownames(dds_macrophage), 1, 15) keep <- rowSums(counts(dds_macrophage) >= 10) >= 6 dds_macrophage <- dds_macrophage[keep, ] # set seed for reproducibility set.seed(2711) # sample randomly for 2k genes selected_genes <- sample(rownames(dds_macrophage), 2000) dds_macrophage <- dds_macrophage[selected_genes, ] # run DESeq dds_macrophage <- DESeq(dds_macrophage) # get de res for 1st contrast IFNg_vs_naive <- results(dds_macrophage, contrast = c("condition", "IFNg", "naive"), lfcThreshold = 1, alpha = 0.05) IFNg_vs_naive <- lfcShrink(dds_macrophage, coef = "condition_IFNg_vs_naive", res = IFNg_vs_naive, type = "apeglm") IFNg_vs_naive$SYMBOL <- rowData(dds_macrophage)$SYMBOL # get de res for 2st contrast SL1344_vs_naive <- results(dds_macrophage, contrast = c("condition", "SL1344", "naive"), lfcThreshold = 1, alpha = 0.05) SL1344_vs_naive <- lfcShrink(dds_macrophage, coef = "condition_SL1344_vs_naive", res = SL1344_vs_naive, type = "apeglm") SL1344_vs_naive$SYMBOL <- rowData(dds_macrophage)$SYMBOL # sort by adjusted p value de_list <- list( IFNg_vs_naive = IFNg_vs_naive, SL1344_vs_naive = SL1344_vs_naive ) de_list <- lapply(de_list, function(df) { df <- df[order(df$padj), ] df <- df[!is.na(df$padj) & df$padj <= 0.05, ] df }) # set background gene list gene_universe <- rownames(dds_macrophage) ``` ## Performing Functional Enrichment Analysis (FEA) Now that we have a list of DE genes for this contrast, we can perform Functional Enrichment Analysis. We do this first *without* the functionality provided by `EMMA`. In the following example, we will use the function `enrichGO()` from `r BiocStyle::Biocpkg("clusterProfiler")` ```{r withoutemma} de_res <- de_list$IFNg_vs_naive fea_res <- enrichGO(gene = rownames(de_res), universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1) ``` # Entering the `EMMA` framework From now on, we will complement the existing workflow with the functionality provided by `EMMA`. ## `EMMA_run()`: Capturing the recorded information `EMMA_run()` accepts a function call (e.g. `enrichGO(...)`) and executes it as it is, while capturing the associated parameters and provenance information: ```{r EMMA_run_1} # perform FEA, but with EMMA! fea_res <- enrichGO(gene = rownames(de_res), universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1) |> EMMA_run() # simply pipe your call to EMMA_run() # check res fea_res ``` ... or you can simply wrap `EMMA_run()` around your call: ```{r EMMA_run_2} # you can also pass the function name and its namespace # e.g. `clusterProfiler::enrichGO(...)` fea_res_nobg <- EMMA_run( clusterProfiler::enrichGO( # no universe set gene = rownames(de_res), keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE ) ) ``` As you can see, `EMMA` returns the FEA results in their **native/standard** format. `EMMA` also warns you about good practices when performing FEA, like in this example, we didn't define a list of background genes (which can influence the results leading e.g. to a larger amount of false discoveries, as a consequence of over-optimistic smaller p-values [@Wijesooriya2022]), so we get warned about that. ## `EMMA_show()`: Summarizing the recorded information To get a quick summary of what `EMMA` captured while we ran the analysis, we use `EMMA_show()`: ```{r EMMA_show} EMMA_show(fea_res) ``` `EMMA` attaches the captured metadata to the attributes of the results object. That's why it is always a good practice to save the original results, and not only the subsets of interest. ## `EMMA_get_record()`: Retrieving the recorded information To access the full recorded information, we use `EMMA_get_record()`: ```{css} /*| echo: false */ .scrollout { max-height: 500px; overflow-y: scroll; } ``` ```{r EMMA_get_record} #| class-output: scrollout emma_record <- EMMA_get_record(fea_res) # get all the record emma_record ``` `EMMA` structures the `EMMA_record` attribute (i.e. the recorded provenance information) into a list of elements: ```{bash EMMA_record_str, eval=FALSE} ├── EMMA_record │ ├── method # how the analysis was performed │ │ ├── call │ │ ├── function_name │ │ ├── package_name │ │ ├── package_version │ │ ├── wrapped_function │ │ ├── wrapped_package │ │ └── wrapper │ ├── input # inputs used for the analysis │ │ └── arguments │ ├── annotation # annotation context │ │ ├── organism │ │ ├── gene_set_db │ │ └── gene_set_db_version │ ├── timestamp # when the analysis was run │ ├── session_info # R session information │ ├── extra # user-defined additions │ └── emma_version ``` ```{r EMMA_method} # get the method record emma_record$method ``` With `EMMA_run()`, we can decide whether we want to save the value of arguments used in our call or not. This can be useful, for example, to avoid unnecessarily increasing the size of the result object. For this, we can use the argument `args_form`: ```{r argument_form} #| class-output: scrollout fea_res_no_param <- enrichGO(gene = rownames(de_res), universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE) |> EMMA_run(args_form = "unevaluated") # when we don't want the values stored # else set to evaluated (default) # check EMMA_get_record(fea_res_no_param) ``` We can also choose whether to save the R session information with the record using the argument `store_session_info`, which defaults to `TRUE`. ## `EMMA_explain()`: Summarizing recorded information into text `EMMA_explain()` generates a human-readable description of the FEA, similar to a Materials and Methods section of a paper, by summarizing the executed call, the parameters, software context, and reference databases used. ```{r EMMA_explain, results = 'asis', message=TRUE} EMMA_explain(fea_res, get_citation = TRUE) ``` ## `EMMA` with custom/wrapper functions You can also use a custom function that you developed, or a wrapper function (from packages such as `r BiocStyle::Biocpkg("mosdef")`). In this case, `EMMA_run()` will attempt to capture as much metadata as possible: ```{r mosdef_eg} #| class-output: scrollout mosdef_fea_res <- mosdef::run_goseq(de_genes = rownames(de_res), bg_genes = gene_universe, mapping = "org.Hs.eg.db", id = "ensGene", genome = "hg19") |> EMMA_run(store_session_info = FALSE, args_form = "unevaluated") # quick inspection EMMA_get_record(mosdef_fea_res) ``` ```{r custom_eg} #| class-output: scrollout # a custom function (not from a package) my_custom_function <- function(gene, universe = NULL, ontology = "BP", id_type = "ENTREZID", org_db_name = "org.Hs.eg.db", organism = "hsapiens") { res1 <- mosdef::run_topGO(de_genes = gene, bg_genes = gene_universe, ontology = ontology, gene_id = id_type, mapping = org_db_name, add_gene_to_terms = TRUE) res2 <- gprofiler2::gost(query = gene, organism = organism, custom_bg = gene_universe) return(list(topGO_res = res1, gost_res = res2 )) } # run analysis with EMMA frankenstein_fea <- my_custom_function( gene = rownames(de_res), universe = gene_universe, ontology = "BP", id_type = "ENSEMBL", org_db_name = "org.Hs.eg.db", organism = "hsapiens" ) |> EMMA_run(store_session_info = FALSE, args_form = "unevaluated") # quick inspection EMMA_get_record(frankenstein_fea) ``` ## `EMMA_add_custom_metadata()`: Adding extra information The user can always attach extra metadata that `EMMA` might not be able to capture automatically. To keep everything organized, we can use `EMMA_add_custom_metadata()` function ```{r add_custom_metadata} #| class-output: scrollout frankenstein_fea2 <- EMMA_add_custom_metadata(res = frankenstein_fea, extra = list( wrapped_function_topGO = "runTest", notes = "any other meaningful info")) EMMA_get_record(frankenstein_fea2)$extra ``` Since the `EMMA_record` is attached as attribute to the original results objects, it can be preserved when integrating results into structured containers such as `r BiocStyle::Biocpkg("DeeDeeExperiment")` [@Abassi2026]. This enables both FEA results and their associated provenance information to be stored and managed together, facilitating reproducibility, organization, and sharing of complex omics analyses. ```{r emma_and_dde} #| class-output: scrollout dde <- DeeDeeExperiment::DeeDeeExperiment(sce = dds_macrophage, de_results = IFNg_vs_naive, enrich_results = list( IFNg_vs_naive = fea_res_no_param) ) fea <- DeeDeeExperiment::getFEA(dde, format = "original") EMMA_get_record(fea) ``` ## `EMMA_freeze()`: Recording the Analysis Environment `EMMA_freeze()` records the R environment at the time of analysis by generating a lockfile using `renv`. By default, the snapshot is created with `force = TRUE`, allowing the environment to be recorded even when inconsistencies (e.g. version mismatches) are present. This behavior reflects the goal of preserving the analysis environment as it was used in practice, rather than attempting to enforce a fully consistent state. ```{r EMMA_freeze} if (requireNamespace("renv", quietly = TRUE)) { project_path <- tempdir() dir.create(project_path) EMMA_freeze(project = project_path, file = "analysis.lock", pkgs = loadedNamespaces(), prompt = FALSE, force = TRUE) } ``` ## `EMMA` with iterative workflows In real-world analyses, multiple FEA results from different contrasts are often analyzed in parallel This can be conveniently handled using iterative approaches such as `lapply()`. ```{r emma_with_lapply} fea_res_list <- lapply(de_list, function(dea){ #### perform fea with EMMA_run #### enrichGO(gene = rownames(dea), universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE) |> EMMA_run() }) #### print the metadata summary for each contrast #### invisible(lapply(names(fea_res_list), function(nm) { cat("\n###", nm) EMMA_show(fea_res_list[[nm]]) })) #### add more info to the metadata #### fea_res_list <- setNames(lapply(names(fea_res_list), function(nm) { EMMA_add_custom_metadata(fea_res_list[[nm]], extra = list(contrast_name = nm)) }), names(fea_res_list)) #### get the raw metadata #### record_list <- list() record_list <- setNames(lapply(names(fea_res_list), function(mn){ EMMA_get_record(fea_res_list[[mn]]) }), names(fea_res_list)) record_list$SL1344_vs_naive$extra record_list$SL1344_vs_naive$method$package_name ``` # Supported methods and future extensions `EMMA` currently provides metadata capture support for several commonly used FEA tools, including methods from `r BiocStyle::Biocpkg("clusterProfiler")`, `r BiocStyle::Biocpkg("gprofiler2")`, and related wrappers/custom functions such as `r BiocStyle::Biocpkg("mosdef")`. `EMMA` is designed to be extensible, and support for additional FEA methods and packages will continue to expand over time. # Session info {.unnumbered} ```{css} /*| echo: false */ .smaller { font-size: 10px } ``` ```{r sessioinfo} #| class-output: smaller sessionInfo() ``` # References {.unnumbered}