Function suba() for enrichment analyses of metabolic subsystems

Author: PiotrTymoszuk

DESCRIPTION

@@ -14,10 +14,10 @@ Authors@R: c(
            comment = c(ORCID = "0000-0002-0398-6034"))
            )
 Imports: 
+    fastTest,
     furrr,
     future,
     ggplot2,
-    microViz,
     purrr,
     Rcpp,
     rlang,

NAMESPACE

@@ -27,6 +27,7 @@ export(plot_mc)
 export(plot_numbers)
 export(reactDB)
 export(select.actiData)
+export(suba)
 export(transform_estimates)
 importFrom(Rcpp,sourceCpp)
 importFrom(dplyr,all_of)
@@ -41,7 +42,9 @@ importFrom(dplyr,relocate)
 importFrom(dplyr,select)
 importFrom(dplyr,summarise)
 importFrom(dplyr,tibble)
+importFrom(dplyr,transmute)
 importFrom(dplyr,ungroup)
+importFrom(fastTest,f_enrichment)
 importFrom(furrr,furrr_options)
 importFrom(furrr,future_map)
 importFrom(future,plan)
@@ -60,7 +63,7 @@ importFrom(ggplot2,labs)
 importFrom(ggplot2,position_jitter)
 importFrom(ggplot2,scale_fill_manual)
 importFrom(ggplot2,scale_x_continuous)
-importFrom(microViz,theme_micro)
+importFrom(ggplot2,theme_classic)
 importFrom(purrr,compact)
 importFrom(purrr,map)
 importFrom(purrr,map2_chr)

R/enrichment_funs.R

@@ -0,0 +1,134 @@
+# functions for subsystem enrichment analyses
+
+#' Subsystem enrichment analysis for differentially regulated reactions.
+#'
+#' @description
+#' Function `suba()` performs enrichment analyses for member reactions of
+#' metabolic subsystems among significantly activated and inhibited reactions in
+#' an \code{\link{actiData}} object.
+#' Please note, that prior to the analyses, significantly activated and
+#' inhibited reactions need to by determined by calling
+#' \code{\link{identify_regulated}}.
+#'
+#' @details
+#' The analyses are done with \code{\link[fastTest]{f_enrichment}} with two methods:
+#' Fisher's exact test (`type = "fisher`) or by random drawing from the entire
+#' reaction pool (`type = "random"`).
+#'
+#' @return a data frame with enrichment analysis results:
+#' numbers of reactions in the activated and inhibited reaction sets and in
+#' the subsystems, enrichment magnitude measured by odds ratio (OR) statistic,
+#' raw p values, and p values adjusted for multiple testing.
+#'
+#' @param x a \code{\link{actiData}} object.
+#' @param type type of the hypothesis test: Fisher's exact test or random drawing
+#' from the whole reaction pool.
+#' @param ... additional arguments passed to \code{\link[fastTest]{f_enrichment}}
+#' which specify e.g. type of confidence intervals for the OR statistic and
+#' number of random draws, and serial/parallel execution of the testing algorithm.
+#'
+#' @seealso [get_regulation()]
+#'
+#' @export
+
+  suba <- function(x,
+                   type = c("fisher", "random"), ...) {
+
+    ## the input control -------
+
+    if(!is_actiData(x)) {
+
+      stop("`x` has to be an `actiData` object created e.g. with `get_regulation()`",
+           call. = FALSE)
+
+    }
+
+    type <- match.arg(type[[1]], c("fisher", "random"))
+
+    if(!"regulation" %in% names(x[["reg"]])) {
+
+      stop(paste("No regulation status in the `x` object.",
+                 "Please call `identify_regulated()`",
+                 "prior to launching `suba()`."),
+           call. = FALSE)
+
+    }
+
+    ## the testing data ---------
+
+    ### vectors of significantly regulated reactions
+
+    signif_reactions <-
+      filter(x$reg,
+             .data[["regulation"]] %in% c("activated", "inhibited"))
+
+    if(nrow(signif_reactions) == 0) {
+
+      stop("No significantly regulated reactions.", call. = FALSE)
+
+    }
+
+    signif_reactions <-
+      split(signif_reactions[["id"]],
+            f = signif_reactions[["regulation"]])
+
+    signif_reactions <-
+      compact(map(signif_reactions,
+                  function(x) if(length(x) == 0) NULL else x))
+
+    ### dictionary of reactions in the subsystems
+    ### vector with identifiers of all reactions
+
+    reaction_dict <- split(x[["reg"]][["id"]],
+                           x[["reg"]][["subsystem"]])
+
+    all_reactions <- unique(x[["reg"]][["id"]])
+
+    ## testing and customizing the results ---------
+
+    result <- list()
+
+    for(i in names(signif_reactions)) {
+
+      result[[i]] <-  f_enrichment(x = signif_reactions[[i]],
+                                   type = type,
+                                   dict = reaction_dict,
+                                   all = all_reactions,
+                                   as_data_frame = TRUE,
+                                   adj_method = "BH", ...)
+
+    }
+
+    reaction_set <- NULL
+
+    result <-
+      map2_dfr(result, names(result),
+               ~mutate(.x,
+                       reaction_set = factor(.y,
+                                             c("activated", "inhibited"))))
+
+    n_total_reaction_set <- NULL
+    subsystem <- NULL
+    n_total_subsystem <- NULL
+    n_intersect <- NULL
+    or <- NULL
+    lower_ci <- NULL
+    upper_ci <- NULL
+    p_value <- NULL
+    p_adjusted <- NULL
+
+    transmute(as_tibble(result),
+              reaction_set = .data[["reaction_set"]],
+              n_total_reaction_set = .data[["n_x_total"]],
+              subsytem = .data[["entry_name"]],
+              n_total_subsystem = .data[["n_entry"]],
+              n_intersect = .data[["n_intersect"]],
+              or = .data[["or"]],
+              lower_ci = if(type == "fisher") NULL else .data[["lower_ci"]],
+              upper_ci = if(type == "fisher") NULL else .data[["upper_ci"]],
+              p_value = .data[["p_value"]],
+              p_adjusted = .data[["p_adjusted"]])
+
+  }
+
+# END ----------

R/imports.R

@@ -2,6 +2,7 @@
 
 #' @importFrom dplyr filter
 #' @importFrom dplyr mutate
+#' @importFrom dplyr transmute
 #' @importFrom dplyr count
 #' @importFrom dplyr tibble
 #' @importFrom dplyr select
@@ -45,6 +46,7 @@
 #' @importFrom ggplot2 position_jitter
 #' @importFrom ggplot2 is.theme
 #' @importFrom ggplot2 labs
+#' @importFrom ggplot2 theme_classic
 #'
 #' @importFrom rlang set_names
 #' @importFrom rlang `.data`
@@ -77,7 +79,7 @@
 #' @importFrom stats pnorm
 #' @importFrom stats na.omit
 #'
-#' @importFrom microViz theme_micro
+#' @importFrom fastTest f_enrichment
 #'
 #' @useDynLib biggrExtra
 

R/plotting_funs.R

@@ -70,7 +70,7 @@
                           hline_color = "black",
                           hline_type = "dashed",
                           hline_alpha = 1,
-                          cust_theme = microViz::theme_micro(),
+                          cust_theme = theme_classic(),
                           plot_title = NULL,
                           plot_subtitle = NULL,
                           x_lab = "regulation estimate",
@@ -241,7 +241,7 @@
                       zero_line_color = "black",
                       zero_line_type = "dashed",
                       zero_line_alpha = 1,
-                      cust_theme = microViz::theme_micro(),
+                      cust_theme = theme_classic(),
                       plot_title = NULL,
                       plot_subtitle = NULL,
                       x_lab = NULL,
@@ -435,7 +435,7 @@
                            show_all = FALSE,
                            labeller_fun = identity,
                            show_n_total = TRUE,
-                           cust_theme = microViz::theme_micro(),
+                           cust_theme = theme_classic(),
                            plot_title = NULL,
                            plot_subtitle = NULL,
                            x_lab = NULL,

inst/examples/development.R

@@ -183,8 +183,10 @@
          point_alpha = 0.15,
          plot_title = "Oxidative phosphorylation")
 
+  tst_mc_estimates <- tst_mc_estimates %>%
+    identify_regulated
+
   tst_mc_estimates %>%
-    identify_regulated %>%
     plot(plot_type = "numbers",
          scale = "percent",
          type = "plus_minus",
@@ -197,5 +199,14 @@
 
 # Subsystem enrichment ---------
 
+  tst_enrichment <- tst_mc_estimates %>%
+    suba(type = "random",
+         n_iter = 1000,
+         .parallel = TRUE)
+
+  tst_enrichment %>%
+    filter(p_adjusted < 0.05,
+           or >= 1.44) %>%
+    blast(reaction_set)
 
 # END --------