seurat_obj <- SetupForWGCNA(scRNA,
  gene_select = "custom", # the gene selection approach
  gene_list = features,  ## custom genes were used total 19996 genes used        
  wgcna_name = "hdWGCNA") # the name of the hdWGCNA experiment

  seurat_obj <- MetacellsByGroups(seurat_obj = seurat_obj,
  group.by = c("seurat_clusters"), # specify the columns in seurat_obj@meta.data to group by
  assay = 'MAGIC_RNA',
  slot = 'data',
  reduction = 'pca', # select the dimensionality reduction to perform KNN on
  k = 25, # nearest-neighbors parameter
  max_shared = 10, # maximum number of shared cells between two metacells
  ident.group = 'seurat_clusters') # set the Idents of the metacell seurat object

  # normalize metacell expression matrix:
seurat_obj <- NormalizeMetacells(seurat_obj, verbose=FALSE)

seurat_obj <- SetDatExpr(seurat_obj,
  group_name = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16"), # the name of the group of interest in the group.by column
  group.by='seurat_clusters', # the metadata column containing the cell type info. This same column should have also been used in MetacellsByGroups
  assay = 'MAGIC_RNA', # using RNA assay
  slot = 'data') # using normalized data

  # Test different soft powers:
seurat_obj <- TestSoftPowers(seurat_obj,
  networkType = 'unsigned') # 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)

## Construct co-expression network

seurat_obj <- ConstructNetwork(seurat_obj, soft_power=4,
  setDatExpr=FALSE,
  networkType = 'unsigned',
  TOMType = "unsigned",
  tom_name = 'TOM') # name of the topoligical overlap matrix written to disk

 # 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
seurat_obj <- ModuleEigengenes(seurat_obj,
 group.by.vars="seurat_clusters")

 # compute eigengene-based connectivity (kME):
seurat_obj <- ModuleConnectivity(seurat_obj,
  group.by = 'seurat_clusters', group_name = '2',sparse=FALSE)

 # get the module assignment table:
modules <- GetModules(seurat_obj)

# show the first 6 columns:
head(modules[,1:6])
dim(modules)
#write.csv(modules, "Modules.csv", row.names=FALSE)

seurat_obj <- RunModuleUMAP(
  seurat_obj,
  n_hubs = 10, # number of hub genes to include for the UMAP embedding
  n_neighbors=15, # neighbors parameter for UMAP
  min_dist=0.1) # min distance between points in UMAP space