@@ -20,14 +20,20 @@

 #' @param fdr False discovery rate (FDR). Numeric. Cutoff value for adjusted

 #'  p-value, terms with FDR below this value are considered significantly

 #'  enriched.

-#' @param ... Ignored. Placeholder to prevent check warning.

 #' @return List of length 'L' where each member contains the significantly

 #'  enriched terms for the corresponding module.

 #' @importFrom enrichR enrichr

 #' @importFrom enrichR listEnrichrDbs

 #' @export

-setGeneric("geneSetEnrich", function(x, ...) {

-    standardGeneric("geneSetEnrich")})

+setGeneric("geneSetEnrich",

+    function(x,

+        celdaModel,

+        useAssay = "counts",

+        altExpName = "featureSubset",

+        databases,

+        fdr = 0.05) {

+

+        standardGeneric("geneSetEnrich")})

 #' @rdname geneSetEnrich

@@ -9,10 +9,10 @@

 #'  with the matrix located in the assay slot under \code{useAssay}.

 #'  Rows represent features and columns represent cells. Rownames of the

 #'  matrix or \linkS4class{SingleCellExperiment} object should be gene names.

-#' @param useAssay A string specifying which \link[SummarizedExperiment]{assay}

+#' @param useAssay A string specifying which \link{assay}

 #'  slot to use if \code{x} is a

 #'  \linkS4class{SingleCellExperiment} object. Default "counts".

-#' @param altExpName The name for the \link[SingleCellExperiment]{altExp} slot

+#' @param altExpName The name for the \link{altExp} slot

 #'  to use. Default "featureSubset".

 #' @param celdaModel Celda object of class \code{celda_G} or \code{celda_CG}.

 #' @param databases Character vector. Name of reference database. Available

@@ -20,6 +20,7 @@

 #' @param fdr False discovery rate (FDR). Numeric. Cutoff value for adjusted

 #'  p-value, terms with FDR below this value are considered significantly

 #'  enriched.

+#' @param ... Ignored. Placeholder to prevent check warning.

 #' @return List of length 'L' where each member contains the significantly

 #'  enriched terms for the corresponding module.

 #' @importFrom enrichR enrichr

@@ -55,7 +55,8 @@ setMethod("geneSetEnrich",

             length = S4Vectors::metadata(altExp)$celda_parameters$L)

         # create dataframe with gene-module associations

-        genes <- data.frame(gene = rownames(altExp), module = celdaModules(x))

+        genes <- data.frame(gene = rownames(altExp),

+            module = celdaModules(x, altExpName = altExpName))

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

         for (i in seq_len(S4Vectors::metadata(altExp)$celda_parameters$L)) {

@@ -12,6 +12,8 @@

 #' @param useAssay A string specifying which \link[SummarizedExperiment]{assay}

 #'  slot to use if \code{x} is a

 #'  \linkS4class{SingleCellExperiment} object. Default "counts".

+#' @param altExpName The name for the \link[SingleCellExperiment]{altExp} slot

+#'  to use. Default "featureSubset".

 #' @param celdaModel Celda object of class \code{celda_G} or \code{celda_CG}.

 #' @param databases Character vector. Name of reference database. Available

 #'  databases can be viewed by \link[enrichR]{listEnrichrDbs}.

@@ -40,17 +42,23 @@ setGeneric("geneSetEnrich", function(x, ...) {

 #' @export

 setMethod("geneSetEnrich",

     signature(x = "SingleCellExperiment"),

-    function(x, useAssay = "counts", databases, fdr = 0.05) {

+    function(x,

+        useAssay = "counts",

+        altExpName = "featureSubset",

+        databases,

+        fdr = 0.05) {

+

+        altExp <- SingleCellExperiment::altExp(x, e = altExpName)

         # initialize list with one entry for each gene module

         modules <- vector("list",

-            length = S4Vectors::metadata(x)$celda_parameters$L)

+            length = S4Vectors::metadata(altExp)$celda_parameters$L)

         # create dataframe with gene-module associations

-        genes <- data.frame(gene = rownames(x), module = celdaModules(x))

+        genes <- data.frame(gene = rownames(altExp), module = celdaModules(x))

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

-        for (i in seq_len(S4Vectors::metadata(x)$celda_parameters$L)) {

+        for (i in seq_len(S4Vectors::metadata(altExp)$celda_parameters$L)) {

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

         }

@@ -23,7 +23,8 @@

 #' @importFrom enrichR enrichr

 #' @importFrom enrichR listEnrichrDbs

 #' @export

-setGeneric("geneSetEnrich", function(x, ...) {standardGeneric("geneSetEnrich")})

+setGeneric("geneSetEnrich", function(x, ...) {

+    standardGeneric("geneSetEnrich")})

 #' @rdname geneSetEnrich

@@ -33,7 +33,7 @@ setGeneric("geneSetEnrich", function(x, ...) {standardGeneric("geneSetEnrich")})

 #' # subset 500 genes for fast clustering

 #' counts <- counts[seq(1501, 2000), ]

 #' # cluster genes into 10 modules for quick demo

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

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

 #' gse <- geneSetEnrich(sce,

 #'   databases = c("GO_Biological_Process_2018", "GO_Molecular_Function_2018"))

 #' @export

@@ -1,65 +1,109 @@

 #' @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`.

+#'  module in a Celda object or a \linkS4class{SingleCellExperiment} object.

+#'  Performs gene set enrichment analysis for Celda

+#'  identified modules using the \link[enrichR]{enrichr}.

+#' @author Ahmed Youssef, Zhe Wang

+#' @param x A numeric \link{matrix} of counts or a

+#'  \linkS4class{SingleCellExperiment}

+#'  with the matrix located in the assay slot under \code{useAssay}.

+#'  Rows represent features and columns represent cells. Rownames of the

+#'  matrix or \linkS4class{SingleCellExperiment} object should be gene names.

+#' @param useAssay A string specifying which \link[SummarizedExperiment]{assay}

+#'  slot to use if \code{x} is a

+#'  \linkS4class{SingleCellExperiment} object. Default "counts".

+#' @param celdaModel Celda object of class \code{celda_G} or \code{celda_CG}.

 #' @param databases Character vector. Name of reference database. Available

-#'  databases can be viewed by running \code{enrichR::listEnrichrDbs()}.

+#'  databases can be viewed by \link[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.

+#' @importFrom enrichR enrichr

+#' @importFrom enrichR listEnrichrDbs

+#' @export

+setGeneric("geneSetEnrich", function(x, ...) {standardGeneric("geneSetEnrich")})

+

+

+#' @rdname geneSetEnrich

 #' @examples

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

-#' # subset 100 genes for fast clustering

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

+#' # subset 500 genes for fast clustering

+#' counts <- counts[seq(1501, 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

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

+#' gse <- geneSetEnrich(sce,

+#'   databases = c("GO_Biological_Process_2018", "GO_Molecular_Function_2018"))

+#' @export

+setMethod("geneSetEnrich",

+    signature(x = "SingleCellExperiment"),

+    function(x, useAssay = "counts", databases, fdr = 0.05) {

+

+        # initialize list with one entry for each gene module

+        modules <- vector("list",

+            length = S4Vectors::metadata(x)$celda_parameters$L)

+

+        # create dataframe with gene-module associations

+        genes <- data.frame(gene = rownames(x), module = celdaModules(x))

+

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

+        for (i in seq_len(S4Vectors::metadata(x)$celda_parameters$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)

+    }

+)

+

+

+#' @rdname geneSetEnrich

 #' @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."

-    )

-  }

+setMethod("geneSetEnrich",

+    signature(x = "matrix"),

+    function(x, 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)

+        # 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

-  )

+        # create dataframe with gene-module associations

+        genes <- data.frame(gene = rownames(x),

+            module = celdaClusters(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"])

-  }

+        # 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"]

-  })

+        # 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)

-}

+        # return results as a list

+        return(enrichment)

+    }

+)

@@ -16,44 +16,50 @@

 #' @examples

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

-#' #subset 100 genes for fast clustering

+#' # subset 100 genes for fast clustering

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

-#' #cluster genes into 10 modules for quick demo

+#' # 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'))

+#'   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.")

-    }

+  # 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)

+  # 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)

+  # 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"])

-    }

+  # 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"]

-    })

+  # 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)

+  # return results as a list

+  return(enrichment)

 }

@@ -34,13 +34,14 @@ geneSetEnrich <- function(counts, celdaModel, databases, fdr = 0.05) {

     }

     #initialize list with one entry for each gene module

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

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

     #create dataframe with gene-module associations

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

+    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(celdaModel@params$L)) {

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

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

     }

@@ -24,6 +24,7 @@

 #'     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

@@ -23,6 +23,7 @@

 #' gse <- geneSetEnrich(counts,

 #'     cm,

 #'     databases = c('GO_Biological_Process_2018','GO_Molecular_Function_2018'))

+#' @importFrom enrichR enrichr

 #' @export

 geneSetEnrich <- function(counts, celdaModel, databases, fdr = 0.05) {

     #check for correct celda object

@@ -17,7 +17,7 @@

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' #subset 100 genes for fast clustering

-#' counts <- counts[sample(seq_len(nrow(counts)), size = 100), ]

+#' 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,

@@ -17,7 +17,7 @@

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' #subset 100 genes for fast clustering

-#' counts <- counts[seq(1200, 2000), ]

+#' counts <- counts[sample(seq_len(nrow(counts)), size = 100), ]

 #' #cluster genes into 10 modules for quick demo

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

 #' gse <- geneSetEnrich(counts,

@@ -39,15 +39,16 @@ geneSetEnrich <- function(counts, celdaModel, databases, fdr = 0.05) {

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

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

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

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

     }

     #enrichment analysis

     enrichment <- lapply(modules, function(module) {

-        invisible(capture.output(table <- enrichR::enrichr(genes = 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']

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

     })

     #return results as a list

@@ -15,7 +15,6 @@

 #'  enriched terms for the corresponding module.

 #' @examples

 #' library(M3DExampleData)

-#' library(celda)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' #subset 100 genes for fast clustering

 #' counts <- counts[seq(1200, 2000), ]

@@ -18,8 +18,7 @@

 #' library(celda)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' #subset 100 genes for fast clustering

-#' counts <- counts[sample(seq_len(nrow(counts)),

-#'     size = 100),]

+#' counts <- counts[seq(1200, 2000), ]

 #' #cluster genes into 10 modules for quick demo

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

 #' gse <- geneSetEnrich(counts,

@@ -15,6 +15,7 @@

 #'  enriched terms for the corresponding module.

 #' @examples

 #' library(M3DExampleData)

+#' library(celda)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' #subset 100 genes for fast clustering

 #' counts <- counts[sample(seq_len(nrow(counts)),

@@ -14,8 +14,6 @@

 #' @return List of length 'L' where each member contains the significantly

 #'  enriched terms for the corresponding module.

 #' @examples

-#' if (!requireNamespace("BiocManager", quietly = TRUE))

-#'     install.packages("M3DExampleData")

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' #subset 100 genes for fast clustering

@@ -19,10 +19,10 @@

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' #subset 100 genes for fast clustering

-#' counts <- counts[sample(seq_len(nrow(countsMatrix)),

+#' counts <- counts[sample(seq_len(nrow(counts)),

 #'     size = 100),]

 #' #cluster genes into 10 modules for quick demo

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

+#' 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'))

@@ -19,7 +19,7 @@

 #' library(M3DExampleData)

 #' counts <- M3DExampleData::Mmus_example_list$data

 #' #subset 100 genes for fast clustering

-#' counts <- countsMatrix[sample(seq_len(nrow(countsMatrix)),

+#' counts <- counts[sample(seq_len(nrow(countsMatrix)),

 #'     size = 100),]

 #' #cluster genes into 10 modules for quick demo

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

@@ -13,7 +13,7 @@

 #'  enriched.

 #' @return List of length 'L' where each member contains the significantly

 #'  enriched terms for the corresponding module.

-#' @example

+#' @examples

 #' if (!requireNamespace("BiocManager", quietly = TRUE))

 #'     install.packages("M3DExampleData")

 #' library(M3DExampleData)

new file mode 100644

@@ -0,0 +1,58 @@

+#' @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.

+#' @example

+#' if (!requireNamespace("BiocManager", quietly = TRUE))

+#'     install.packages("M3DExampleData")

+#' library(M3DExampleData)

+#' counts <- M3DExampleData::Mmus_example_list$data

+#' #subset 100 genes for fast clustering

+#' counts <- countsMatrix[sample(seq_len(nrow(countsMatrix)),

+#'     size = 100),]

+#' #cluster genes into 10 modules for quick demo

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

+#' gse <- geneSetEnrich(counts,

+#'     cm,

+#'     databases = c('GO_Biological_Process_2018','GO_Molecular_Function_2018'))

+#' @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 = celdaModel@params$L)

+

+    #create dataframe with gene-module associations

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

+

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

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

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

+    }

+

+    #enrichment analysis

+    enrichment <- lapply(modules, function(module) {

+        invisible(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)

+}