Update GO analysis

R/GSEA.R

@@ -108,11 +108,12 @@ runGSEA <- function(
     return(gsea)
 }
 
-#' Run Gene Ontology enrichment analysis on metagenes
+#' Run Gene Ontology enrichment analysis on differentially expressed genes.
 #' @description
 #' This function forms genesets basing on the differential expression result,
 #' and calls gene ontology (GO) analysis method provided by gprofiler2.
-#' @param result Data frame of unfiltered output from \code{\link{runWilcoxon}}.
+#' @param result Data frame of unfiltered output from \code{\link{runMarkerDEG}}
+#' or \code{\link{runPairwiseDEG}}.
 #' @param group Selection of one group available from \code{result$group}.
 #' Default \code{NULL} uses all groups involved in DE \code{result} table.
 #' @param useBg Logical, whether to set all genes involved in DE analysis
@@ -124,27 +125,38 @@ runGSEA <- function(
 #' @param logFCThresh The log2FC threshold above which the genes will be used.
 #' Default \code{1}.
 #' @param padjThresh The adjusted p-value threshold less than which the genes
-#' will be used. Default \code{0.01}.
-#' @param ... Additional arguments passed to \code{\link[gprofiler2]{gost}}.
+#' will be used. Default \code{0.05}.
+#' @param splitReg Whether to have queries of both up-regulated and
+#' down-regulated genes for each group. Default \code{FALSE} only queries
+#' up-regulated genes and should be preferred when \code{result} comes from
+#' marker detection test. When \code{result} comes from group-to-group DE test,
+#' it is recommended to set \code{splitReg = TRUE}.
+#' @param ... Additional arguments passed to
+#' \code{gprofiler2::\link[gprofiler2]{gost}}. Arguments \code{query},
+#' \code{custom_bg}, \code{domain_scope}, and \code{ordered_query} are
+#' pre-specified by this wrapper function. Users must set
+#' \code{organism = "mmusculus"} when working on mouse data.
 #' @references Kolberg, L. et al, 2020 and Raudvere, U. et al, 2019
-#' @return A list object with the following entries
+#' @return A list object where each element is a result list for a group. Each
+#' result list contains two elements:
 #' \item{result}{data.frame of main GO analysis result.}
 #' \item{meta}{Meta information for the query.}
 #'
 #' See \code{\link[gprofiler2]{gost}}. for detailed explanation.
 #' @export
 #' @examples
-#' res <- runWilcoxon(pbmcPlot)
+#' res <- runMarkerDEG(pbmcPlot)
 #' # Setting `significant = FALSE` because it's hard for a gene list obtained
 #' # from small test dataset to represent real-life biology.
 #' go <- runGOEnrich(res, group = 0, significant = FALSE)
 runGOEnrich <- function(
         result,
         group = NULL,
         useBg = TRUE,
-        orderBy = "logFC",
+        orderBy = "padj",
         logFCThresh = 1,
-        padjThresh = 0.01,
+        padjThresh = 0.05,
+        splitReg = FALSE,
         ...
 ) {
     if (!requireNamespace("gprofiler2", quietly = TRUE))
@@ -164,9 +176,12 @@ runGOEnrich <- function(
         domain_scope <- "custom"
     }
     filter <- result$group %in% group &
-        result$logFC > logFCThresh &
+        abs(result$logFC) > logFCThresh &
         result$padj < padjThresh
+    filter[is.na(filter)] <- FALSE
     result <- result[filter, ]
+    resultUp <- result[result$logFC > 0,]
+    if (isTRUE(splitReg)) resultDown <- result[result$logFC < 0,]
 
     ordered_query <- FALSE
     if (!is.null(orderBy)) {
@@ -176,18 +191,25 @@ runGOEnrich <- function(
             stop("`orderBy` should be one of 'logFC', 'pval' or 'padj'.")
         }
         if (orderBy == "logFC") {
-            result <- result[order(result$logFC, decreasing = TRUE),]
+            resultUp <- resultUp[order(resultUp$logFC, decreasing = TRUE),]
+            if (isTRUE(splitReg))
+                resultDown <- resultDown[order(resultDown$logFC),]
         } else {
-            result <- result[order(result[[orderBy]], decreasing = FALSE),]
+            resultUp <- resultUp[order(resultUp[[orderBy]]),]
+            if (isTRUE(splitReg))
+                resultDown <- resultDown[order(resultDown[[orderBy]]),]
         }
     }
+    resultList <- list()
+    for (g in unique(result$group)) {
+        query <- list(Up = resultUp$feature[resultUp$group == g])
+        if (splitReg) query$Down <- resultDown$feature[resultDown$group == g]
+        resultList[[g]] <- gprofiler2::gost(
+            query = query, custom_bg = bg, domain_scope = domain_scope,
+            ordered_query = ordered_query, ...
+        )
+    }
 
-    gsLists <- split(result$feature, droplevels(result$group))
-    output <- gprofiler2::gost(
-        query = gsLists, custom_bg = bg, domain_scope = domain_scope,
-        ordered_query = ordered_query, ...
-    )
-
-    return(output)
+    return(resultList)
 }