bugfix for ModulEigengenes with Seurat v5 assay

R/ModuleConnectivity.R

@@ -64,6 +64,7 @@ ModuleConnectivity <- function(
   genes_use <- as.character(modules$gene_name)
   params <- GetWGCNAParams(seurat_obj, wgcna_name)
 
+  # get the assay
   if(is.null(assay)){assay <- DefaultAssay(seurat_obj)}
 
   if(!is.null(group.by)){
@@ -79,10 +80,8 @@ ModuleConnectivity <- function(
       seurat_obj, assay=assay, layer=layer
     )[genes_use,cells.use]
   } else{
-    exp_mat <- GetAssayData(
-      seurat_obj,
-      assay=assay,
-      slot=slot
+    exp_mat <- Seurat::GetAssayData(
+      seurat_obj, assay=assay, slot=slot
     )[genes_use,cells.use]
   }
 

R/ModuleEigengenes.R

@@ -32,25 +32,35 @@ ComputeModuleEigengene <- function(
 
   # get the assay
   if(is.null(assay)){assay <- DefaultAssay(seurat_obj)}
-  if(dim(seurat_obj@assays[[assay]]@data)[1] == 0){
-    stop(paste0("Normalized data slot not found in selected assay ", assay))
-  }
 
   # get genes in this module:
   cur_genes <- modules %>% subset(module == cur_mod) %>% .$gene_name %>% as.character()
 
-  # subset seurat object by these genes only:
-  X_dat <- GetAssayData(seurat_obj, slot='data', assay = assay)[cur_genes,]
-  if(dim(seurat_obj@assays[[assay]]@counts)[1] == 0){
-    X <- X_dat
+  # get the expression matrix:
+  if(CheckSeurat5()){
+    X <- SeuratObject::LayerData(seurat_obj, layer='counts', assay=assay)[cur_genes,]
+    X_dat <- SeuratObject::LayerData(seurat_obj, layer='data', assay=assay)[cur_genes,]
   } else{
-    X <- GetAssayData(seurat_obj, slot='counts', assay = assay)[cur_genes,]
+    X <- Seurat::GetAssayData(seurat_obj, slot='counts', assay=assay)[cur_genes,]
+    X_dat <- Seurat::GetAssayData(seurat_obj, slot='data', assay=assay)[cur_genes,]
   }
 
+  # subset seurat object by these genes only:
+  # X_dat <- GetAssayData(seurat_obj, slot='data', assay = assay)[cur_genes,]
+  # if(dim(seurat_obj@assays[[assay]]@counts)[1] == 0){
+  #   X <- X_dat
+  # } else{
+  #   X <- GetAssayData(seurat_obj, slot='counts', assay = assay)[cur_genes,]
+  # }
+
   # create seurat obj with just these genes
   cur_seurat <- CreateSeuratObject(X, assay = assay, meta.data = seurat_obj@meta.data)
-  cur_seurat <- SetAssayData(cur_seurat, slot='data', new.data=X_dat, assay=assay)
-
+
+  if(CheckSeurat5()){
+    cur_seurat <- SetAssayData(cur_seurat, layer='data', new.data=X_dat, assay=assay)
+  } else{
+    cur_seurat <- SetAssayData(cur_seurat, slot='data', new.data=X_dat, assay=assay)
+  }
   # scale the subsetted expression dataset:
   if(is.null(vars.to.regress)){
     cur_seurat <- ScaleData(cur_seurat, features=rownames(cur_seurat), model.use=scale.model.use)
@@ -61,7 +71,11 @@ ComputeModuleEigengene <- function(
   }
 
   # compute average expression of each gene
-  cur_expr <- GetAssayData(cur_seurat, slot='data')
+  if(CheckSeurat5()){
+    cur_expr <- SeuratObject::GetAssayData(cur_seurat, layer='data')
+  } else{
+    cur_expr <- Seurat::GetAssayData(cur_seurat, slot='data')
+  }
   expr <- Matrix::t(cur_expr)
   averExpr <- Matrix::rowSums(expr) / ncol(expr)
 
@@ -152,7 +166,6 @@ ComputeModuleEigengene <- function(
 #' @param scale.model.use model to scale data when running ScaleData choices are "linear", "poisson", or "negbinom"
 #' @param pc_dim Which PC to use as the module eigengene? Default to 1.
 #' @param assay Assay in seurat_obj to compute module eigengenes. Default is DefaultAssay(seurat_obj)
-#' @param exclude_grey logical determining whether to compute MEs for the grey module
 #' @param wgcna_name name of the WGCNA experiment
 #'
 #' @details
@@ -179,7 +192,6 @@ ModuleEigengenes <- function(
   verbose=TRUE,
   assay = NULL,
   pc_dim = 1,
-  exclude_grey = FALSE,
   wgcna_name=NULL, ...
 ){
 
@@ -197,14 +209,6 @@ ModuleEigengenes <- function(
   # get the assay
   if(is.null(assay)){assay <- DefaultAssay(seurat_obj)}
 
-  # check for the data slot in this assay
-  if(dim(seurat_obj@assays[[assay]]@data)[1] == 0){
-    stop(paste0("Normalized data slot not found in selected assay ", assay))
-  }
-
-  # exclude grey doesn't work yet:
-  if(exclude_grey){exclude_grey <- FALSE}
-
   me_list <- list()
   harmonized_me_list <- list()
 
@@ -236,11 +240,6 @@ ModuleEigengenes <- function(
   mods <- levels(modules$module)
   mods_loop <- mods
 
-  # exclude grey:
-  if(exclude_grey){
-    mods_loop <- mods[mods != 'grey']
-  }
-
   # loop over modules:
   for(cur_mod in mods_loop){
 

R/SelectNetworkGenes.R

@@ -6,8 +6,10 @@
 #' @param seurat_obj A Seurat object
 #' @param type How to select genes? Select "variable", "fraction", "all", or "custom".
 #' @param fraction A numeric that determines the minimum cells that a gene must be expressed in order to be included. For example, fraction = 0.05 means that 5% of cells must express a gene (count > 0) for it to be included.
+#' @param assay Assay in seurat_obj to compute module eigengenes. Default is DefaultAssay(seurat_obj)
 #' @param gene_list A character string of gene names, only used if type = "custom"
-#' 
+#' @param wgcna_name name of the WGCNA experiment
+#'
 #' @details 
 #' SelectNetworkGenes allows us to specify the genes that will be used for co-expression network analysis.
 #' This function is called by SetupForWGCNA. By default, the variable features in VariableFeatures(seurat_obj) are used.
@@ -17,35 +19,38 @@
 #' For example, by setting fraction=0.1 and group.by='seurat_clusters', this function will identify the set of genes 
 #' that are expressed in 10% of cells in at least one of the clusters.
 #' 
-#' 
 #' @export
 SelectNetworkGenes <- function(
   seurat_obj,
   gene_select="variable", fraction=0.05,
   group.by=NULL, # should be a column in the Seurat object, eg clusters
   gene_list=NULL,
+  assay = NULL,
   wgcna_name=NULL
  ){
 
+  # set as active assay if wgcna_name is not given
+  if(is.null(wgcna_name)){wgcna_name <- seurat_obj@misc$active_wgcna}
+
   # validate inputs:
   if(!(gene_select %in% c("variable", "fraction", "all", "custom"))){
     stop(paste0("Invalid selection gene_select: ", gene_select, '. Valid gene_selects are variable, fraction, all, or custom.'))
   }
 
-
-  # get assay
-  assay <- DefaultAssay(seurat_obj)
-  tryCatch(
-   tmp <- dim(seurat_obj[[assay]]@counts),
-   error = function(e){
-     print("Assay must contain counts slot.")
-   })
+  # get the assay
+  if(is.null(assay)){assay <- DefaultAssay(seurat_obj)}
+
+  # get the expression matrix
+  if(CheckSeurat5()){
+    expr_mat <- SeuratObject::LayerData(seurat_obj, layer='counts', assay=assay)
+  } else{
+    expr_mat <- Seurat::GetAssayData(seurat_obj, slot='counts', assay=assay)
+  }
 
   # handle different selection strategies
   if(gene_select == "fraction"){
 
     # binarize counts matrix in chunks to save memory
-    expr_mat <- GetAssayData(seurat_obj, slot='counts')
     n_chunks <- ceiling(ncol(expr_mat) / 10000)
 
     if(n_chunks == 1){

R/getters_and_setters.R

@@ -350,6 +350,7 @@ GetDatExpr <- function(seurat_obj, wgcna_name=NULL){
 #' @param multi_groups A character vecrtor containing the names of groups to select
 #' @param assay The name of the assay in the Seurat object
 #' @param slot The name of the slot in the Seurat object (counts, data)
+#' @param layer Layer to extract data for aggregation. Default = 'counts'. Layer is used with Seurat v5 instead of slot.
 #' @param mat A Matrix containing gene expression data. Supplying a matrix using this parameter ignores other options. This is almost exclusively used for pseudobulk analysis.
 #' @param wgcna_name A string containing the name of the WGCNA slot in seurat_obj@misc. Default = NULL which retrieves the currently active WGCNA data
 #' @keywords scRNA-seq
@@ -363,6 +364,7 @@ SetMultiExpr <- function(
   multi_groups = NULL,
   assay=NULL,
   slot='data',
+  layer = 'data',
   mat=NULL,
   mat_group_delim=3,
   wgcna_name=NULL,
@@ -416,7 +418,8 @@ SetMultiExpr <- function(
         use_metacells = use_metacells,
         wgcna_name = wgcna_name,
         assay = assay,
-        slot = slot
+        slot = slot,
+        layer = layer
       ) 
       as.matrix(cur_datExpr)
     })