scRNA-hdWGCNA共表达网络分析【1】:（一文掌握）详细讲解教程更新及网络图可视化函数

library(Seurat)
setwd("/home/data_analysis")
GSE282955_ARH_Sex_by_Nutr <- readRDS("~/data_analysis/GSE282955_ARH_Sex_by_Nutr.rds")
dim(GSE282955_ARH_Sex_by_Nutr)
#[1] 35318 93282

#抽样
set.seed(123)
# 从 meta.data 中抽样
cells_keep <- GSE282955_ARH_Sex_by_Nutr@meta.data %>%
  rownames_to_column("cell") %>%
  group_by(cell_type) %>%  # 替换为你注释的细胞类型列名
  sample_frac(size = 0.1) %>%
  ungroup() %>%
  pull("cell")
# 创建子集对象
GSE282955_ARH_Sex_by_Nutr_sub <- subset(GSE282955_ARH_Sex_by_Nutr, cells = cells_keep)
DimPlot(GSE282955_ARH_Sex_by_Nutr_sub,label = T)+NoLegend()
dim(GSE282955_ARH_Sex_by_Nutr_sub)
#[1] 35318  9330
save(GSE282955_ARH_Sex_by_Nutr_sub,file='GSE282955_ARH_Sex_by_Nutr_sub.RData')

#devtools::install_github('smorabit/hdWGCNA', ref='dev')
library(hdWGCNA)
package.version('hdWGCNA')

library(hdWGCNA)
library(WGCNA)
library(Seurat) 
library(tidyverse) 
library(igraph)
library(cowplot)
library(patchwork)
library(dplyr)
library(ggplot2)
library(stringr)

setwd('/Users/ks_ts/Documents/公众号文章/hdWGCNA复现/')
load("./GSE282955_ARH_Sex_by_Nutr_sub.RData")

#将原始数据赋予一个新的对象，以保留原来数据
seurat_obj <- GSE282955_ARH_Sex_by_Nutr_sub

#setup 
seurat_obj <- SetupForWGCNA(
  seurat_obj,#seurat对象
  gene_select = "fraction", #这里使用所有基因
  fraction = 0.05, #基因表达百分比，0.05表示这个基因必在5%的细胞中表达
  wgcna_name = "ARH" #此次WGCNA分析项目名称
)

# construct metacells  in each group
seurat_obj <- MetacellsByGroups(
  seurat_obj = seurat_obj,
  group.by = c("cell_type3", "Batch"), 
  reduction = 'pca', # select the dimensionality reduction to perform KNN on
  k = 25, 
  ident.group = 'cell_type3' 
)

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

##set up the expression matrix
seurat_obj <- SetDatExpr(
  seurat_obj,
  group_name = "Agrp", #需要进行分析的celltype
  group.by='cell_type3', #metadata中包含celltype的列
  assay = 'SCT', 
  layer = 'data' 
)

#Select soft-power threshold
# Test different soft powers:
seurat_obj <- TestSoftPowers(
  seurat_obj,
  networkType = 'signed' #对于生物学数据，网络type更适合signed，明确区分正负调控
)

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

wrap_plots(plot_list, ncol=2)

# construct co-expression network:
seurat_obj <- ConstructNetwork(
  seurat_obj, soft_power=15,
  setDatExpr=FALSE,
  tom_name = 'Agrp' 
)

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

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="Batch"
)

hMEs <- GetMEs(seurat_obj)

# compute eigengene-based connectivity (kME):
seurat_obj<- ModuleConnectivity(
  seurat_obj,
  group.by = 'cell_type3', 
  group_name = 'Agrp'
)

# rename the modules
seurat_obj <- ResetModuleNames(
  seurat_obj,
  new_name = "Agrp_M"#原始的module名称是阿拉伯数字1，2，3. new_name相当于在他们前面添加了前缀Agrp_M1，Agrp_M2.。。
)

# get the module assignment table:
modules <- GetModules(seurat_obj) %>% subset(module != 'grey')

#查看目前的module
modules <- GetModules(seurat_obj)

#获取module level对应的module color，也就是颜色顺序与levels(modules$module)一致
mod_colors <- dplyr::select(modules, c(module, color)) %>%
  distinct %>% arrange(module) %>% .$color
mod_colors

#这里一共有14个module需要重新分配颜色
my_color = c("#D52126", "#88CCEE","#FEE52C","#117733","#CC61B0","#99C945","#2F8AC4",
          "#332288","#E68316", "#661101", "#F97B72", "#DDCC77", "#11A579", "#89288F")
seurat_obj <- ResetModuleColors(seurat_obj, my_color)

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

PlotKMEs(seurat_obj, ncol=4)

# compute gene scoring for the top 25 hub genes by kME for each module
# with UCell method
library(UCell)
seurat_obj <- ModuleExprScore(
  seurat_obj,
  n_genes = 25,
  method='UCell'
)

plot_list <- ModuleFeaturePlot(
  seurat_obj,
  features='scores', # plot the hub gene scores
  order='shuffle',
  reduction = 'tsne',# order so cells are shuffled
  ucell = TRUE,# depending on Seurat vs UCell for gene scoring
  point_size =0.1)

wrap_plots(plot_list, ncol=4)+ plot_annotation(title = "Hub genes signature scores")

#plot所有module的网络，并将文件导出储存在Agrp_all_ModuleNetworks文件夹
ModuleNetworkPlot(
  seurat_obj,
  outdir = './Agrp_all_ModuleNetworks'
)

ModuleNetworkPlot(
  seurat_obj, 
  mods = 'Agrp_M4',
  n_inner = 10, # number of genes in inner ring
  n_outer = 15, # number of genes in outer ring
  n_conns = Inf, # show all of the connections
  vertex.label.cex=1 # font size
)

HubGeneNetworkPlot(
  seurat_obj,
  n_hubs = 10, #每个module需要展示的hub基因数量
  n_other=5,#每个module需要展示的非hub基因数量
  edge_prop = 0.75,#对网络中的边数进行降采样
  mods = c('Agrp_M1','Agrp_M2','Agrp_M4','Agrp_M6','Agrp_M7','Agrp_M8','Agrp_M9') #为了演示选择几个module,可以选择任意module
)

##Applying UMAP to co-expression networks
seurat_obj <- RunModuleUMAP(
  seurat_obj,
  n_hubs = 10, 
  n_neighbors=15, 
  min_dist=0.3 # min distance between points in UMAP space
)

ModuleUMAPPlot(
  seurat_obj,
  edge.alpha=0.5,
  sample_edges=TRUE,
  edge_prop=0.1, # proportion of edges to sample (20% here)
  label_hubs=3 ,# how many hub genes to plot per module?
  keep_grey_edges=FALSE
)

ks_hdWGCNA_UMAPnet(seurat_obj=seurat_obj,
                   wgcna_name = 'ARH',
                   edge_prop=0.05,
                   top_hub=3,
                   r_text=10)