🤩 hdWGCNA | 单细胞数据怎么做WGCNA呢！？（一）（小试牛刀）

rm(list = ls())
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(WGCNA)
library(hdWGCNA)

theme_set(theme_cowplot())
set.seed(123)
enableWGCNAThreads(nThreads = 8)

seurat_obj <- readRDS('Zhou_2020.rds')

p <- DimPlot(seurat_obj, 
             group.by='cell_type',
             label=T) 

p

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

seurat_obj <- MetacellsByGroups(
  seurat_obj = seurat_obj,
  group.by = c("cell_type", "Sample"), 
  reduction = 'harmony', 
  k = 25, 
  max_shared = 10, 
  ident.group = 'cell_type'
)

# normalize metacell
seurat_obj <- NormalizeMetacells(seurat_obj)

seurat_obj <- SetDatExpr(
  seurat_obj,
  group_name = "INH",
  group.by='cell_type',
  assay = 'RNA', 
  slot = 'data'
)

seurat_obj <- TestSoftPowers(
  seurat_obj,
  networkType = 'signed'
  )

plot_list <- PlotSoftPowers(seurat_obj)

wrap_plots(plot_list, ncol=2)

power_table <- GetPowerTable(seurat_obj)
head(power_table)

seurat_obj <- ConstructNetwork(
  seurat_obj,
  tom_name = 'INH'
)

PlotDendrogram(seurat_obj, main='INH hdWGCNA Dendrogram')

seurat_obj <- ScaleData(seurat_obj, features=VariableFeatures(seurat_obj))

seurat_obj <- ModuleEigengenes(
 seurat_obj,
 group.by.vars="Sample"
 )

# harmonized module eigengenes:
hMEs <- GetMEs(seurat_obj)

# module eigengenes:
MEs <- GetMEs(seurat_obj, harmonized=FALSE)

seurat_obj <- ModuleConnectivity(
  seurat_obj,
  group.by = 'cell_type', group_name = 'INH'
  )

seurat_obj <- ResetModuleNames(
  seurat_obj,
  new_name = "INH-M"
)

p <- PlotKMEs(seurat_obj, ncol=5)

p

modules <- GetModules(seurat_obj) %>% subset(module != 'grey')

head(modules[,1:6])

hub_df <- GetHubGenes(seurat_obj, n_hubs = 10)

head(hub_df)

saveRDS(seurat_obj, file='hdWGCNA_object.rds')