#' @title Gene set enrichment
 #' @description Identify and return significantly-enriched terms for each gene
 #'  module in a Celda object. Performs gene set enrichment analysis for Celda
 #'  identified modules using the enrichR package.
 #' @author Ahmed Youssef
 #' @param counts Integer count matrix. Rows represent genes and columns
 #'  represent cells. Row names of the matrix should be gene names.
 #' @param celdaModel Celda object of class `celda_G` or `celda_CG`.
 #' @param databases Character vector. Name of reference database. Available
 #'  databases can be viewed by running \code{enrichR::listEnrichrDbs()}.
 #' @param fdr False discovery rate (FDR). Numeric. Cutoff value for adjusted
 #'  p-value, terms with FDR below this value are considered significantly
 #'  enriched.
 #' @return List of length 'L' where each member contains the significantly
 #'  enriched terms for the corresponding module.

 #' @examples

 #' library(M3DExampleData)
 #' counts <- M3DExampleData::Mmus_example_list$data
 #' #subset 100 genes for fast clustering

 #' counts <- counts[1500:2000, ]

 #' #cluster genes into 10 modules for quick demo

 #' cm <- celda_G(counts = as.matrix(counts), L = 10, verbose = FALSE)

 #' gse <- geneSetEnrich(counts,
 #'     cm,
 #'     databases = c('GO_Biological_Process_2018','GO_Molecular_Function_2018'))

 #' @importFrom enrichR enrichr

 #' @importFrom enrichR listEnrichrDbs

 #' @export
 geneSetEnrich <- function(counts, celdaModel, databases, fdr = 0.05) {
     #check for correct celda object
     if (!(class(celdaModel) %in% c("celda_G", "celda_CG"))) {
         stop("No gene modules in celda object. ",
             "Please provide object of class celda_G or celda_CG.")
     }
 
     #initialize list with one entry for each gene module

     modules <- vector("list", length = params(celdaModel)$L)

 
     #create dataframe with gene-module associations

     genes <- data.frame(gene = rownames(counts),
         module = clusters(celdaModel)$y)

 
     #iterate over each module, get genes in that module, add to list

     for (i in seq_len(params(celdaModel)$L)) {

         modules[[i]] <- as.character(genes[genes$module == i, "gene"])

     }
 
     #enrichment analysis
     enrichment <- lapply(modules, function(module) {

         invisible(utils::capture.output(table <- enrichR::enrichr(
             genes = module,

             databases = databases)))
         table <- Reduce(f = rbind, x = table)

         table[table$Adjusted.P.value < fdr, "Term"]

     })
 
     #return results as a list
     return(enrichment)
 }