Is there a way to perform Differential ME analysis using two independent seurat objects?

[ozambadam]

I'm comparing a Setaria Viridis WT vs a mutant (one gene knocked-down) to look for the changes on the coexpression network caused by the lack of this specific gene. I have been using many functions of hdWGCNA like projecting modules and module preservation analysis, but I think DME analysis could be one of the most interesting features to explore. However, given that I have been using the other functions my datasets are separated in seurat_ref and seurat_query. Given this, and aiming to follow the same line in this analysis, I would like to know if it is possible to apply de DME function to these separated objects or if it is necessary to create a combinatory one where I have both the mutant and the wildtype cells in order to perform this analysis.

[smorabit]

Hi,

Thanks for your interest in hdWGCNA. To clarify, you can only run DME analysis on a single Seurat object, so you will need to combine them. However, merging and subsetting objects after running hdWGCNA is not supported, so you will run into more issues. I recommend starting your analysis with a merge object. If you are using Seurat v5, you will need to run the function JoinLayers in order to get a merged expression matrix.

You can still identify modules in your WT vs Mutant groups separately by specifying the expression matrix in the SetDatExpr function, for example, you code might look something like this:

seurat_obj <- SetupForWGCNA(
  seurat_obj,
  gene_select = "fraction", 
  fraction = 0.05, 
  wgcna_name = "mutant" 
)

# construct metacells  in each group
seurat_obj <- MetacellsByGroups(
  seurat_obj = seurat_obj,
  group.by = c("cell_type", "Sample", "genotype"), 
  reduction = 'harmony', 
  k = 25, 
  max_shared = 10, 
  ident.group = 'cell_type'
)

# normalize metacell expression matrix:
seurat_obj <- NormalizeMetacells(seurat_obj)

# set up expression matrix
seurat_obj <- SetDatExpr(
  seurat_obj,
  group_name = "mutant", 
  group.by="genotype"
)

# construct co-expression network:
seurat_obj <- TestSoftPowers(seurat_obj)
seurat_obj <- ConstructNetwork(seurat_obj)

# need to run ScaleData first or else harmony throws an error:
#seurat_obj <- ScaleData(seurat_obj, features=VariableFeatures(seurat_obj))

# compute all MEs in the full single-cell dataset, and compute connectivity
seurat_obj <- ModuleEigengenes(seurat_obj)
seurat_obj <- ModuleConnectivity(
  seurat_obj,
  group.by = 'genotype', group_name = 'mutant'
)

You can then run hdWGCNA again for the WT group.

# get genes used for the other analysis
genes_use <- GetWGCNAGenes(seurat_obj)

# setup another WGCNA experiment
seurat_obj <- SetupForWGCNA(
  seurat_obj,
  gene_select = "custom",
  gene_list = genes_use,
  metacell_location = 'mutant', # this allows us to use the same metacell object, so we don't need to re-make it.
  wgcna_name = "WT" 
)

# set up expression matrix
seurat_obj <- SetDatExpr(
  seurat_obj,
  group_name = "WT", 
  group.by="genotype"
)

# construct co-expression network:
seurat_obj <- TestSoftPowers(seurat_obj)
seurat_obj <- ConstructNetwork(seurat_obj)

# need to run ScaleData first or else harmony throws an error:
#seurat_obj <- ScaleData(seurat_obj, features=VariableFeatures(seurat_obj))

# compute all MEs in the full single-cell dataset, and compute connectivity
seurat_obj <- ModuleEigengenes(seurat_obj)
seurat_obj <- ModuleConnectivity(
  seurat_obj,
  group.by = 'genotype', group_name = 'WT'
)

This will give you a set of modules from the mutant group and the WT group separately, but you will have gene expression scores (MEs) for the modules in the entire merged dataset, enabling you to do DME analysis downstream.

In this case where the same seurat object has more than one WGCNA experiment, when you run FindDMEs, make sure that you specify the correct wgcna_name.

group1 <- seurat_obj@meta.data %>% subset(genotype == 'mutant') %>% rownames
group2 <- seurat_obj@meta.data %>% subset(genotype == 'WT) %>% rownames

DMEs_mutant_modules <- FindDMEs(
  seurat_obj,
  barcodes1 = group1,
  barcodes2 = group2,
  test.use='wilcox',
  wgcna_name='mutant'
)

DMEs_WT_modules <- FindDMEs(
  seurat_obj,
  barcodes1 = group1,
  barcodes2 = group2,
  test.use='wilcox',
  wgcna_name='WT'
)

Let me know if this answers your question, or if you have any doubts.

[ozambadam]

Thanks for your descriptive answer. I´m trying this in my dataset. I´ll let you know if I have any doubt.