pca error in MetacellsByGroups

[ShaolongCao]

Hi Sam,

I am interested in using your hdWGCNA package for snRNAseq analysis. However, when I tried your data and code in tutorial. I got this error for the MetacellsByGroups function: "Error in MetacellsByGroups(seurat_obj = seurat_obj, group.by = c("cell_type", : Invalid reduction (pca). Reductions in Seurat object: harmony, umap"

Here is the code I used:

load the Zhou et al snRNA-seq dataset

seurat_obj <- readRDS('Zhou_2020.rds')

seurat_obj <- SetupForWGCNA( seurat_obj, gene_select = "fraction", # the gene selection approach fraction = 0.05, # fraction of cells that a gene needs to be expressed in order to be included wgcna_name = "tutorial" # the name of the hdWGCNA experiment )

construct metacells in each group

seurat_obj <- MetacellsByGroups( seurat_obj = seurat_obj, group.by = c("cell_type", "Sample"), # specify the columns in seurat_obj@meta.data to group by k = 25, # nearest-neighbors parameter max_shared = 10, # maximum number of shared cells between two metacells ident.group = 'cell_type' # set the Idents of the metacell seurat object )

Do you have any idea what was happened?

Thanks, Shaolong

[smorabit]

Hi Shaolong,

The error message that you got is essentially telling you what to do.

"Invalid reduction (pca). Reductions in Seurat object: harmony, umap"

So you either have to run PCA on your dataset with RunPCA, or you can just set the reduction to harmony reduction="harmony" inside of the MetacellsByGroups function.

Let me know if this solves your issue.

[KangchengHou]

Related to this thread, I am noting that the output from the tutorial is different from the run I had in the following:

# single-cell analysis package
library(Seurat)

# plotting and data science packages
library(tidyverse)
library(cowplot)
library(patchwork)

# co-expression network analysis packages:
library(WGCNA)
library(hdWGCNA)

# using the cowplot theme for ggplot
theme_set(theme_cowplot())

# set random seed for reproducibility
set.seed(12345)

# optionally enable multithreading
enableWGCNAThreads(nThreads = 8)

# load the Zhou et al snRNA-seq dataset
seurat_obj <- readRDS('data/Zhou_2020.rds')

seurat_obj <- seurat_obj %>%
  NormalizeData() %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA()

p <- DimPlot(seurat_obj, group.by='cell_type', label=TRUE) +
   umap_theme() + ggtitle('Zhou et al Control Cortex') + NoLegend()

p

seurat_obj <- SetupForWGCNA(
  seurat_obj,
  gene_select = "fraction", # the gene selection approach
  fraction = 0.05, # fraction of cells that a gene needs to be expressed in order to be included
  wgcna_name = "tutorial" # the name of the hdWGCNA experiment
)

# construct metacells  in each group
seurat_obj <- MetacellsByGroups(
  seurat_obj = seurat_obj,
  group.by = c("cell_type", "Sample"), # specify the columns in seurat_obj@meta.data to group by
  k = 25, # nearest-neighbors parameter
  max_shared = 10, # maximum number of shared cells between two metacells
  ident.group = 'cell_type', # set the Idents of the metacell seurat object
)

seurat_obj <- NormalizeMetacells(seurat_obj)
seurat_obj <- ScaleMetacells(seurat_obj, features=VariableFeatures(seurat_obj))
seurat_obj <- RunPCAMetacells(seurat_obj, features=VariableFeatures(seurat_obj))
seurat_obj <- RunHarmonyMetacells(seurat_obj, group.by.vars='Sample')
seurat_obj <- RunUMAPMetacells(seurat_obj, reduction='harmony', dims=1:15)

p1 <- DimPlotMetacells(seurat_obj, group.by='cell_type') + umap_theme() + ggtitle("Cell Type")
p2 <- DimPlotMetacells(seurat_obj, group.by='Sample') + umap_theme() + ggtitle("Sample")

p1 | p2

The above code gives the following figure,

And

seurat_obj <- SetDatExpr(
  seurat_obj,
  group_name = "INH", # the name of the group of interest in the group.by column
  group.by='cell_type', # the metadata column containing the cell type info. This same column should have also been used in MetacellsByGroups
  assay = 'RNA', # using RNA assay
  slot = 'data' # using normalized data
)

# Test different soft powers:
seurat_obj <- TestSoftPowers(
  seurat_obj,
  networkType = 'signed' # you can also use "unsigned" or "signed hybrid"
)

# plot the results:
plot_list <- PlotSoftPowers(seurat_obj)

# assemble with patchwork
wrap_plots(plot_list, ncol=2)

gives

Thanks in advance!