diff --git a/.Rbuildignore b/.Rbuildignore index b7aa75c..2a7eb9b 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -7,3 +7,5 @@ ^_pkgdown\.yml$ ^docs$ ^pkgdown$ +^\.positai$ +^\.claude$ diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 5bef95c..5631f06 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -1,28 +1,146 @@ -# Workflow derived from https://github.com/r-lib/actions/tree/master/examples -# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help on: push: - branches: [main, master] pull_request: - branches: [main, master] + branches: + - devel + schedule: + - cron: '0 8 * * 5' name: R-CMD-check jobs: R-CMD-check: - runs-on: ubuntu-latest + runs-on: ${{ matrix.config.os }} + container: ${{ matrix.config.image }} + + name: ${{ matrix.config.os }} (${{ matrix.config.bioc }} - ${{ matrix.config.image }}) + + strategy: + fail-fast: false + matrix: + config: + - { os: windows-latest, bioc: 'devel'} + - { os: macOS-latest, bioc: 'devel', curlConfigPath: '/usr/bin/'} + - { os: ubuntu-latest, bioc: 'devel'} + # - { os: ubuntu-latest, bioc: 'devel', cran: "https://demo.rstudiopm.com/all/__linux__/xenial/latest"} + - { os: ubuntu-latest, image: 'bioconductor/bioconductor_docker:devel'} + env: + R_REMOTES_NO_ERRORS_FROM_WARNINGS: true + CRAN: ${{ matrix.config.cran }} + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - R_KEEP_PKG_SOURCE: yes + CURL_CONFIG: ${{ matrix.config.curlConfigPath }}curl-config + steps: - - uses: actions/checkout@v2 + - name: Check out repo + uses: actions/checkout@v2 - - uses: r-lib/actions/setup-r@v1 + - name: Set up R and install BiocManager + uses: grimbough/bioc-actions/setup-bioc@v1 + if: matrix.config.image == null with: - use-public-rspm: true + bioc-version: ${{ matrix.config.bioc }} + + - name: Set up pandoc + uses: r-lib/actions/setup-pandoc@v2 + if: matrix.config.image == null + + - name: Install remotes + run: | + install.packages('remotes') + shell: Rscript {0} - - uses: r-lib/actions/setup-r-dependencies@v1 + - name: Query dependencies + run: | + saveRDS(remotes::dev_package_deps(dependencies = TRUE, repos = c(getOption('repos'), BiocManager::repositories())), 'depends.Rds', version = 2) + shell: Rscript {0} + + - name: Cache R packages + if: runner.os != 'Windows' && matrix.config.image == null + uses: actions/cache@v4 + with: + path: ${{ env.R_LIBS_USER }} + key: ${{ runner.os }}-bioc-${{ matrix.config.bioc }}-${{ hashFiles('depends.Rds') }} + restore-keys: ${{ runner.os }}-bioc-${{ matrix.config.bioc }}- + + - name: Install system dependencies + if: runner.os == 'Linux' + env: + RHUB_PLATFORM: linux-x86_64-ubuntu-gcc + uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: rcmdcheck + extra-packages: any::rcmdcheck + pak-version: devel + + - name: Install system dependencies (macOS) + if: runner.os == 'macOS' + run: | + brew install cairo + brew install --cask xquartz + brew install harfbuzz + brew install fribidi + brew install fftw + + # - name: Set up gfortran symlinks (macOS) + # if: runner.os == 'macOS' + # run: | + # set -x + # sudo ln -s /usr/local/Cellar/gcc@11/*/lib/gcc/11 /usr/local/gfortran/lib + # gfortran --version + + - name: Install dependencies + run: | + local_deps <- remotes::local_package_deps(dependencies = TRUE) + deps <- remotes::dev_package_deps(dependencies = TRUE, repos = BiocManager::repositories()) + BiocManager::install(local_deps[local_deps %in% deps$package[deps$diff != 0]], Ncpu = 2L) + remotes::install_cran('rcmdcheck', Ncpu = 2L) + shell: Rscript {0} + + - name: Session info + run: | + options(width = 100) + pkgs <- installed.packages()[, "Package"] + sessioninfo::session_info(pkgs, include_base = TRUE) + shell: Rscript {0} + + - name: Build, Install, Check + uses: grimbough/bioc-actions/build-install-check@v1 + + - name: Run BiocCheck + uses: grimbough/bioc-actions/run-BiocCheck@v1 + with: + arguments: '--no-check-bioc-views --no-check-bioc-help' + error-on: 'error' + + # - name: Upload check results + # if: failure() + # uses: actions/upload-artifact@master + # with: + # name: ${{ runner.os }}-r${{ matrix.config.r }}-results + # path: check + + - name: Show testthat output + if: always() + run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true + shell: bash + + - name: Upload check results + if: failure() + uses: actions/upload-artifact@master + with: + name: ${{ runner.os }}-bioc-${{ matrix.config.bioc }}-results + path: check + + - name: Test coverage + if: matrix.config.os == 'macOS-latest' + run: | + install.packages("covr") + covr::codecov(token = "${{secrets.CODECOV_TOKEN}}") + shell: Rscript {0} - - uses: r-lib/actions/check-r-package@v1 + - name: Deploy + if: github.event_name == 'push' && github.ref == 'refs/heads/devel' && matrix.config.os == 'macOS-latest' + run: | + R CMD INSTALL . + Rscript -e "remotes::install_dev('pkgdown'); pkgdown::deploy_to_branch(new_process = FALSE)" diff --git a/DESCRIPTION b/DESCRIPTION index b44fb89..da41875 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: EMMA -Title: EMMA: Enrichment Methods Matter for enabeling fully reproducible and +Title: EMMA: Enrichment Methods Matter for enabling fully reproducible and provenance-aware pathway analysis -Version: 0.0.1 +Version: 0.99.0 Authors@R: c( person( @@ -28,15 +28,27 @@ Description: EMMA is a package that provides a provenance-aware execution License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.3 +RoxygenNote: 8.0.0 +Imports: + AnnotationDbi, + cli Suggests: knitr, rmarkdown, + macrophage, + DESeq2, + clusterProfiler, + gprofiler2, org.Hs.eg.db, + GO.db, + mosdef, + topGO, + renv, + DeeDeeExperiment, testthat (>= 3.0.0) VignetteBuilder: knitr URL: https://github.com/imbeimainz/EMMA BugReports: https://github.com/imbeimainz/EMMA/issues -biocViews: Pathways, GO, KEGG, GeneSetEnrichment, Annotation, ImmunoOncology, - ReproducibleResearch, Transcriptomics, SingleCell, Software +biocViews: Software, Pathways, GO, KEGG, GeneSetEnrichment, ImmunoOncology, + Transcriptomics, SingleCell, GeneExpression, DifferentialExpression Config/testthat/edition: 3 diff --git a/NAMESPACE b/NAMESPACE index 4ee2689..b47fb89 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,12 @@ # Generated by roxygen2: do not edit by hand +export(EMMA_add_custom_metadata) +export(EMMA_explain) +export(EMMA_freeze) +export(EMMA_get_record) export(EMMA_run) export(EMMA_show) +importFrom(AnnotationDbi,metadata) +importFrom(cli,cli_alert_info) importFrom(utils,packageVersion) +importFrom(utils,sessionInfo) diff --git a/NEWS.md b/NEWS.md index bbd341d..667354d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,39 @@ +# EMMA 0.99.0 + +* Ready for Bioconductor submission! + +# EMMA 0.3.0 + +* Added the initial implementation of `EMMA_freeze()`. +* `EMMA_explain()` now returns also citations of the used packages. +* Renamed `getEMMARecord()` to `EMMA_get_record()` for consistency. +* Resizing toy data to speed up examples. +* Updated the vignette. +* Updating unit tests. + +# EMMA 0.2.0 + +* `EMMA_run()` can accept custom functions and wrappers and collect metadata +based on which function was used in the wrapper. +* Added `EMMA_add_custom_metadata()` to give the user manual/easy access to modify +`user_metdata` field in `EMMA_record`. +* Added a fully runnable vignette. + + +# EMMA 0.1.0 + +* `EMMA_run()` captures a function call, executes the FEA analysis, and returns +the results in their native format while attaching structured metadata to the +result object as attribute. +* `EMMA_run()` can now capture information from functions in `mosdef`, +`gprofiler2`, and 5 commonly used functions from `clusterProfiler`. +* `EMMA_show()` prints the captured metadata in a more user-friendly format. +* Added `getEMMARecord()` to fetch all metadata stored in the attributes of an +object returned by `EMMA_run()`. +* Added toy data. +* Added the initial implementation of `EMMA_explain`. + + # EMMA 0.0.1 * Added the initial implementation of `EMMA_run()` and `EMMA_show()`. diff --git a/R/EMMA-data.R b/R/EMMA-data.R new file mode 100644 index 0000000..d7ad654 --- /dev/null +++ b/R/EMMA-data.R @@ -0,0 +1,70 @@ +#' A sample `data.frame` containing Differential Expression Analysis, generated +#' with `DESeq2` +#' +#' @details This `data.frame` object contains the results of a Differential +#' Expression Analysis performed on data from the `macrophage` package, more +#' precisely contrasting the counts from naive macrophage to those associated +#' with IFNg. +#' +#' The code to create said object can be found in the folder `/inst/scripts` in +#' the EMMA package, the file is called `create_datasets_examples.R`. +#' +#' @return A sample `data.frame` object, extracted from `DESeq2` results +#' +#' @format A `data.frame` object +#' +#' @references Alasoo, et al. "Shared genetic effects on chromatin and gene +#' expression indicate a role for enhancer priming in immune response", +#' Nature Genetics, January 2018 doi: 10.1038/s41588-018-0046-7. +#' +#' @name de_res_IFNg_vs_naive +#' @docType data +NULL + +#' A sample `character vector` containing the background gene list used to +#' perform FEA on the `macrophage` dataset +#' +#' @details This `character vector` object contains the assay's `rownames` +#' of the `macrophage` data +#' +#' The code to create said object can be found in the folder `/inst/scripts` in +#' the EMMA package, the file is called `create_datasets_examples.R`. +#' +#' @return A sample `character vector` containing the assay's `rownames` +#' of the `macrophage` data +#' +#' @format A `character vector` +#' +#' +#' @references Alasoo, et al. "Shared genetic effects on chromatin and gene +#' expression indicate a role for enhancer priming in immune response", +#' Nature Genetics, January 2018 doi: 10.1038/s41588-018-0046-7. +#' +#' @name gene_universe +#' @docType data +NULL + +#' A sample `list` containing Functional Enrichment Analysis results, +#' generated with `gprofiler2` +#' +#' @details This `list` object contains the result table and metadata of the +#' functional enrichment analysis (FEA) performed on the `macrophage` data, +#' specifically using the `gost()` function from the `gprofiler2` package, and +#' wrapped in `EMMA_run()` +#' +#' The code to create said object can be found in the folder `/inst/scripts` in +#' the EMMA package, the file is called `create_datasets_examples.R`. +#' +#' @return A sample `list` containing the FEA results `result` and `metadata`. +#' This results object has the `EMMA_record` attribute. +#' +#' @format A `list` +#' +#' +#' @references Alasoo, et al. "Shared genetic effects on chromatin and gene +#' expression indicate a role for enhancer priming in immune response", +#' Nature Genetics, January 2018 doi: 10.1038/s41588-018-0046-7. +#' +#' @name fea_res +#' @docType data +NULL diff --git a/R/EMMA-pkg.R b/R/EMMA-pkg.R index 328fa03..1537c72 100644 --- a/R/EMMA-pkg.R +++ b/R/EMMA-pkg.R @@ -3,8 +3,9 @@ #' EMMA stands for Enrichment Methods Matter. #' And EMMA stands to help you in realizing it. #' -#' -#' @importFrom utils packageVersion +#' @importFrom AnnotationDbi metadata +#' @importFrom utils packageVersion sessionInfo +#' @importFrom cli cli_alert_info #' @name EMMA-pkg #' @docType package "_PACKAGE" diff --git a/R/EMMA_add_custom_metadata.R b/R/EMMA_add_custom_metadata.R new file mode 100644 index 0000000..675cd0b --- /dev/null +++ b/R/EMMA_add_custom_metadata.R @@ -0,0 +1,39 @@ +#' EMMA_add_custom_metadata +#' +#' Append or replace the `extra` field in the `EMMA_record` attribute +#' of a result object returned by `EMMA_run()`. This allows users to manually +#' provide additional annotation or contextual information that could not be +#' captured automatically +#' +#' @param res A functional enrichment analysis results object as returned by +#' `EMMA_run()` +#' @param extra A named list of user-defined metadata elements to store in +#' the `extra` field +#' +#' @returns The input result object with updated `EMMA_record` attribute +#' @export +#' +#' @examples +#' data("fea_res", package = "EMMA") +#' fea_res <- EMMA_add_custom_metadata(fea_res, extra = +#' list(note = "The background gene set list was all expressed genes in the assay")) +EMMA_add_custom_metadata <- function(res, + extra = list()) { + + if (!is.list(extra)) { + stop("`extra` must be a list!") + } + + if (length(extra) > 0L && is.null(names(extra))) { + stop("`extra` must be a named list!") + } + + emma_rec <- EMMA_get_record(res = res) + + emma_rec$extra <- extra + + # update + attr(res, "EMMA_record") <- emma_rec + + return(res) +} diff --git a/R/EMMA_explain.R b/R/EMMA_explain.R index 3a32aea..c741cfa 100644 --- a/R/EMMA_explain.R +++ b/R/EMMA_explain.R @@ -4,6 +4,124 @@ #' to a Materials and Methods section of a paper, by summarizing the executed #' call, the parameters, software context, and reference databases used. #' -EMMA_explain <- function(res){ - # code -} \ No newline at end of file +#' @param res A functional enrichment analysis results object as returned by +#' `EMMA_run()`. Its attributes contain `EMMA_record`, which contains all +#' provenance information of the performed FEA +#' @param get_citation Logical indicating whether to display the citations of +#' the packages used in the FEA. It only prints the citations in an interactive +#' session (e.g console). Defaults to `TRUE` +#' @returns A character string describing how the FEA was performed using the +#' recorded metadata +#' +#' @export +#' @examples +#' data("fea_res", package = "EMMA") +#' EMMA_explain(fea_res) +#' +EMMA_explain <- function(res, get_citation = TRUE){ + + emma_rec <- EMMA_get_record(res) + + function_name <- emma_rec$method$function_name + pkg_name <- emma_rec$method$package_name + pkg_version <- emma_rec$method$package_version + db <- emma_rec$annotation$gene_set_db + db_version <- emma_rec$annotation$gene_set_db_version + + args <- emma_rec$input$arguments + arg_names <- names(args) + + cli::cli_alert_info( + "You can always complete your text with additional information from `EMMA_get_record()`!" + ) + + + if (isTRUE(emma_rec$method$wrapper)) { + + text <- paste0("Functional Enrichment Analysis was performed using a wrapper function ", + function_name, "()") + } else { + text <- paste0("Functional Enrichment Analysis was performed using the ", + function_name, "() function") + } + + # checks to avoid text with NA + if (!is.null(pkg_name) && !is.na(pkg_name)) { + text <- paste0(text, " from the ", pkg_name, " package") + } + + if (!is.null(pkg_version) && !is.na(pkg_version)) { + text <- paste0(text, " (version ", pkg_version, ")") + } + + if (!is.null(db) && !all(is.na(db))) { + text <- paste0(text, " with the ", + paste(db, collapse = ", "), " database") + } + + if (!is.null(db_version) && !all(is.na(db_version))) { + text <- paste0(text, " (version ", + paste(db_version, collapse = ", "),")") + } + + text <- paste0(text, ".") + + ### info abt bg genes + bg_arg <- intersect(c("universe", "background", "custom_bg", "bg_genes"), + arg_names) + fdr_arg <- intersect(c("correction_method", "pAdjustMethod"), arg_names) + + if (length(bg_arg) == 1) { + bg_value <- args[[bg_arg]] + + if (length(bg_value) > 1) { + # evaluated arguments + text <- paste0(text, + " A custom background gene set was provided (n = ", + length(bg_value),").") + + } else if (length(bg_value) == 1) { + # unevaluated arguments. + text <- paste0(text, " A custom background gene set was provided.") + } + + } else { + text <- paste0(text, " No custom background gene set was recorded.") + } + + + ### info abt the fdr correction + if (length(fdr_arg) == 1) { + fdr_value <- args[[fdr_arg]] + text <- paste0(text, + " Multiple testing correction was performed using the ", + fdr_value, " method." + ) + } else if ("do_padj" %in% arg_names) { + if (isTRUE(args[["do_padj"]])) { + text <- paste0(text, " Multiple testing correction was applied.") + } else if (identical(args[["do_padj"]], FALSE)) { + text <- paste0(text, " Multiple testing correction was not applied.") + } + } + + ### get citations + + if (get_citation) { + pkgs <- unique(stats::na.omit(c( + emma_rec$method$package_name, + emma_rec$method$wrapped_package + ))) + + if (length(pkgs) > 0L) { + cli::cli_alert_info("References:") + for (pkg in pkgs) { + cli::cli_verbatim(paste(format(utils::citation(pkg)), + collapse = "\n")) + } + } + + } + + return(text) +} diff --git a/R/EMMA_freeze.R b/R/EMMA_freeze.R index 0b2efb0..37f6475 100644 --- a/R/EMMA_freeze.R +++ b/R/EMMA_freeze.R @@ -1,8 +1,113 @@ #' EMMA_freeze -#' +#' #' This function records the R environment during analysis runtime and generates -#' a lockfile that can be used with tools such as `renv` +#' a lockfile that can be used with tools such as `renv`. +#' By default, all currently loaded namespaces are recorded. +#' +#' @param project Character string corresponding to the path to the project +#' directory where the lockfile should be written. If the directory does not +#' exist, it will be created. It defaults to the current working directory +#' @param file Character string referring to the name of the lockfile to generate. +#' It defaults to "renv.lock" +#' @param pkgs Character vector of package names to snapshot. It defaults to all +#' currently loaded namespaces via `loadedNamespaces()` +#' @param prompt Logical indicating whether to prompt before taking actions. +#' Defaults to `interactive()` +#' @param force Logical indicating whether to force creation of the lockfile. +#' Defaults to `TRUE` +#' +#' @details +#' This function calls `renv::snapshot()` with the specified packages. +#' The resulting lockfile can later be restored with `renv::restore()` to +#' recreate the same package environment. +#' +#' By default, the lockfile is created with `force = TRUE`, allowing snapshot +#' creation even if inconsistencies are detected in the environment. +#' +#' @returns Invisibly returns the path to the generated lockfile. +#' The lockfile is written in JSON format and can be used with `renv::restore()` +#' to recreate the package environment +#' +#' @seealso \code{\link[renv]{snapshot}}, \code{\link[renv]{restore}} +#' +#' @export +#' +#' @examples +#' # create a lockfile +#' if (requireNamespace("renv", quietly = TRUE)) { +#' tmp <- tempfile("emma_env") +#' dir.create(tmp) +#' +#' EMMA_freeze(project = tmp) +#' +#' # inspect generated files +#' list.files(tmp) #' -EMMA_freeze <- function(){ - #code -} \ No newline at end of file +#' # to restore the environment later +#' # renv::restore(project = tmp) +#' } +EMMA_freeze <- function(project = getwd(), + file = "renv.lock", + pkgs = loadedNamespaces(), + prompt = interactive(), + force = TRUE){ + + # when using EMMA_freeze, we only need renv loaded (without attaching) + # to avoid extra heavy dependency. + if (!requireNamespace("renv", quietly = TRUE)) { + stop( + "The 'renv' package is required for EMMA_freeze(). Please install it." + ) + } + + # checks on args + if (!is.character(project) || length(project) != 1L || is.na(project)) { + stop("`project` must be a single character string.") + } + + if (!is.character(file) || length(file) != 1L || is.na(file)) { + stop("`file` must be a single character string.") + } + + if (!is.character(pkgs)) { + stop("`pkgs` must be a character vector of package names.") + } + + if (!is.logical(prompt) || length(prompt) != 1L || is.na(prompt)) { + stop("`prompt` must be TRUE or FALSE.") + } + + if (!is.logical(force) || length(force) != 1L || is.na(force)) { + stop("`force` must be TRUE or FALSE.") + } + + + # create project dir if needed + if (!dir.exists(project)) { + dir.create(project, recursive = TRUE) + } + + lockfile <- file.path(project, file) + + prompt <- isTRUE(prompt) # if interactive it s TRUE, else takes what the user passes + + # no prompt in non interactive sessions + if (!interactive()) { + prompt <- FALSE + } + + ##### snapshot environment + renv::snapshot( + project = project, + lockfile = lockfile, + packages = pkgs, + prompt = prompt, + force = force + ) + + cli::cli_alert_info( + "Environment snapshot saved to: {.val {lockfile}}. \nTo recreate this environment, use `renv::restore()`.") + + invisible(lockfile) + +} diff --git a/R/EMMA_get_record.R b/R/EMMA_get_record.R new file mode 100644 index 0000000..f105ef2 --- /dev/null +++ b/R/EMMA_get_record.R @@ -0,0 +1,28 @@ +#' EMMA_get_record +#' +#' @param res Functional Enrichment Analysis results (enrichResult, gseaResult ...) +#' generated by `EMMA_run()` +#' +#' @return list of metadata recorded during FEA runtime +#' @export +#' +#' @seealso [enrichGO()], [groupGO()], [gseGO()], [gseKEGG()], [enrichKEGG()], +#' [gost()], [run_cluPro()], [run_topGO()], [run_goseq()] +#' +#' @examples +#' data("fea_res", package = "EMMA") +#' EMMA_get_record(fea_res) +EMMA_get_record <- function(res){ + + rec <- attr(res, "EMMA_record") + + if (is.null(rec)) { + stop("No `EMMA_record` was found. Try running `EMMA_run()` first.") + } + + if (!is.list(rec)) { + stop("Invalid structure. `EMMA_record` must be a `list`") + } + + attr(res, "EMMA_record") +} \ No newline at end of file diff --git a/R/EMMA_internal-utils.R b/R/EMMA_internal-utils.R new file mode 100644 index 0000000..62e8a0e --- /dev/null +++ b/R/EMMA_internal-utils.R @@ -0,0 +1,647 @@ +# metadata capture ------------------------------------------------------------- +#' This function creates a standardized metadata structure used to store +#' annotation info for FEAs. This ensures a consistent structure across +#' different enrichment methods and packages +#' +#' @noRd +.EMMA_empty_metadata <- function() { + list( + organism = NA_character_, + gene_set_db = NA_character_, + gene_set_db_version = NA_character_ +) +} + + +#' This function acts as a dispatcher, routing to package-specific helpers +#' (e.g. for `clusterProfiler`, `gprofiler2` ...) to standardize metadata +#' collection across different enrichment tools +#' +#' +#' @param call_class A named list returned by `.EMMA_classify_call()`. It must +#' contain an `info_call` list with at least `function_name`, `package_name`, +#' and the original `call` +#' @param args A list of evaluated arguments passed to the enrichment function +#' @param envir The environment in which the original call was evaluated +#' +#' @return A list containing annotation metadata(organism, gene set database and +#' its version), depending on the originating package. It returns an empty +#' metadata structure if no package-specific method is available +#' +#' @noRd +.EMMA_get_metadata <- function(call_class, + args, + envir = parent.frame()) { + info_call <- call_class$info_call + function_name <- info_call$function_name + package_name <- info_call$package_name + call <- info_call$call + + meta <- .EMMA_empty_metadata() + + if (is.null(package_name) || is.na(package_name) || package_name == "") { + package_name <- "custom" + } + + meta <- switch( + package_name, + clusterProfiler = .EMMA_get_clusterprofiler_metadata(function_name, args), + gprofiler2 = .EMMA_get_gprofiler2_metadata(args), + mosdef = .EMMA_cp_GO_metadata(args$mapping), + custom = .EMMA_get_custom_metadata(call_class, + args, + envir = parent.frame()) + ) + + return(meta) + +} + +#' This function retrieves metadata associated with FEAs performed using +#' functions from the `clusterProfiler` package. The specific metadata +#' extraction method is selected based on the enrichment function used +#' +#' @param function_name A character string specifying the name of the +#' `clusterProfiler` function used to perform FEA +#' @param args A list containing the evaluated arguments passed into the +#' function call to perform FEA +#' +#' @return A list containing annotation metadata (organism, gene set database +#' and its version), depending on the underlying database used (e.g. GO vs KEGG) +#' +#' @noRd +.EMMA_get_clusterprofiler_metadata <- function(function_name, args) { + switch( + function_name, + enrichGO = .EMMA_cp_GO_metadata(args$OrgDb), + gseGO = .EMMA_cp_GO_metadata(args$OrgDb), + groupGO = .EMMA_cp_GO_metadata(args$OrgDb), + enrichKEGG = .EMMA_cp_KEGG_metadata(args), + gseKEGG = .EMMA_cp_KEGG_metadata(args), + .EMMA_empty_metadata() + # more to be added + ) +} + +#' This function extracts the organism from OrgDb objects +#' +#' @param orgdb An `OrgDb` object or a single character string giving the name +#' of such an object (e.g. `"org.Hs.eg.db"`) +#' +#' @return A character string containing the organism name (e.g. `"Homo sapiens"`) +#' and returns `NA_character_` if the organism cannot be determined +#' +#' @noRd +.EMMA_get_organism_from_OrgDb <- function(orgdb) { + + if (is.null(orgdb)) { + return(NA_character_) + } + + if (is.character(orgdb) && length(orgdb) == 1) { + if (!exists(orgdb, mode = "S4")) { + stop(orgdb, + " was provided as a string, but no loaded object with that name was found. ", + "Please load the corresponding org.*.db library") + } + orgdb <- get(orgdb) + } + + md <- AnnotationDbi::metadata(orgdb) + organism <- md$value[md$name == "ORGANISM"] + + if (length(organism) != 1L) { + return(NA_character_) + } + return(organism) +} + + +#' This function assembles standardized metadata describing GO enrichment +#' analyses performed with `clusterProfiler`. The metadata includes the +#' organism, the geneset database used, and its version when available +#' +#' @param org An `OrgDb` object or a single character string giving the name of +#' such an object (e.g. `"org.Hs.eg.db"`) +#' +#' @return A named list containing annotation metadata, mainly the organism, the +#' geneset database and its version +#' +#' @noRd +.EMMA_cp_GO_metadata <- function(org) { + meta <- .EMMA_empty_metadata() + + meta$organism <- .EMMA_get_organism_from_OrgDb(org) + meta$gene_set_db <- "GO" + meta$gene_set_db_version <- if (requireNamespace("GO.db", quietly = TRUE)) { + as.character(utils::packageVersion("GO.db")) + } else { + NA_character_ + } + + return(meta) +} + +#' This function assembles standardized metadata describing KEGG enrichment +#' analyses performed with `clusterProfiler`. The metadata includes the +#' organism, the geneset database used. The database version is currently not +#' captured +#' +#' @param args A list containing the evaluated arguments passed into the +#' function call to perform FEA +#' +#' @return A named list containing annotation metadata, mainly the organism, the +#' geneset database and its version +#' +#' @noRd +.EMMA_cp_KEGG_metadata <- function(args) { + meta <- .EMMA_empty_metadata() + + meta$organism <- args$organism + meta$gene_set_db <- "KEGG" + meta$gene_set_db_version <- NA_character_ + + return(meta) +} + + +#' This function retrieves metadata associated with FEAs performed using +#' functions from the `gprofiler2` package +#' +#' @param args A list containing the evaluated arguments passed into the +#' function call to perform FEA +#' +#' @return A named list containing annotation metadata, mainly the organism, the +#' geneset database and its version +#' +#' @noRd +.EMMA_get_gprofiler2_metadata <- function(args) { + meta <- .EMMA_empty_metadata() + + version_info <- gprofiler2::get_version_info() + + meta$organism <- if (!is.null(args$organism)) args$organism else NA_character_ + + sources_used <- if (!is.null(args$sources)) { + args$sources + } else { + names(version_info[["sources"]]) + } + + # in case the source is GO, we extract all the dbs ... + if ("GO" %in% sources_used) { + sources_used <- unique(c( + setdiff(sources_used, "GO"), + "GO:BP", "GO:CC", "GO:MF" + )) + } + + meta$gene_set_db <- sources_used + + meta$gene_set_db_version <- vapply( + version_info[["sources"]][sources_used], function(x) x[["version"]], + character(1) + ) # needs more work when only GO, no need to print 3 times the same version + # also, differentiate between which version corresponds to which db + + return(meta) +} + + +#' This function retrieves metadata associated with FEAs performed using +#' user-defined wrapper functions or fully custom functions +#' +#' @param call_class A named list returned by `.EMMA_classify_call()` +#' @param args A list of evaluated arguments passed to the enrichment function +#' @param envir The environment in which the original call was evaluated +#' +#' @return A named list +#' +#' @noRd +.EMMA_get_custom_metadata <- function(call_class, + args, + envir = parent.frame()) { + + meta <- .EMMA_empty_metadata() + + if (call_class$type == "wrapper") { + # get all known functions used in the wrapper + wrapped_fun <- unique(call_class$wrapped_original$fun) + + go_funs <- c("enrichGO", "gseGO", "groupGO") + kegg_funs <- c("enrichKEGG", "gseKEGG") + # dont use .EMMA_get_clusterprofiler_metadata cause the switch wont work with + # more than 1 wrapped fun + meta_list <- list() + if (any(wrapped_fun %in% go_funs)) { + meta_list$GO <- .EMMA_cp_GO_metadata(args$OrgDb) + } + + if (any(wrapped_fun %in% kegg_funs)) { + meta_list$KEGG <- .EMMA_cp_KEGG_metadata(args) + } + + if (any(wrapped_fun %in% c("gost"))) { + meta_list$gprofiler2 <- .EMMA_get_gprofiler2_metadata(args) + } + + if (length(meta_list) == 0L) { + return(meta) + } + + meta$organism <- unique(unlist(lapply(meta_list, `[[`, "organism"))) + meta$gene_set_db <- unique(unlist(lapply(meta_list, `[[`, "gene_set_db"))) + meta$gene_set_db_version <- unique(unlist(lapply(meta_list, `[[`, "gene_set_db_version"))) + + return(meta) + + } + + if (call_class$type == "custom") { + # if 100% custom + cli::cli_alert_info( + "You used a custom function, so `EMMA` wasn't able to record annotation-related information. Please consider adding the `organism`, `geneset database` and its `version` manually into the `extra` field in EMMA_record using EMMA_add_custom_metadata.") + + return(meta) + } + + return(meta) +} + + + +# call info capture ------------------------------------------------------------ + +#' This function extracts call related metadata. It handles two call +#' formats: bare function calls (e.g. `fun(...)`) and namespace-qualified +#' calls (e.g. `pkg::fun(...)`) +#' +#' @param call A call object +#' @param envir The environment in which to look up the function when a bare +#' call is used +#' +#' @return A named list +#' @noRd +.EMMA_capture_call_info <- function(call, envir = parent.frame()) { + # param checks + if (!is.call(call)) { + stop("`call` must be a function call") + } + # capture function name + call_name <- call[[1]] + # when we only use function name e.g. enrichGO(...) + if (is.symbol(call_name)) { + function_name <- as.character(call_name) + + fun <- get(function_name, envir = envir, mode = "function") + + pkg <- utils::packageName(environment(fun)) + + package_name <- if (is.null(pkg) || pkg == "" ) NA_character_ else pkg + + pkg_version <- if (!is.na(package_name)) { + as.character(packageVersion(package_name))} + + else NA_character_ + + # capture args (unevaluated) + arg_list <- as.list(call)[-1] + + return(list( + call = call, + function_name = function_name, + package_name = package_name, + package_version = pkg_version, + arg_list = arg_list + )) + } + + # when we use function name with namespace e.g. clusterProfiler::enrichGO(...) + if (is.call(call_name) && + length(call_name) == 3L && + identical(call_name[[1]], as.symbol("::"))) { + package_name <- as.character(call_name[[2]]) + function_name <- as.character(call_name[[3]]) + pkg_version <- if (!is.na(package_name)) { + as.character(utils::packageVersion(package_name)) + } else { + NA_character_ + } + # capture args (unevaluated) + arg_list <- as.list(call)[-1] + + return(list( + call = call, + function_name = function_name, + package_name = package_name, + package_version = pkg_version, + arg_list = arg_list + )) + } + + stop( + "Unsupported call format. Use a direct function call like `fun(...)` or `pkg::fun(...)`") +} + + +#' This function, used in in `.EMMA_walk()`, converts the head of a call into +#' a character string representing the function being called. It supports both +#' bare calls (e.g. `fun`) and namespace-qualified calls (e.g. `pkg::fun`) +#' +#' @param x The head of a call (the function being called) +#' +#' @return A character string representing the function name +#' +#' @noRd +.EMMA_call_name <- function(x) { + if (is.symbol(x)) { + return(as.character(x)) + } + if ( + is.call(x) && + identical(x[[1]], as.name("::")) && + length(x) == 3 + ) { + return(paste0(as.character(x[[2]]), "::", as.character(x[[3]]))) + } + NULL +} + + +#' This function recursively collects function names/operators from a call. It +#' is used to inspect the body of a wrapper function and to identify whether it +#' calls any known functional enrichment functions (e.g. `enrichGO()`) +#' +#' @param x A call object +#' +#' @return A character vector containing the names of all functions and +#' operators found in `x`, including nested calls +#' +#' @noRd +.EMMA_walk <- function(x) { + out <- character() + + if (is.call(x)) { + # extract function/operator name + nm <- .EMMA_call_name(x[[1]]) + if (!is.null(nm)) { + out <- c(out, nm) + } + + } + if (is.call(x) || is.pairlist(x) || is.expression(x)) { + #check children/nested elements + for (i in seq_along(x)) { + out <- c(out, .EMMA_walk(x[[i]])) + } + } + + return(out) +} + + +#' This function inspects the body of a function call and determines whether +#' it wraps one or more known FEA functions (e.g. from `clusterProfiler` or +#' any other package). This is used internally by `EMMA_run()` when the top-level +#' call is not itself a known enrichment function, but may be a user-defined +#' wrapper around one +#' +#' @param call An unevaluated call +#' @param envir An environment in which to evaluate `call` +#' +#' @return `NULL` if no known enrichment function is detected in the wrapper body. +#' Otherwise, a list with two components: `pkg`, a character vector of package +#' names, and `fun`, character vector of corresponding function names used in the +#' wrapper +#' +#' @noRd +.EMMA_find_original_wrapped_fun <- function(call, envir = parent.frame()) { + + if (!is.call(call)) { + # do i want it to fail here ? + return(NULL) + } + + available <- list( + clusterProfiler = c("enrichGO", "gseGO", "groupGO", "enrichKEGG", "gseKEGG"), + gprofiler2 = c("gost"), + goseq = c("goseq"), + mosdef = c("run_topGO", "run_cluPro", "run_goseq"), + topGO = c("runTest") + # will add more later here + ) + + # flatten available into full names + full_targets <- mapply(function(pkg, funs){ + paste0(pkg, "::", funs)}, + names(available), #pkg + available, #funs + SIMPLIFY = FALSE) |> unlist(use.names = FALSE) + + # get wrapper name + call_head <- call[[1]] + + # get function object, from which we'll extract the body + # if we have a bare wrapper name + if (is.symbol(call_head)) { + fun <- get(as.character(call_head), envir = envir, mode = "function") + } + # if the wrapper is namespace qualified function name, + # extract the fun name to get the body + else if (is.call(call_head) && + length(call_head) == 3L && + identical(call_head[[1]], as.symbol("::"))) { + pkg <- as.character(call_head[[2]]) + fn <- as.character(call_head[[3]]) + fun <- getExportedValue(pkg, fn) + } else { + stop( + "Only `fun(...)` and `pkg::fun(...)` are supported" + ) + } + + + found <- .EMMA_walk(body(fun)) |> unique() + + matched <- unique(c( + intersect(found, full_targets), + #remove the pkg name (what s before ::) + full_targets[sub("^.*::", "", full_targets) %in% found] + )) + + if (length(matched) == 0) { + return(NULL) + } + + parts <- strsplit(matched, "::", fixed = TRUE) + + return(list( + pkg = vapply(parts, `[`, character(1), 1), + fun = vapply(parts, `[`, character(1), 2) + )) +} + + +#' decide if top-level function passed to EMMA_run is a known fun or a wrapper +#' +#' @param call A call object passed to `EMMA_run()` +#' @param envir The environment in which `call` should be evaluated +#' +#' @return A named list +#' +#' @noRd +.EMMA_classify_call <- function(call, envir = parent.frame()) { + # listing all the functions that are not wrappers + originals <- list( + clusterProfiler = c("enrichGO", "gseGO", "groupGO", "enrichKEGG", "gseKEGG"), + gprofiler2 = c("gost"), + goseq = c("goseq"), + topGO = c("runTest") + ) + + full_originals <- mapply( function(pkg, funs) paste0(pkg, "::", funs), + names(originals), + originals, + SIMPLIFY = FALSE) |> unlist(use.names = FALSE) + + # capture call information + info_call <- .EMMA_capture_call_info(call = call, envir = envir) + function_name <- info_call$function_name + package_name <- info_call$package_name + + top_name <- if (!is.na(package_name) && !is.null(package_name) + && package_name != "" && length(package_name) == 1L) { + paste0(package_name, "::", function_name) + } else { + function_name + } + + is_original_top_level <- if (!is.na(package_name) && !is.null(package_name) + && package_name != "" && length(package_name) == 1L) { + paste0(package_name, "::", function_name) %in% full_originals + } else { + function_name %in% sub("^.*::", "", full_originals) + } + + wrapped_original <- NULL + wrapper <- FALSE # default + + if (!is_original_top_level) { + # inspect body of expr function, because this could be a wrapper + wrapped_original <- .EMMA_find_original_wrapped_fun(call = call, + envir = envir) + wrapper <- !is.null(wrapped_original) + } + + type <- if (is_original_top_level) { + "original" + } else if (wrapper) { + "wrapper" + } else { + "custom" + } + + return(list( + type = type, + is_original_top_level = is_original_top_level, + wrapper = wrapper, + wrapped_original = wrapped_original, + info_call = info_call + )) + +} + + +# build EMMA record ------------------------------------------------------------ + + +#' This function assembles the structured provenance record that is stored as +#' an attribute on the FEA results object +#' +#' @param info_call A list returned by `.EMMA_capture_call_info()` +#' @param args_form A character string, either `"evaluated"` or `"unevaluated"` +#' to decide how to store the arguments +#' @param metadata A list returned by `.EMMA_get_metadata()` +#' @param start_time A timestamp marking when the enrichment analysis started +#' @param store_session_info Logical. If `TRUE`, `sessionInfo()` is captured +#' and stored in the record; if `FALSE` the `session_info` slot is `NULL` +#' +#' @return A named list of the recorded metadata +#' +#' @noRd +.EMMA_build_record <- function(info_call, args_form, metadata, + wrapped_original, wrapper, + start_time, store_session_info) { + emma_rec <- list( + method = list( + call = info_call$call, + function_name = info_call$function_name, + package_name = info_call$package_name, + package_version = info_call$package_version, + wrapped_function = if (!is.null(wrapped_original)) + wrapped_original$fun else NULL, + wrapped_package = if (!is.null(wrapped_original)) + wrapped_original$pkg else NULL, + wrapper = wrapper + ), + input = list( + arguments = if (args_form == "evaluated") info_call$args + else info_call$arg_list + ), + annotation = list( + organism = metadata$organism, + gene_set_db = metadata$gene_set_db, + gene_set_db_version = metadata$gene_set_db_version + ), + timestamp = start_time, + session_info = if (isTRUE(store_session_info)) sessionInfo() + else NULL, + extra = list(),# free field for extra metadata added by user + emma_version = as.character(packageVersion(pkg = "EMMA")) + ) + + return(emma_rec) + +} + +# good practice warnings ------------------------------------------------------- + + +#' This function warns about missing good-practice arguments in enrichment calls. +#' A warning is raised if none of the synonyms for a given category appear +#' in `arg_names` +#' +#' @param arg_names A character vector of argument names as written in the +#' user's call, obtained within `EMMA_run()` +#' @param function_name A character corresponding to the name of the enrichment +#' function called, used only for constructing the warning message +#' +#' @return `base::invisible()` +#' @noRd +.EMMA_warnings <- function(arg_names, function_name){ + checks <- list( + list( + params = c("pAdjustMethod", "correction_method", "do_padj"), + label = "multiple-testing correction method" + ), + list( + params = c("universe", "background", "custom_bg", "bg_genes"), + label = "background gene set" + ) + ) + + for (check in checks) { + if (!any(check$params %in% arg_names)) { + warning( + sprintf( + "No %s was specified for %s(). +Consider using the corresponding parameter for your call.", + check$label, + function_name + ), + call. = FALSE + ) + } + } + invisible() +} diff --git a/R/EMMA_run.R b/R/EMMA_run.R index 05a978c..5926519 100644 --- a/R/EMMA_run.R +++ b/R/EMMA_run.R @@ -1,20 +1,48 @@ #' EMMA_run #' -#' This function executes functional enrichment analysis using existing tools -#' and captures the associated parameters and provenance information for the -#' analysis during runtime. +#' This function executes any functional enrichment analysis function and +#' attaches a provenance record (`EMMA_record`) describing the analysis. The +#' captured record includes the original call, metadata derived from the call +#' and its arguments, runtime information, and optionally the current session +#' information. +#' +#' `EMMA_run()` accepts both direct calls to known enrichment functions +#' (`enrichGO()`, `GSEA()`, `fgsea()` ...) and calls to wrapper functions that +#' internally invoke a know enrichment function. +#' +#' #' @param expr A function call that performs functional enrichment analysis. #' The call is captured and executed by EMMA to record analysis parameters and -#' provenance information +#' provenance information. Both bare calls (`enrichGO(...)`) and +#' namespace-qualified calls (`clusterProfiler::enrichGO(...)`) are supported. #' @param envir An environment in which to evaluate `expr` +#' @param store_session_info Logical, indicating whether to store the output of +#' `sessionInfo()` or not. If `TRUE` (default), the session is stored in the +#' provenance record +#' @param args_form A character string indicating whether to store the evaluated +#' or the unevaluated arguments in the provenance record. It default to +#' `"evaluated"` #' -#' @returns Functional enrichment analysis results in the native format -#' returned by the original `expr` +#' @returns The result object returned by the enrichment function in `expr`, +#' in standard format, with an additional `EMMA_record` attribute containing the +#' provenance information. Use `EMMA_get_record()` to retrieve this record +#' #' @export #' #' @examples -#' EMMA_run(enrichGO(gene = rownames(de_res_IFNg_vs_naive), universe = universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP")) -EMMA_run <- function(expr, envir = parent.frame()) { +#' data("de_res_IFNg_vs_naive", package = "EMMA") +#' data("gene_universe", package = "EMMA") +#' library(gprofiler2) +#' +#' EMMA_run(gost(query = de_res_IFNg_vs_naive$SYMBOL, organism = "hsapiens", +#' correction_method = "fdr", custom_bg = gene_universe, sources = "GO:BP"), +#' store_session_info = FALSE, args_form = "unevaluated") +EMMA_run <- function(expr, + envir = parent.frame(), + store_session_info = TRUE, + args_form = c("evaluated", "unevaluated")) { + + args_form <- match.arg(args_form) # capture call call <- substitute(expr) @@ -23,77 +51,62 @@ EMMA_run <- function(expr, envir = parent.frame()) { if (!is.call(call)) { stop("`expr` must be a function call, e.g. enrichGO(...)!") } - - # capture function name - function_name <- paste(deparse(call[[1]]), collapse = "") - - # capture args (unevaluated) - arg_list <- as.list(call)[-1] - - arg_names <- names(arg_list) - # now we'll put some warning on important args if - # the user did not define them (like the fdr correction ...) + if (!is.environment(envir)) { + stop("`envir` must be an environment!") + } - checks <- list( - list(params = c("pAdjustMethod", "correction_method"), - label = "multiple-testing correction method"), - list(params = c("universe", "background"), - label = "background gene set") - ) + if (!is.logical(store_session_info) || length(store_session_info) != 1L || is.na(store_session_info)) { + stop("`store_session_info` must be a single TRUE or FALSE value") + } + # decide whether we are dealing with a known function or a wrapper + call_class <- .EMMA_classify_call(call, envir = envir) + info_call <- call_class$info_call + function_name <- info_call$function_name + package_name <- info_call$package_name + wrapper <- call_class$wrapper + wrapped_original <- call_class$wrapped_original - for (check in checks) { - if (!any(check$params %in% arg_names)) { - warning( - sprintf( - "No %s was specified for %s(). -Consider using the corresponding parameter for your call.", - check$label, - function_name - ), - call. = FALSE - ) - } - } + # capture args (unevaluated) + arg_list <- info_call$arg_list + arg_names <- names(arg_list) - # maybe warn when the organism and the gene names dont match? like u have - # mouse data but u defined organism as human + # some good practice warning, i.e. when multiple testing correction is skipped + # or bg geneset not set + .EMMA_warnings(arg_names = arg_names, + function_name = function_name) #capture analysis time start_time <- Sys.time() - message("Running Enrichment Analysis with ", function_name, " ...") # maybe we can print the list of arg used also in the msg? + + cli::cli_alert_info( + "Running Enrichment Analysis with {.val {function_name}} ...") # capture the value of the arguments args <- lapply(arg_list, eval, envir = envir) + info_call$args <- args # get the function fun <- eval(call[[1]], envir = envir) - # run the analysis + ##### run the analysis #### results <- do.call(fun, args) - # capture pkg used - enrich_function <- match.fun(function_name) - pkg <- environmentName(environment(enrich_function)) - - - # store everything - EMMA_record <- list( - call = call, - package = pkg, - package_version = if (!is.na(pkg)) - as.character(packageVersion(pkg)) else NA, - arguments = args, - gene_set_library = NA, #placeholder - gene_set_library_version = NA, - organism = NA, - runtime = start_time, - session = sessionInfo() + # capture metadata from the used function and arguments + metadata <- .EMMA_get_metadata( + call_class, + args, + envir = envir ) - # store the emma records as attribute of the results obj + # record everything in EMMA_record + EMMA_record <- .EMMA_build_record(info_call, args_form, metadata, + wrapped_original,wrapper, + start_time, store_session_info) + + # attach EMMA_record as attribute of the results obj attr(results, "EMMA_record") <- EMMA_record - + return(results) } diff --git a/R/EMMA_show.R b/R/EMMA_show.R index b8ad0a2..58c5657 100644 --- a/R/EMMA_show.R +++ b/R/EMMA_show.R @@ -1,23 +1,62 @@ #' EMMA_show #' -#' This function prints the EMMA record associated with a functional enrichment -#' analysis object, including the executed call, parameters, and provenance -#' information. +#' This function displays a human-readable summary of the `EMMA_record` attribute +#' attached to a result object produced by `EMMA_run()` +#' #' @param res A functional enrichment analysis results object as returned by #' `EMMA_run()` #' -#' @returns A list of the recorded information +#' @returns `base::invisible()` #' @export #' #' @examples -#' res <- EMMA_run(enrichGO(gene = rownames(de_res_IFNg_vs_naive), universe = universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP")) -#' EMMA_show(res) +#' data("fea_res", package = "EMMA") +#' EMMA_show(fea_res) EMMA_show <- function(res){ + if ("EMMA_record" %in% names(attributes(res))) { - message("Found EMMA record!!") - emma_rec <- attr(res, "EMMA_record") - # should think of a better way to show the records - str(emma_rec) + cli::cli_alert_info("Found EMMA record!") + + emma_rec <- EMMA_get_record(res) + + if (is.list(res) && !is.data.frame(res)) { + # e.g. case of gost, returns a list but it's 1 FEA (result) + if ("result" %in% names(res)) { + cat("Number of Pathways:", NROW(res$result), "\n") + } else { + # let's say if we have of list of FEAs (returned by custom function) + cat("Number of FEAs:", length(res), "\n") + + nms <- names(res) + if (is.null(nms) || any(nms == "")) { + nms <- paste0("FEA_", seq_along(res)) + } + + for (i in seq_along(res)) { + # check the number of pathways for each element of the list + cat(" -", nms[i], ":", NROW(res[[i]]), " pathways\n") + } + } + } else { + cat("Number of Pathways:", NROW(res), "\n") + } + + method_info <- emma_rec$method + db_info <- emma_rec$annotation + + cat("Call:", paste(deparse(method_info$call), collapse = " "), " \n") + cat("Wrapper:", method_info$wrapper, " \n") + cat("Package:", paste0(method_info$package_name , " v. ", + method_info$package_version), " \n") + cat("Organism:", db_info$organism, " \n") + cat("Gene set library:", paste(db_info$gene_set_db, collapse = ","), " \n") + cat("Gene set library version:",db_info$gene_set_db_version, " \n") + cat("\n") + + } else { + warning("No `EMMA_record` attribute was found") } -} \ No newline at end of file + invisible() + +} diff --git a/README.md b/README.md index 55e4485..685fc19 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,10 @@ -# EMMA +# EMMA - Enrichment Methods MAtter -EMMA aims to explicitly capture analytical parameters during functional -enrichment analysis runtime, while returning native enrichment results together -with structured metadata. +EMMA enables the execution of Functional Enrichment Analyses using a wide range +of existing tools (e.g. `clusterProfiler`, `topGO`, `gprofiler2` among others) +while systematically capturing analysis parameters and provenance information +during runtime, and returning enrichment results in their standard format +alongside structured and reusable metadata. ## Installation @@ -21,11 +23,11 @@ remotes::install_github("imbeimainz/EMMA", library(EMMA) # load data data("de_res_IFNg_vs_naive", package = "EMMA") -data("universe", package = "EMMA") +data("gene_universe", package = "EMMA") # run analysis fea_results <- enrichGO(gene = rownames(de_res_IFNg_vs_naive), - universe = universe, + universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP") |> @@ -33,6 +35,11 @@ fea_results <- enrichGO(gene = rownames(de_res_IFNg_vs_naive), ``` +## Usage Overview + +You can find the rendered version of the documentation of `EMMA` at the project +website https://imbeimainz.github.io/EMMA/ + ## Development If you encounter a bug, have usage questions, or want to share ideas and @@ -42,7 +49,7 @@ functionality to make this package better, feel free to file an ## Code of Conduct Please note that the EMMA project is released with a [Contributor Code -of Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html). +of Conduct](https://contributor-covenant.org/version/3/0/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms. ## License diff --git a/_pkgdown.yml b/_pkgdown.yml new file mode 100644 index 0000000..f753f9e --- /dev/null +++ b/_pkgdown.yml @@ -0,0 +1,21 @@ +authors: + Najla Abassi: + href: https://najlaabassi.github.io/ + Annekathrin Nedwed: + href: https://github.com/AnnekathrinSilvia + Federico Marini: + href: https://federicomarini.github.io + +url: ~ +template: + bootstrap: 5 + bslib: + preset: "bootstrap" + font_scale: 1.0 + base_font: + google: "Atkinson Hyperlegible" + code_font: + google: "IBM Plex Mono" + +navbar: + type: light diff --git a/data/de_res_IFNg_vs_naive.RData b/data/de_res_IFNg_vs_naive.RData index 8c20de5..1528601 100644 Binary files a/data/de_res_IFNg_vs_naive.RData and b/data/de_res_IFNg_vs_naive.RData differ diff --git a/data/fea_res.RData b/data/fea_res.RData new file mode 100644 index 0000000..4c5a56a Binary files /dev/null and b/data/fea_res.RData differ diff --git a/data/gene_universe.RData b/data/gene_universe.RData new file mode 100644 index 0000000..4a0f021 Binary files /dev/null and b/data/gene_universe.RData differ diff --git a/data/universe.RData b/data/universe.RData deleted file mode 100644 index 74e9c0d..0000000 Binary files a/data/universe.RData and /dev/null differ diff --git a/inst/scripts/create_datasets_examples.R b/inst/scripts/create_datasets_examples.R index 1603ca9..d388a9c 100644 --- a/inst/scripts/create_datasets_examples.R +++ b/inst/scripts/create_datasets_examples.R @@ -7,15 +7,10 @@ rownames(dds_macrophage) <- substr(rownames(dds_macrophage), 1, 15) dds_macrophage # DE run -keep <- rowSums(counts(dds_macrophage) >= 10) >= 6 -dds_macrophage <- dds_macrophage[keep, ] -dds_macrophage - - # set seed for reproducibility -set.seed(42) -# sample randomly for 1k genes -selected_genes <- sample(rownames(dds_macrophage), 2000) +set.seed(2711) +# sample randomly for 500 genes +selected_genes <- sample(rownames(dds_macrophage), 500) dds_macrophage <- dds_macrophage[selected_genes, ] @@ -38,10 +33,21 @@ de_res_IFNg_vs_naive <- de_res_IFNg_vs_naive[order(de_res_IFNg_vs_naive$padj), ] de_res_IFNg_vs_naive <- de_res_IFNg_vs_naive[!(is.na(de_res_IFNg_vs_naive$padj)) & de_res_IFNg_vs_naive$padj <= 0.05, ] -universe <- rownames(dds_macrophage) +# define gene universe +gene_universe <- rownames(dds_macrophage) -save(de_res_IFNg_vs_naive, file = "de_res_IFNg_vs_naive.RData", compress = "xz") -save(universe, file = "universe.RData", compress = "xz") +# perform FEA +library("gprofiler2") +fea_res <- gprofiler2::gost(query = de_res_IFNg_vs_naive$SYMBOL, + organism = "hsapiens", + correction_method = "fdr", + custom_bg = gene_universe, + sources = "GO:BP") |> EMMA_run( + store_session_info = FALSE, + args_form = "unevaluated") +save(de_res_IFNg_vs_naive, file = "de_res_IFNg_vs_naive.RData", compress = "xz") +save(gene_universe, file = "gene_universe.RData", compress = "xz") +save(fea_res, file = "fea_res.RData", compress = "xz") diff --git a/man/EMMA-pkg.Rd b/man/EMMA-pkg.Rd index c4f1850..e3590ec 100644 --- a/man/EMMA-pkg.Rd +++ b/man/EMMA-pkg.Rd @@ -23,6 +23,7 @@ Useful links: Authors: \itemize{ + \item Najla Abassi \email{abassi.nejla96@gmail.com} (\href{https://orcid.org/0000-0001-8357-0938}{ORCID}) \item Annekathrin Nedwed \email{anneludt@uni-mainz.de} (\href{https://orcid.org/0000-0002-2475-4945}{ORCID}) \item Federico Marini \email{marinif@uni-mainz.de} (\href{https://orcid.org/0000-0003-3252-7758}{ORCID}) } diff --git a/man/EMMA_add_custom_metadata.Rd b/man/EMMA_add_custom_metadata.Rd new file mode 100644 index 0000000..782521b --- /dev/null +++ b/man/EMMA_add_custom_metadata.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/EMMA_add_custom_metadata.R +\name{EMMA_add_custom_metadata} +\alias{EMMA_add_custom_metadata} +\title{EMMA_add_custom_metadata} +\usage{ +EMMA_add_custom_metadata(res, extra = list()) +} +\arguments{ +\item{res}{A functional enrichment analysis results object as returned by +\code{EMMA_run()}} + +\item{extra}{A named list of user-defined metadata elements to store in +the \code{extra} field} +} +\value{ +The input result object with updated \code{EMMA_record} attribute +} +\description{ +Append or replace the \code{extra} field in the \code{EMMA_record} attribute +of a result object returned by \code{EMMA_run()}. This allows users to manually +provide additional annotation or contextual information that could not be +captured automatically +} +\examples{ +data("fea_res", package = "EMMA") +fea_res <- EMMA_add_custom_metadata(fea_res, extra = + list(note = "The background gene set list was all expressed genes in the assay")) +} diff --git a/man/EMMA_explain.Rd b/man/EMMA_explain.Rd index 03e177e..b3f5d4b 100644 --- a/man/EMMA_explain.Rd +++ b/man/EMMA_explain.Rd @@ -4,10 +4,28 @@ \alias{EMMA_explain} \title{EMMA_explain} \usage{ -EMMA_explain(res) +EMMA_explain(res, get_citation = TRUE) +} +\arguments{ +\item{res}{A functional enrichment analysis results object as returned by +\code{EMMA_run()}. Its attributes contain \code{EMMA_record}, which contains all +provenance information of the performed FEA} + +\item{get_citation}{Logical indicating whether to display the citations of +the packages used in the FEA. It only prints the citations in an interactive +session (e.g console). Defaults to \code{TRUE}} +} +\value{ +A character string describing how the FEA was performed using the +recorded metadata } \description{ This function 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. } +\examples{ +data("fea_res", package = "EMMA") +EMMA_explain(fea_res) + +} diff --git a/man/EMMA_freeze.Rd b/man/EMMA_freeze.Rd index 7f7c69f..496b7fd 100644 --- a/man/EMMA_freeze.Rd +++ b/man/EMMA_freeze.Rd @@ -4,9 +4,64 @@ \alias{EMMA_freeze} \title{EMMA_freeze} \usage{ -EMMA_freeze() +EMMA_freeze( + project = getwd(), + file = "renv.lock", + pkgs = loadedNamespaces(), + prompt = interactive(), + force = TRUE +) +} +\arguments{ +\item{project}{Character string corresponding to the path to the project +directory where the lockfile should be written. If the directory does not +exist, it will be created. It defaults to the current working directory} + +\item{file}{Character string referring to the name of the lockfile to generate. +It defaults to "renv.lock"} + +\item{pkgs}{Character vector of package names to snapshot. It defaults to all +currently loaded namespaces via \code{loadedNamespaces()}} + +\item{prompt}{Logical indicating whether to prompt before taking actions. +Defaults to \code{interactive()}} + +\item{force}{Logical indicating whether to force creation of the lockfile. +Defaults to \code{TRUE}} +} +\value{ +Invisibly returns the path to the generated lockfile. +The lockfile is written in JSON format and can be used with \code{renv::restore()} +to recreate the package environment } \description{ This function records the R environment during analysis runtime and generates -a lockfile that can be used with tools such as \code{renv} +a lockfile that can be used with tools such as \code{renv}. +By default, all currently loaded namespaces are recorded. +} +\details{ +This function calls \code{renv::snapshot()} with the specified packages. +The resulting lockfile can later be restored with \code{renv::restore()} to +recreate the same package environment. + +By default, the lockfile is created with \code{force = TRUE}, allowing snapshot +creation even if inconsistencies are detected in the environment. +} +\examples{ +# create a lockfile +if (requireNamespace("renv", quietly = TRUE)) { +tmp <- tempfile("emma_env") +dir.create(tmp) + +EMMA_freeze(project = tmp) + +# inspect generated files +list.files(tmp) + +# to restore the environment later +# renv::restore(project = tmp) +} +} +\seealso{ +\code{\link[renv]{snapshot}}, \code{\link[renv]{restore}} } diff --git a/man/EMMA_get_record.Rd b/man/EMMA_get_record.Rd new file mode 100644 index 0000000..7678801 --- /dev/null +++ b/man/EMMA_get_record.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/EMMA_get_record.R +\name{EMMA_get_record} +\alias{EMMA_get_record} +\title{EMMA_get_record} +\usage{ +EMMA_get_record(res) +} +\arguments{ +\item{res}{Functional Enrichment Analysis results (enrichResult, gseaResult ...) +generated by \code{EMMA_run()}} +} +\value{ +list of metadata recorded during FEA runtime +} +\description{ +EMMA_get_record +} +\examples{ +data("fea_res", package = "EMMA") +EMMA_get_record(fea_res) +} +\seealso{ +\code{\link[clusterProfiler:enrichGO]{enrichGO()}}, \code{\link[clusterProfiler:groupGO]{groupGO()}}, \code{\link[clusterProfiler:gseGO]{gseGO()}}, \code{\link[clusterProfiler:gseKEGG]{gseKEGG()}}, \code{\link[clusterProfiler:enrichKEGG]{enrichKEGG()}}, +\code{\link[gprofiler2:gost]{gost()}}, \code{\link[mosdef:run_cluPro]{run_cluPro()}}, \code{\link[mosdef:run_topGO]{run_topGO()}}, \code{\link[mosdef:run_goseq]{run_goseq()}} +} diff --git a/man/EMMA_run.Rd b/man/EMMA_run.Rd index bebb97d..c503f82 100644 --- a/man/EMMA_run.Rd +++ b/man/EMMA_run.Rd @@ -4,24 +4,52 @@ \alias{EMMA_run} \title{EMMA_run} \usage{ -EMMA_run(expr, envir = parent.frame()) +EMMA_run( + expr, + envir = parent.frame(), + store_session_info = TRUE, + args_form = c("evaluated", "unevaluated") +) } \arguments{ \item{expr}{A function call that performs functional enrichment analysis. The call is captured and executed by EMMA to record analysis parameters and -provenance information} +provenance information. Both bare calls (\code{enrichGO(...)}) and +namespace-qualified calls (\code{clusterProfiler::enrichGO(...)}) are supported.} \item{envir}{An environment in which to evaluate \code{expr}} + +\item{store_session_info}{Logical, indicating whether to store the output of +\code{sessionInfo()} or not. If \code{TRUE} (default), the session is stored in the +provenance record} + +\item{args_form}{A character string indicating whether to store the evaluated +or the unevaluated arguments in the provenance record. It default to +\code{"evaluated"}} } \value{ -Functional enrichment analysis results in the native format -returned by the original \code{expr} +The result object returned by the enrichment function in \code{expr}, +in standard format, with an additional \code{EMMA_record} attribute containing the +provenance information. Use \code{EMMA_get_record()} to retrieve this record } \description{ -This function executes functional enrichment analysis using existing tools -and captures the associated parameters and provenance information for the -analysis during runtime. +This function executes any functional enrichment analysis function and +attaches a provenance record (\code{EMMA_record}) describing the analysis. The +captured record includes the original call, metadata derived from the call +and its arguments, runtime information, and optionally the current session +information. +} +\details{ +\code{EMMA_run()} accepts both direct calls to known enrichment functions +(\code{enrichGO()}, \code{GSEA()}, \code{fgsea()} ...) and calls to wrapper functions that +internally invoke a know enrichment function. } \examples{ -EMMA_run(enrichGO(gene = rownames(de_res_IFNg_vs_naive), universe = universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP")) +data("de_res_IFNg_vs_naive", package = "EMMA") +data("gene_universe", package = "EMMA") +library(gprofiler2) + +EMMA_run(gost(query = de_res_IFNg_vs_naive$SYMBOL, organism = "hsapiens", + correction_method = "fdr", custom_bg = gene_universe, sources = "GO:BP"), + store_session_info = FALSE, args_form = "unevaluated") } diff --git a/man/EMMA_show.Rd b/man/EMMA_show.Rd index 14f6ee9..cd30ba3 100644 --- a/man/EMMA_show.Rd +++ b/man/EMMA_show.Rd @@ -11,14 +11,13 @@ EMMA_show(res) \code{EMMA_run()}} } \value{ -A list of the recorded information +\code{base::invisible()} } \description{ -This function prints the EMMA record associated with a functional enrichment -analysis object, including the executed call, parameters, and provenance -information. +This function displays a human-readable summary of the \code{EMMA_record} attribute +attached to a result object produced by \code{EMMA_run()} } \examples{ -res <- EMMA_run(enrichGO(gene = rownames(de_res_IFNg_vs_naive), universe = universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP")) -EMMA_show(res) +data("fea_res", package = "EMMA") +EMMA_show(fea_res) } diff --git a/man/de_res_IFNg_vs_naive.Rd b/man/de_res_IFNg_vs_naive.Rd new file mode 100644 index 0000000..dae7c8f --- /dev/null +++ b/man/de_res_IFNg_vs_naive.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/EMMA-data.R +\docType{data} +\name{de_res_IFNg_vs_naive} +\alias{de_res_IFNg_vs_naive} +\title{A sample \code{data.frame} containing Differential Expression Analysis, generated +with \code{DESeq2}} +\format{ +A \code{data.frame} object +} +\value{ +A sample \code{data.frame} object, extracted from \code{DESeq2} results +} +\description{ +A sample \code{data.frame} containing Differential Expression Analysis, generated +with \code{DESeq2} +} +\details{ +This \code{data.frame} object contains the results of a Differential +Expression Analysis performed on data from the \code{macrophage} package, more +precisely contrasting the counts from naive macrophage to those associated +with IFNg. + +The code to create said object can be found in the folder \verb{/inst/scripts} in +the EMMA package, the file is called \code{create_datasets_examples.R}. +} +\references{ +Alasoo, et al. "Shared genetic effects on chromatin and gene +expression indicate a role for enhancer priming in immune response", +Nature Genetics, January 2018 doi: 10.1038/s41588-018-0046-7. +} diff --git a/man/fea_res.Rd b/man/fea_res.Rd new file mode 100644 index 0000000..b08224c --- /dev/null +++ b/man/fea_res.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/EMMA-data.R +\docType{data} +\name{fea_res} +\alias{fea_res} +\title{A sample \code{list} containing Functional Enrichment Analysis results, +generated with \code{gprofiler2}} +\format{ +A \code{list} +} +\value{ +A sample \code{list} containing the FEA results \code{result} and \code{metadata}. +This results object has the \code{EMMA_record} attribute. +} +\description{ +A sample \code{list} containing Functional Enrichment Analysis results, +generated with \code{gprofiler2} +} +\details{ +This \code{list} object contains the result table and metadata of the +functional enrichment analysis (FEA) performed on the \code{macrophage} data, +specifically using the \code{gost()} function from the \code{gprofiler2} package, and +wrapped in \code{EMMA_run()} + +The code to create said object can be found in the folder \verb{/inst/scripts} in +the EMMA package, the file is called \code{create_datasets_examples.R}. +} +\references{ +Alasoo, et al. "Shared genetic effects on chromatin and gene +expression indicate a role for enhancer priming in immune response", +Nature Genetics, January 2018 doi: 10.1038/s41588-018-0046-7. +} diff --git a/man/gene_universe.Rd b/man/gene_universe.Rd new file mode 100644 index 0000000..aaa6359 --- /dev/null +++ b/man/gene_universe.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/EMMA-data.R +\docType{data} +\name{gene_universe} +\alias{gene_universe} +\title{A sample \verb{character vector} containing the background gene list used to +perform FEA on the \code{macrophage} dataset} +\format{ +A \verb{character vector} +} +\value{ +A sample \verb{character vector} containing the assay's \code{rownames} +of the \code{macrophage} data +} +\description{ +A sample \verb{character vector} containing the background gene list used to +perform FEA on the \code{macrophage} dataset +} +\details{ +This \verb{character vector} object contains the assay's \code{rownames} +of the \code{macrophage} data + +The code to create said object can be found in the folder \verb{/inst/scripts} in +the EMMA package, the file is called \code{create_datasets_examples.R}. +} +\references{ +Alasoo, et al. "Shared genetic effects on chromatin and gene +expression indicate a role for enhancer priming in immune response", +Nature Genetics, January 2018 doi: 10.1038/s41588-018-0046-7. +} diff --git a/pkgdown/extra.css b/pkgdown/extra.css new file mode 100644 index 0000000..c51e05b --- /dev/null +++ b/pkgdown/extra.css @@ -0,0 +1,97 @@ +/* +Developed and maintained by Kevin Rue-Albrecht (@kevinrue) +*/ + +/* +#0092ac blue +#00758a darker blue (active menu) +#c4d931 green (on blue) +#87b13f green (on white) +*/ + +.headroom { + background-color: #0092ac; +} + +.navbar-default .navbar-link { + color: #ffffff; +} + +.navbar-default .navbar-link:hover { + color: #c4d931; +} + +.navbar-default .navbar-nav>.active>a, +.navbar-default .navbar-nav>.active>a:hover, +.navbar-default .navbar-nav>.active>a:focus { + color: #c4d931; + background-color: #00758a; +} + +.navbar-default .navbar-nav>.open>a, +.navbar-default .navbar-nav>.open>a:hover, +.navbar-default .navbar-nav>.open>a:focus { + color: #c4d931; + background-color: #00758a; +} + +.dropdown-menu>.active>a, +.dropdown-menu>.active>a:hover, +.dropdown-menu>.active>a:focus { + color: #c4d931; + background-color: #00758a; +} + +.navbar-default .navbar-nav>li>a:hover, +.navbar-default .navbar-nav>li>a:focus { + color: #c4d931; +} + +.dropdown-menu>li>a:hover { + color: #87b13f; + background-color: #ffffff; +} + +.navbar-default .navbar-nav>li>a { + color: #ffffff; +} + +h1 { + color: #87b13f; +} + +h2 { + color: #1a81c2; +} + +h3 { + color: #1a81c2; + font-weight: bold; +} + +.btn-copy-ex { + color: #ffffff; + background-color: #0092ac; + border-color: #0092ac; +} + +.btn-copy-ex:hover { + color: #ffffff; + background-color: #00758a; + border-color: #00758a; +} + +.btn-copy-ex:active:focus { + color: #c4d931; + background-color: #00758a; + border-color: #0092ac; +} + +p>.fa, +p>.fas { + color: #0092ac; +} + +img { + width: auto; +} diff --git a/tests/testthat/setuptests_EMMA.R b/tests/testthat/setuptests_EMMA.R index 9d480ed..104f1ae 100644 --- a/tests/testthat/setuptests_EMMA.R +++ b/tests/testthat/setuptests_EMMA.R @@ -1,2 +1,12 @@ +suppressPackageStartupMessages({ + library("clusterProfiler") + library("org.Hs.eg.db") + library("mosdef") + library("topGO") + library("gprofiler2") +}) + data("de_res_IFNg_vs_naive", package = "EMMA") -data("universe", package = "EMMA") +data("gene_universe", package = "EMMA") +data("fea_res", package = "EMMA") + diff --git a/tests/testthat/test-EMMA.R b/tests/testthat/test-EMMA.R deleted file mode 100644 index 713438a..0000000 --- a/tests/testthat/test-EMMA.R +++ /dev/null @@ -1,26 +0,0 @@ -test_that("EMMA_run", { - - fea_res <- EMMA_run(enrichGO(gene = rownames(de_res_IFNg_vs_naive), - keyType = "ENSEMBL", - OrgDb = org.Hs.eg.db, - pAdjustMethod = "BH", - pvalueCutoff = 0.05, - qvalueCutoff = 0.1, - universe = universe, - readable = TRUE)) - expect_s4_class(fea_res, "enrichResult") - - expect_true("EMMA_record" %in% names(attributes(fea_res))) - - expect_type(attr(fea_res, "EMMA_record"), "list") - - expect_error(EMMA_run(enrichGO,gene = rownames(de_res_IFNg_vs_naive), - keyType = "ENSEMBL", - OrgDb = org.Hs.eg.db, - pAdjustMethod = "BH", - pvalueCutoff = 0.05, - qvalueCutoff = 0.1, - universe = universe, - readable = TRUE)) - -}) \ No newline at end of file diff --git a/tests/testthat/test-EMMA_explain.R b/tests/testthat/test-EMMA_explain.R new file mode 100644 index 0000000..f063cb6 --- /dev/null +++ b/tests/testthat/test-EMMA_explain.R @@ -0,0 +1,57 @@ +test_that("EMMA_explain", { + res <- data.frame(ID = "GO:0000001") + + attr(res, "EMMA_record") <- list( + method = list( + function_name = "enrichGO", + package_name = "clusterProfiler", + package_version = "4.10.0", + wrapper = FALSE + ), + input = list( + arguments = list( + gene = c("g1", "g2"), + universe = c("g1", "g2", "g3", "g4"), + pAdjustMethod = "BH" + ) + ), + annotation = list( + gene_set_db = "GO", + gene_set_db_version = "3.18.0" + ) + ) + + expect_message( + txt <- EMMA_explain(res), + "You can always complete your text" + ) + + expect_match(txt, "Functional Enrichment Analysis was performed using the enrichGO\\(\\) function") + expect_match(txt, "from the clusterProfiler package") + expect_match(txt, "with the GO database") + expect_match(txt, "A custom background gene set was provided \\(n = 4\\)") + expect_match(txt, "Multiple testing correction was performed using the BH method") + expect_false(grepl("NA", txt)) + + de_res <- de_res_IFNg_vs_naive[!(is.na(de_res_IFNg_vs_naive$padj)) & de_res_IFNg_vs_naive$padj <= 0.05, ] + + + res_topGO <- run_topGO(de_genes = rownames(de_res), + bg_genes = gene_universe, + ontology = "BP", + gene_id = "ENSEMBL", + mapping = "org.Hs.eg.db", + add_gene_to_terms = TRUE, + do_padj = TRUE) |> EMMA_run() + + expect_message( + txt <- EMMA_explain(res_topGO), + "You can always complete your text" + ) + + expect_match(txt, + "Functional Enrichment Analysis was performed using a wrapper function run_topGO\\(\\)") + + expect_match(txt,"Multiple testing correction was applied") + +}) diff --git a/tests/testthat/test-EMMA_freeze.R b/tests/testthat/test-EMMA_freeze.R new file mode 100644 index 0000000..2e54fb0 --- /dev/null +++ b/tests/testthat/test-EMMA_freeze.R @@ -0,0 +1,34 @@ +test_that("EMMA_freeze", { + + tmp <- tempfile("emma_env") + dir.create(tmp) + + lockfile <- file.path(tmp, "renv.lock") + + res_path <- EMMA_freeze( + project = tmp, + file = "renv.lock", + pkgs = "stats", + prompt = TRUE, + force = TRUE + ) + + expect_true(dir.exists(tmp)) + + expect_true(file.exists(lockfile)) + + expect_identical(res_path, lockfile) + + expect_error(EMMA_freeze(project = 2)) + + expect_error(EMMA_freeze(file = NULL)) + + expect_error(EMMA_freeze(file = NA_character_)) + + expect_error(EMMA_freeze(pkgs = 1:3)) + + expect_error(EMMA_freeze(promt = "no")) + + expect_error(EMMA_freeze(force = "yes")) + +}) diff --git a/tests/testthat/test-EMMA_get_record.R b/tests/testthat/test-EMMA_get_record.R new file mode 100644 index 0000000..358897f --- /dev/null +++ b/tests/testthat/test-EMMA_get_record.R @@ -0,0 +1,27 @@ +test_that("EMMA_get_record", { + expect_true(is.list(EMMA_get_record(fea_res))) + expect_length(EMMA_get_record(fea_res), 7) + + fea_no_emma <- mosdef::run_cluPro(de_genes = + rownames(de_res_IFNg_vs_naive), + bg_genes = gene_universe, + mapping = "org.Hs.eg.db", + keyType = "ENSEMBL", + ont = "BP", + pAdjustMethod = "BH") + + expect_error(EMMA_get_record(fea_no_emma)) + + attr(fea_no_emma, "EMMA_record") <- c() + + expect_error(EMMA_get_record(fea_no_emma)) + + expect_error(EMMA_get_record("guiga")) + + attr(fea_res, "EMMA_record") <- NULL + expect_error(EMMA_get_record(fea_res)) + + attr(fea_res, "EMMA_record") <- "i am a corrupted rec" + expect_error(EMMA_get_record(fea_res)) + +}) \ No newline at end of file diff --git a/tests/testthat/test-EMMA_metadata.R b/tests/testthat/test-EMMA_metadata.R new file mode 100644 index 0000000..023a21f --- /dev/null +++ b/tests/testthat/test-EMMA_metadata.R @@ -0,0 +1,129 @@ +test_that("test metadata content & structure", { + + fea_res <- EMMA_run(mosdef::run_cluPro(de_genes = + rownames(de_res_IFNg_vs_naive), + bg_genes = gene_universe, + mapping = "org.Hs.eg.db", + keyType = "ENSEMBL", + ont = "BP", + pAdjustMethod = "BH")) + + emma_rec <- attr(fea_res, "EMMA_record") + + org <- emma_rec$annotation$organism + + db <- emma_rec$annotation$gene_set_db + + expect_true(is.list(emma_rec$annotation)) + expect_true(is.list(emma_rec$method)) + expect_true(is.list(emma_rec$input)) + + expect_identical(org, "Homo sapiens") + expect_identical(db, "GO") + + expect_true(emma_rec$method$wrapper) + + expect_null(.EMMA_find_original_wrapped_fun("mosdef::run_cluPro(de_genes = rownames(de_res_IFNg_vs_naive), + bg_genes = gene_universe, + mapping = 'org.Hs.eg.db', + keyType = 'ENSEMBL', + ont = 'BP')")) + + + custom <- .EMMA_classify_call(substitute(summary(EMMA_get_record(fea_res)))) + + expect_equal(custom$type, "custom") + + + expect_warning(res <- EMMA_run(mosdef::run_goseq(de_genes = rownames(de_res_IFNg_vs_naive), + bg_genes = gene_universe, + mapping = "org.Hs.eg.db", + id = "ensGene", + genome = "hg19"))) + + + res <- EMMA_add_custom_metadata(res, + extra = list( + note = "The background gene set list was all expressed genes in the assay")) + + rec <- EMMA_get_record(res) + + expect_equal( + rec$extra$note, + "The background gene set list was all expressed genes in the assay" + ) + + expect_null(rec$extra$notfoud) + + expect_error(res <- EMMA_add_custom_metadata(res, + extra = list("The background gene set list was all expressed genes in the assay")) + ) + + expect_error(res <- EMMA_add_custom_metadata(res, + extra = "The background gene set list was all expressed genes in the assay")) + + + fea_res <- gprofiler2::gost(query = de_res_IFNg_vs_naive$SYMBOL, + organism = "hsapiens", + correction_method = "fdr", + custom_bg = gene_universe) |> EMMA_run() + + rec <- EMMA_get_record(fea_res) + + expect_true(length(rec$annotation$gene_set_db_version) != 1) + + + custom_fun <- function(gene, OrgDb) { + clusterProfiler::groupGO( + gene = gene, + OrgDb = OrgDb, + keyType = "ENSEMBL", + ont = "BP", + level = 2, + readable = FALSE + ) + } + + expect_warning(wrapper <- EMMA_run(custom_fun(rownames(de_res_IFNg_vs_naive), + org.Hs.eg.db))) + + expect_equal(EMMA_get_record(wrapper)$method$function_name, "custom_fun") + expect_equal(EMMA_get_record(wrapper)$method$wrapped_package, "clusterProfiler") + expect_equal(EMMA_get_record(wrapper)$method$wrapped_function, "groupGO") + expect_true(EMMA_get_record(wrapper)$method$wrapper) + + + empty <- EMMA_run(summary(rec$method)) + + expect_null(EMMA_get_record(empty)$annotation$organism) + expect_null(EMMA_get_record(empty)$annotation$gene_set_db) + expect_null(EMMA_get_record(empty)$annotation$gene_set_db_version) + + + geneList <- de_res_IFNg_vs_naive$log2FoldChange + names(geneList) <- rownames(de_res_IFNg_vs_naive) + geneList <- sort(geneList, decreasing = TRUE) + + expect_warning(gse_res <- gseGO( geneList = geneList, keyType = "ENSEMBL", + OrgDb = org.Hs.eg.db, ont = "BP", + minGSSize = 100, maxGSSize = 500, + pvalueCutoff = 0.05, verbose = FALSE, + pAdjustMethod = "BH") |> EMMA_run()) + + expect_equal(EMMA_get_record(gse_res)$method$function_name, "gseGO") + expect_null(EMMA_get_record(gse_res)$method$wrapped_function) + expect_null(EMMA_get_record(gse_res)$method$wrapped_package) + expect_equal(EMMA_get_record(gse_res)$annotation$gene_set_db, "GO") + + # pretending a function that doesn't exist in any package, so fallback + # to empty metadata + expect_warning(user_only_fea <- summary(as.data.frame(gse_res)) |> + EMMA_run(store_session_info = FALSE, args_form = "unevaluated")) + + expect_null(EMMA_get_record(user_only_fea)$annotation$organism) + expect_null(EMMA_get_record(user_only_fea)$annotation$gene_set_db) + expect_null(EMMA_get_record(user_only_fea)$annotation$gene_set_db_version) + +}) + + diff --git a/tests/testthat/test-EMMA_run.R b/tests/testthat/test-EMMA_run.R new file mode 100644 index 0000000..ef70d52 --- /dev/null +++ b/tests/testthat/test-EMMA_run.R @@ -0,0 +1,84 @@ +test_that("EMMA_run", { + + fea_res <- EMMA_run(enrichGO(gene = rownames(de_res_IFNg_vs_naive), + keyType = "ENSEMBL", + OrgDb = org.Hs.eg.db, + pAdjustMethod = "BH", + pvalueCutoff = 0.05, + qvalueCutoff = 0.1, + universe = gene_universe, + readable = TRUE)) + expect_s4_class(fea_res, "enrichResult") + + expect_true("EMMA_record" %in% names(attributes(fea_res))) + + expect_true(is.list(attr(fea_res, "EMMA_record"))) + + expect_true(length(attr(fea_res, "EMMA_record")) == 7) + + info <- .EMMA_capture_call_info(substitute(enrichGO(gene = rownames(de_res_IFNg_vs_naive), + keyType = "ENSEMBL", + OrgDb = org.Hs.eg.db, + pAdjustMethod = "BH", + pvalueCutoff = 0.05, + qvalueCutoff = 0.1, + universe = gene_universe, + readable = TRUE))) + + expect_type(info, "list") + expect_length(info, 5) + expect_length(info$arg_list, 8) + + expect_error(EMMA_run(enrichGO,gene = rownames(de_res_IFNg_vs_naive), + keyType = "ENSEMBL", + OrgDb = org.Hs.eg.db, + pAdjustMethod = "BH", + pvalueCutoff = 0.05, + qvalueCutoff = 0.1, + universe = gene_universe, + readable = TRUE)) + + expect_error(EMMA_run("enrichGO(gene = rownames(de_res_IFNg_vs_naive), + keyType = 'ENSEMBL', + OrgDb = org.Hs.eg.db)")) + + expect_error(EMMA_run(gene = rownames(de_res_IFNg_vs_naive), + keyType = "ENSEMBL", + OrgDb = org.Hs.eg.db, + pAdjustMethod = "BH", + pvalueCutoff = 0.05, + qvalueCutoff = 0.1, + universe = gene_universe, + readable = TRUE)) + + expect_warning(EMMA_run(enrichGO(gene = rownames(de_res_IFNg_vs_naive), + keyType = "ENSEMBL", + OrgDb = org.Hs.eg.db, + pvalueCutoff = 0.05, + qvalueCutoff = 0.1, + universe = gene_universe, + readable = TRUE))) + + expect_warning(EMMA_run(enrichGO(gene = rownames(de_res_IFNg_vs_naive), + keyType = "ENSEMBL", + OrgDb = org.Hs.eg.db, + pvalueCutoff = 0.05, + qvalueCutoff = 0.1, + pAdjustMethod = "BH", + readable = TRUE))) + + expect_error(EMMA_run(groupGO(gene = rownames(de_res_IFNg_vs_naive), + OrgDb = org.Hs.eg.db, + keyType = "ENSEMBL", + ont = "CC", + level = 2), envir = getwd())) + + expect_error(EMMA_run(groupGO(gene = rownames(de_res_IFNg_vs_naive), + OrgDb = org.Hs.eg.db, + keyType = "ENSEMBL", + ont = "CC", + level = 2), store_session_info = "yes")) + +}) + + diff --git a/tests/testthat/test-EMMA_show.R b/tests/testthat/test-EMMA_show.R new file mode 100644 index 0000000..b5fb514 --- /dev/null +++ b/tests/testthat/test-EMMA_show.R @@ -0,0 +1,82 @@ +test_that("EMMA_show", { + + fea_res <- data.frame(ID = "term1") + + expect_warning(EMMA_show(fea_res)) + + + res <- data.frame(ID = c("term1", "term2")) + + attr(res, "EMMA_record") <- list( + method = list( + call = substitute(fake_enrich(gene = genes)), + wrapper = FALSE, + package_name = "fakepkg", + package_version = "1.0.0" + ), + annotation = list( + organism = "Homo sapiens", + gene_set_db = "GO", + gene_set_db_version = "3.22.0" + ) + ) + + expect_message( + expect_output( + EMMA_show(res), + "Number of Pathways" + ), + "Found EMMA record" + ) + + expect_output(EMMA_show(res), "Package") + expect_output(EMMA_show(res), "Organism") + expect_output(EMMA_show(res), "Gene set library") + + + fea2 <- list(fea_res, + res) + + attr(fea2, "EMMA_record") <- list( + method = list( + call = substitute(fake_enrich(gene = genes)), + wrapper = TRUE, + package_name = NA, + package_version = NA + ), + annotation = list( + organism = "Homo sapiens", + gene_set_db = "KEGG", + gene_set_db_version = NA + ) + ) + expect_output(EMMA_show(fea2), "Number of FEAs: 2") + expect_output(EMMA_show(fea2), "- FEA_1") + expect_output(EMMA_show(fea2),"FEA_2") + + + fea_gost <- fea2 + names(fea_gost)[2] <- "result" + + attr(fea_gost, "EMMA_record") <- list( + method = list( + call = substitute(gost(query = de_res_IFNg_vs_naive$SYMBOL, + organism = "hsapiens", + correction_method = "fdr", + custom_bg = gene_universe, + sources = "GO:BP")), + wrapper = FALSE, + package_name = "gprofiler2", + package_version = NA + ), + annotation = list( + organism = "Homo sapiens", + gene_set_db = "GO", + gene_set_db_version = NA + ) + ) + + expect_output(EMMA_show(fea_gost), "Number of Pathways: 2") + + +}) diff --git a/vignettes/EMMA_bibliography.bib b/vignettes/EMMA_bibliography.bib index 1c783b0..6238d17 100644 --- a/vignettes/EMMA_bibliography.bib +++ b/vignettes/EMMA_bibliography.bib @@ -13,4 +13,95 @@ @article{Alasoo2018 url = {http://dx.doi.org/10.1038/s41588-018-0046-7}, volume = {50}, year = {2018} -} \ No newline at end of file +} + +@article{Khatri2012, + title = {Ten Years of Pathway Analysis: Current Approaches and Outstanding Challenges}, + volume = {8}, + ISSN = {1553-7358}, + url = {http://dx.doi.org/10.1371/journal.pcbi.1002375}, + DOI = {10.1371/journal.pcbi.1002375}, + number = {2}, + journal = {PLoS Computational Biology}, + publisher = {Public Library of Science (PLoS)}, + author = {Khatri, Purvesh and Sirota, Marina and Butte, Atul J.}, + editor = {Ouzounis, Christos A.}, + year = {2012}, + month = Feb, + pages = {e1002375} +} + +@article{Subramanian2005, + title = {Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles}, + volume = {102}, + ISSN = {1091-6490}, + url = {http://dx.doi.org/10.1073/pnas.0506580102}, + DOI = {10.1073/pnas.0506580102}, + number = {43}, + journal = {Proceedings of the National Academy of Sciences}, + publisher = {Proceedings of the National Academy of Sciences}, + author = {Subramanian, Aravind and Tamayo, Pablo and Mootha, Vamsi K. and Mukherjee, Sayan and Ebert, Benjamin L. and Gillette, Michael A. and Paulovich, Amanda and Pomeroy, Scott L. and Golub, Todd R. and Lander, Eric S. and Mesirov, Jill P.}, + year = {2005}, + month = Sept, + pages = {15545–15550} +} + +@article{Wijesooriya2022, + title = {Urgent need for consistent standards in functional enrichment analysis}, + volume = {18}, + ISSN = {1553-7358}, + url = {http://dx.doi.org/10.1371/journal.pcbi.1009935}, + DOI = {10.1371/journal.pcbi.1009935}, + number = {3}, + journal = {PLOS Computational Biology}, + publisher = {Public Library of Science (PLoS)}, + author = {Wijesooriya, Kaumadi and Jadaan, Sameer A. and Perera, Kaushalya L. and Kaur, Tanuveer and Ziemann, Mark}, + editor = {Kemp, Melissa L.}, + year = {2022}, + month = Mar, + pages = {e1009935} +} + +@article{Brazma2001, + title = {Minimum information about a microarray experiment (MIAME)—toward standards for microarray data}, + volume = {29}, + ISSN = {1546-1718}, + url = {http://dx.doi.org/10.1038/ng1201-365}, + DOI = {10.1038/ng1201-365}, + number = {4}, + journal = {Nature Genetics}, + publisher = {Springer Science and Business Media LLC}, + author = {Brazma, Alvis and Hingamp, Pascal and Quackenbush, John and Sherlock, Gavin and Spellman, Paul and Stoeckert, Chris and Aach, John and Ansorge, Wilhelm and Ball, Catherine A. and Causton, Helen C. and Gaasterland, Terry and Glenisson, Patrick and Holstege, Frank C.P. and Kim, Irene F. and Markowitz, Victor and Matese, John C. and Parkinson, Helen and Robinson, Alan and Sarkans, Ugis and Schulze-Kremer, Steffen and Stewart, Jason and Taylor, Ronald and Vilo, Jaak and Vingron, Martin}, + year = {2001}, + month = Dec, + pages = {365–371} +} + +@article{Abassi2026, + title = {DeeDeeExperiment: building an infrastructure for integrating and managing omics data analysis results in R/Bioconductor}, + volume = {42}, + ISSN = {1367-4811}, + url = {http://dx.doi.org/10.1093/bioinformatics/btag157}, + DOI = {10.1093/bioinformatics/btag157}, + number = {4}, + journal = {Bioinformatics}, + publisher = {Oxford University Press (OUP)}, + author = {Abassi, Najla and Schwarz, Lea and Filippi, Edoardo and Marini, Federico}, + editor = {Schwartz, Russell}, + year = {2026}, + month = Mar +} + + + + + + + + + + + + + + diff --git a/vignettes/EMMA_workflow.png b/vignettes/EMMA_workflow.png new file mode 100644 index 0000000..dc8ba88 Binary files /dev/null and b/vignettes/EMMA_workflow.png differ diff --git a/vignettes/Introduction_to_EMMA.Rmd b/vignettes/Introduction_to_EMMA.Rmd index 3e45a11..b9ac952 100644 --- a/vignettes/Introduction_to_EMMA.Rmd +++ b/vignettes/Introduction_to_EMMA.Rmd @@ -1,6 +1,6 @@ --- title: > - The `EMMA` User's Guide + The `EMMA` User's Guide - Enrichment Method MAtters author: - name: Najla Abassi affiliation: @@ -34,32 +34,90 @@ editor_options: bibliography: EMMA_bibliography.bib --- -```{r, include = FALSE} + +```{r knitr, include = FALSE} knitr::opts_chunk$set( - comment = "#>", - error = FALSE, - warning = FALSE, - eval = FALSE, - message = FALSE + comment = "#>", + error = FALSE, + warning = FALSE, + eval = TRUE, + message = FALSE ) ``` # Introduction {#introduction} -This vignette describes how to use the `r BiocStyle::Biocpkg("EMMA")` -package to perform Functional Enrichment Analysis using existing tools (e.g. -`clusterProfiler`, `topGO`, `Enrichr`, `gprofiler2`), while systematically -capturing analysis parameters and provenance information during runtime. +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) # Getting started {#gettingstarted} -To install the development version of this package, start R and enter: +To install this package, start R and enter: ```{r install, eval = FALSE} -library("remotes") -remotes::install_github("imbeimainz/EMMA", - dependencies = TRUE, - build_vignettes = TRUE) +if (!requireNamespace("BiocManager", quietly = TRUE)) { + install.packages("BiocManager")} + +BiocManager::install("EMMA") ``` Once installed, the package can be loaded and attached to the current workspace @@ -69,13 +127,15 @@ as follows: library("EMMA") ``` -# `EMMA` on the `macrophage` dataset +# 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] [doi:10.1038/s41588-018-0046-7](https://doi.org/10.1038/s41588-018-0046-7). +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 @@ -95,17 +155,20 @@ 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. -# Get a list of Differentially Expressed Genes +## 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 dds} +```{r differential_expression} # load data data(gse, "macrophage") # set up design @@ -113,24 +176,18 @@ 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, ] -dds_macrophage -``` -We then run the main `DESeq()` function and check the `resultsNames` being -generated: +# 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, ] -```{r run_deseq} # run DESeq dds_macrophage <- DESeq(dds_macrophage) -# contrasts -resultsNames(dds_macrophage) -``` - -Let's extract the DE results for each contrast. We will do three contrasts: -```{r get_de_res} # get de res for 1st contrast IFNg_vs_naive <- results(dds_macrophage, contrast = c("condition", "IFNg", "naive"), @@ -139,131 +196,394 @@ 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 -IFNg_vs_naive -# get de res for 2nd contrast -Salm_vs_naive <- results(dds_macrophage, +# get de res for 2st contrast +SL1344_vs_naive <- results(dds_macrophage, contrast = c("condition", "SL1344", "naive"), lfcThreshold = 1, alpha = 0.05) -Salm_vs_naive <- lfcShrink(dds_macrophage, coef = "condition_SL1344_vs_naive", - res = Salm_vs_naive, +SL1344_vs_naive <- lfcShrink(dds_macrophage, coef = "condition_SL1344_vs_naive", + res = SL1344_vs_naive, type = "apeglm") -Salm_vs_naive$SYMBOL <- rowData(dds_macrophage)$SYMBOL -Salm_vs_naive +SL1344_vs_naive$SYMBOL <- rowData(dds_macrophage)$SYMBOL -# get de res for 3rd contrast -Both_vs_naive <- results(dds_macrophage, - contrast = c("condition", "IFNg_SL1344", "naive"), - lfcThreshold = 1, alpha = 0.05) -Both_vs_naive <- lfcShrink(dds_macrophage, - coef = "condition_IFNg_SL1344_vs_naive", - res = Both_vs_naive, - type = "apeglm") -Both_vs_naive$SYMBOL <- rowData(dds_macrophage)$SYMBOL -Both_vs_naive +# 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) ``` -Now let's sort the results based on the adjusted p-value: - -```{r} -# have a list of de results per contrast -contrast_list <- list(IFNg_vs_naive = IFNg_vs_naive, - Salm_vs_naive = Salm_vs_naive, - Both_vs_naive = Both_vs_naive) - -de_res_list <- list() -# sort -for (i in names(contrast_list)) { - de_res <- as.data.frame(contrast_list[[i]]) - de_res <- de_res[order(de_res$padj), ] - de_res <- de_res[!(is.na(de_res$padj)) & de_res$padj <= 0.05, ] - - de_res_list[[i]] <- de_res -} +## 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 -# set universe -universe <- rownames(dds_macrophage) +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) ``` -# Perform Functional Enrichment Analysis (FEA) +# Entering the `EMMA` framework + +From now on, we will complement the existing workflow with the functionality +provided by `EMMA`. -Now that we have a list of DE genes for each contrast, we can perform Functional -Enrichment Analysis. For the first contrast `IFNg_vs_naive`, we will use the -function `enrichGO()` from `r BiocStyle::Biocpkg("clusterProfiler")` + +## `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} -# get de res from the list -de_res_IFNg_vs_naive <- de_res_list$IFNg_vs_naive - # perform FEA, but with EMMA! -fea_res <- enrichGO(gene = de_res_IFNg_vs_naive$SYMBOL, - keyType = "SYMBOL", +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} -fea_res <- EMMA_run(enrichGO(gene = de_res_IFNg_vs_naive$SYMBOL, - keyType = "ENSEMBL", - OrgDb = org.Hs.eg.db, - ont = "BP", - pAdjustMethod = "BH", - pvalueCutoff = 0.05, - qvalueCutoff = 0.1, - readable = TRUE)) - -# check res -fea_res +# 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, so we get warned about that. +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 check what EMMA captured while we ran the analysis, we use `EMMA_show()`: +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 aims to capture the following elements: - - * the full call used to perform the FEA - * the package used and its version - * the evaluated arguments used - * the gene set library (GO, KEGG, Reactome ...) and its version - * the organism type - * the runtime, i.e. when the analysis started - * Information about the current R session (`sessionInfo`) +`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_freeze()` records the R environment during analysis runtime and generates -a lockfile that can be used with tools such as `renv`, to facilitate -reproducible reconstruction of the analysis environment. +## `EMMA_get_record()`: Retrieving the recorded information -```{r EMMA_freeze} -EMMA_freeze() +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} -EMMA_explain(fea_res) +```{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 ``` -# Session info {.unnumbered .smaller} +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() ```