2026版单细胞测序scRNA-seq入门实战（十三）：hdWGCNA共表达网络分析——解码基因的“社交网络”

本节概览：

WGCNA/hdWGCNA简介：理解基因共表达网络分析的核心思想及其在单细胞数据中的价值。
为什么单细胞需要hdWGCNA：解析单细胞数据稀疏性、异质性背景下传统WGCNA的局限与hdWGCNA的解决方案（元细胞Metacell）。
数据准备：从已注释的Seurat对象出发，筛选适合构建网络的基因。
核心概念：软阈值、模块（Module）、特征基因（Module Eigengene）的深度解读。
全流程代码实战：从SetupForWGCNA到模块识别、Hub基因鉴定与可视化。
模块-性状关联：将基因模块与细胞类型、样本分组等“性状”建立关联。

在前几节中，我们完成了细胞类型注释、差异表达分析、功能富集和细胞通讯分析。这些分析帮助我们回答了“细胞是什么”、“它们表达了什么差异基因”、“参与了什么通路”等问题。

但有一个更深层的问题被忽略了：基因之间是如何协同工作的？

在生物学中，基因很少“单打独斗”——它们往往以共表达模块的形式协同发挥作用，共同执行某个生物学功能。WGCNA（Weighted Gene Co‑expression Network Analysis）正是通过分析基因之间的共表达模式，识别出那些“步调一致”的基因模块。然而，传统的WGCNA是为Bulk RNA-seq设计的，直接应用于单细胞数据会面临数据稀疏性和细胞异质性两大挑战。

hdWGCNA（high-dimensional WGCNA）正是为解决这些问题而生的。它通过元细胞（Metacell） 策略聚合相似的单细胞，有效缓解稀疏性问题，同时支持在特定细胞类型或状态中构建网络，捕捉细胞类型特异性的共表达模式。

一、WGCNA的核心概念：一张图读懂

在开始写代码之前，先理清几个WGCNA中的核心术语：

术语	英文	含义
共表达网络	Co-expression Network	以基因为节点、以表达相关性为边的网络
软阈值	Soft Threshold	将相关性矩阵转换为加权邻接矩阵的幂指数β，用于放大强相关、压制弱相关
模块	Module	一组高度共表达的基因，通常用颜色命名（如“蓝色模块”）
特征基因	Module Eigengene (ME)	模块中所有基因表达的第一主成分，代表整个模块的“浓缩”表达特征
Hub基因	Hub Gene	模块内部连接度最高的基因，通常被认为是模块的核心调控者
元细胞	Metacell	将多个相似的单细胞合并为一个“超级细胞”，用于缓解单细胞数据的稀疏性

💡 一句话理解：WGCNA就像给基因做“社交网络分析”——把表达模式相似的基因拉到一个群里（模块），再找出群里的“意见领袖”（Hub基因）。

二、准备工作：安装hdWGCNA

在R中安装hdWGCNA：

devtools::install_github('smorabit/hdWGCNA', ref='dev')

三、全流程代码实战（可直接运行）

阶段一：环境准备与数据载入

# 清空环境，设置基础选项
rm(list = ls())
options(stringsAsFactors = FALSE)

# 加载所有必需的R包
library(Seurat)
library(tidyverse)
library(cowplot)
library(patchwork)
library(WGCNA)          # 传统WGCNA核心算法
library(hdWGCNA)        # 单细胞版本WGCNA
library(UCell)          # 用于后续模块评分的快速计算

# 设置图形主题和随机种子
theme_set(theme_cowplot())
set.seed(12345)

# 启用多线程加速（根据你的CPU核心数调整，这里使用12线程）
enableWGCNAThreads(nThreads = 12)

# 创建输出文件夹
main_dir <- getwd()
if (!dir.exists("hdWGCNA")) dir.create("hdWGCNA")
setwd("hdWGCNA")

# 读取保存的含有细胞类型注释的Seurat对象
seurat_obj <- readRDS("../4_scRNA_celltype.rds"))
# 查看数据中的细胞类型，确认分组信息
print(unique(seurat_obj$cell_type))

阶段二：基因筛选与元细胞构建

# ============================================================
# 2. 设置 WGCNA 实验（筛选用于构建网络的基因）
# ============================================================
seurat_obj <- SetupForWGCNA(
  seurat_obj,
  gene_select   = "fraction",      # 基于表达比例筛选基因
  fraction      = 0.05,            # 至少在 5% 的细胞中表达才纳入
  wgcna_name    = "P5"             # 本次实验的标签（会存入seurat_obj@misc$P5）
)
# 查看入选的基因数量（通常为几千个）
cat("入选基因数：", length(seurat_obj@misc$P5$wgcna_genes), "\n")

# ============================================================
# 3. 构建元细胞（Metacell）解决数据稀疏性问题
# ============================================================
seurat_obj <- MetacellsByGroups(
  seurat_obj,
  group.by     = "cell_type",     # 按细胞类型分别构建元细胞
  reduction    = "umap",          # 使用UMAP空间坐标进行邻近聚类
  k            = 25,              # 每个元细胞聚合25个单细胞（15-30之间均可）
  max_shared   = 10,              # 防止过度共享，避免元细胞之间重叠过多
  ident.group  = "cell_type"      # 元细胞的身份标记
)
# 对元细胞的表达矩阵进行标准化（方法同常规Seurat）
seurat_obj <- NormalizeMetacells(seurat_obj)

阶段三：软阈值选择与网络构建

软阈值β是WGCNA最关键的参数，它决定了网络对相关性强弱的敏感度。

# ============================================================
# 4. 设置表达矩阵（指定用哪部分数据来计算共表达）
# ============================================================
# 这里我们使用全部细胞类型，设定 group.by = "cell_type" 
# 使算法在每种细胞类型内部分别评估基因的相关性
seurat_obj <- SetDatExpr(
  seurat_obj,
  group_name   = unique(seurat_obj$cell_type),
  group.by     = "cell_type",
  assay        = "RNA",
  slot         = "data"           # 使用 log-normalized 表达值
)

# ============================================================
# 5. 自动评估并选择软阈值（Soft-thresholding power）
# ============================================================
# 测试一系列软阈值，评估其是否能使网络接近“无标度拓扑”
seurat_obj <- TestSoftPowers(
  seurat_obj,
  networkType = "signed"          # 有符号网络（保留正负相关信息）
)
# 可视化结果，寻找 R² 首次达到 0.8 以上对应的软阈值
pdf("WGCNA_soft_power.pdf", width = 12, height = 8)
wrap_plots(PlotSoftPowers(seurat_obj), ncol = 2)
dev.off()

# 查看详细的阈值评估表格
power_table <- GetPowerTable(seurat_obj)
print(head(power_table))

# 手动指定软阈值（此处以9为例，请根据上一步图表中的R²>=0.8来选择）
# 如果不指定，ConstructNetwork 会根据评估结果自动选择一个
soft_power <- 9

# ============================================================
# 6. 构建共表达网络并识别基因模块
# ============================================================
seurat_obj <- ConstructNetwork(
  seurat_obj,
  soft_power   = soft_power,       # 使用刚才选定的软阈值
  tom_name     = "P5_TOM"          # 拓扑重叠矩阵（TOM）的存储名称
)

# 绘制模块聚类树状图（树状图的不同颜色分支代表不同模块）
pdf("PlotDendrogram.pdf", width = 12, height = 6)
PlotDendrogram(seurat_obj, main = "hdWGCNA Dendrogram")
dev.off()

# ============================================================
# 7. 导出模块成员信息
# ============================================================
# 提取每个基因所属的模块（注意存储路径在 seurat_obj@misc$P5 下）
module_genes <- seurat_obj@misc$P5$wgcna_modules
# 查看各模块的基因数量
table(module_genes$module)
# 保存到本地
write.csv(module_genes, "module_genes.csv", row.names = FALSE)
saveRDS(module_genes, "module_genes.rds")

阶段四：模块特征基因与Hub基因鉴定

这是生物学解读的黄金步骤——找出每个模块的“代言人”（ME）和“意见领袖”（Hub基因）。

# ============================================================
# 8. 计算模块特征基因（MEs）与连通性（kME）
# ============================================================
# 提示：ModuleEigengenes 需要数据已经过 ScaleData
# 对筛选出的高变基因进行缩放（加速计算）
# seurat_obj <- ScaleData(seurat_obj, features = VariableFeatures(seurat_obj))

# 计算模块特征基因（即模块内所有基因的第一主成分）
seurat_obj <- ModuleEigengenes(
  seurat_obj,
  group.by.vars = NULL            # 若存在批次效应可设为 "orig.ident"
)

# 获取模块特征基因矩阵（行=细胞，列=模块）
MEs <- GetMEs(seurat_obj, harmonized = FALSE)

# # 计算每个基因与其所属模块特征基因的相关性（kME值）
# # kME值越高，说明该基因在模块内越“核心”
# seurat_obj <- ModuleConnectivity(
#   seurat_obj,
#   group.by   = "cell_type",       # 按细胞类型分组计算
#   harmonized = FALSE
# )
# 
# # 重命名模块，使其更清晰（例如 M1, M2, M3...）
# seurat_obj <- ResetModuleNames(
#   seurat_obj,
#   new_name = "M"
# )
# 
# # 绘制各模块的 kME 分布图（展示每个模块中基因的重要性排名）
# pdf("PlotKMEs.pdf", width = 15, height = 2.5)
# PlotKMEs(seurat_obj, ncol = 5)
# dev.off()
# 1. 获取模块分配表（包含每个基因属于哪个模块）
modules <- GetModules(seurat_obj)
modules <- modules[modules$module != "grey", ]  # 剔除未分配的基因

# 2. 获取元细胞的标准化表达矩阵（基因 × 元细胞）
expr_matrix <- GetAssayData(seurat_obj, assay = "RNA", layer = "data")

# 3. 获取模块特征基因矩阵（元细胞 × 模块）
MEs <- GetMEs(seurat_obj, harmonized = FALSE)  # 注意：行名是元细胞，列名是模块

# 4. 手动计算每个基因与其模块 ME 的 Pearson 相关系数（kME）
library(WGCNA)   # 如果不加载，可以用 base::cor

kME_list <- list()
for (mod in unique(modules$module)) {
  cat("正在处理模块:", mod, "\n")
  
  # 该模块的基因列表
  genes_in_mod <- modules$gene_name[modules$module == mod]
  if (length(genes_in_mod) < 3) next
  
  # 提取该模块基因的表达数据（基因 × 元细胞）
  expr_mod <- expr_matrix[genes_in_mod, , drop = FALSE]
  
  # 提取该模块的特征基因值（一个向量，长度 = 元细胞数）
  me_vec <- MEs[, mod]
  
  # 对每个基因计算与 ME 的相关性
  cor_values <- apply(expr_mod, 1, function(gene_expr) {
    cor(gene_expr, me_vec, use = "pairwise.complete.obs")
  })
  
  # 整理为数据框
  temp_df <- data.frame(
    gene_name = names(cor_values),
    module = mod,
    kME = as.numeric(cor_values),
    stringsAsFactors = FALSE
  )
  kME_list[[mod]] <- temp_df
}

# 合并所有模块
kME_df <- do.call(rbind, kME_list)

# 查看结果（应该没有 NaN）
summary(kME_df$kME)  # 确认有数值

# 保存（可选）
write.csv(kME_df, "kME_manual.csv", row.names = FALSE)

# 5. 提取每个模块的 Top 20 Hub 基因
library(dplyr)

hub_df <- kME_df %>%
  group_by(module) %>%
  arrange(desc(kME)) %>%
  slice_head(n = 20) %>%
  ungroup()

# 查看
head(hub_df)

# 6. 绘制条形图（替代 PlotKMEs）
library(ggplot2)

p <- ggplot(hub_df, aes(x = reorder(gene_name, kME), y = kME, fill = module)) +
  geom_bar(stat = "identity", width = 0.7) +
  facet_wrap(~module, scales = "free_x", ncol = 5) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, size = 6),
    strip.background = element_rect(fill = "lightgrey"),
    legend.position = "none"
  ) +
  labs(x = "", y = "kME (Correlation with Module Eigengene)",
       title = "Top 20 Hub Genes per Module")

# 保存
ggsave("HubGenes_Top20_Manual.pdf", plot = p, width = 16, height = 10)
ggsave("HubGenes_Top20_Manual.png", plot = p, width = 16, height = 10, dpi = 300)


# 也可以只画 Top 10（如果图太密）
hub_top10 <- kME_df %>%
  group_by(module) %>%
  arrange(desc(kME)) %>%
  slice_head(n = 10) %>%
  ungroup()

p10 <- ggplot(hub_top10, aes(x = reorder(gene_name, kME), y = kME, fill = module)) +
  geom_bar(stat = "identity", width = 0.7) +
  facet_wrap(~module, scales = "free_x", ncol = 5) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, size = 7),
    strip.background = element_rect(fill = "lightgrey"),
    legend.position = "none"
  ) +
  labs(x = "", y = "kME", title = "Top 10 Hub Genes per Module")

ggsave("HubGenes_Top10_Manual.pdf", plot = p10, width = 14, height = 8)


# ============================================================
# 9. 提取分析结果并保存
# ============================================================
# 获取完整的模块分配表（包含基因、模块、kME等）
modules <- GetModules(seurat_obj) %>% subset(module != "grey")  # 剔除未分配的基因
write.csv(modules, "hdWGCNA_module_full.csv", row.names = FALSE)

# 提取每个模块的 Top 20 Hub 基因（最重要的核心基因）
hub_df <- GetHubGenes(seurat_obj, n_hubs = 20)
write.csv(hub_df, "hub_genes_top20.csv", row.names = FALSE)

# 保存完整的Seurat对象（包含所有WGCNA计算结果），方便后续直接画图
saveRDS(seurat_obj, "hdWGCNA_object.rds")
# 同时保存当前R工作空间镜像
save.image("hdWGCNA_workspace.RData")

阶段五：发表级可视化

hdWGCNA提供了丰富的可视化函数，帮你把结果直接变成论文可用的图片。

# ============================================================
# 10. 结果可视化
# ============================================================

# 10.1 模块特征基因在UMAP上的空间分布
# 首先基于Hub基因计算每个细胞的模块评分（使用UCell快速算法）
seurat_obj <- ModuleExprScore(
  seurat_obj,
  n_genes = 20,           # 使用每个模块前20个Hub基因计算评分
  method  = "UCell"
)

# 绘制 MEs 的 FeaturePlot（展示模块整体的激活区域）
plot_list_me <- ModuleFeaturePlot(
  seurat_obj,
  features = "MEs",
  order    = TRUE         # 将高表达细胞置于顶层，避免遮挡
)
pdf("ModuleFeaturePlot_MEs.pdf", width = 15, height = 10)
wrap_plots(plot_list_me, ncol = 3)
dev.off()

# 10.2 雷达图：展示不同细胞类型中模块的富集程度
# 准备分组变量（通常用细胞类型或聚类群）
pdf("ModuleRadarPlot.pdf", width = 9, height = 6)
ModuleRadarPlot(
  seurat_obj,
  group.by = "cell_type",
  axis.label.size = 3,
  grid.label.size = 3
)
dev.off()

# 10.3 模块间相关性热图（识别相似功能的模块）
pdf("ModuleCorrelogram.pdf", width = 6, height = 6)
ModuleCorrelogram(seurat_obj)
dev.off()

# 10.4 模块在不同细胞类型中的平均表达（DotPlot）
# 将模块特征基因添加到元数据中，方便使用Seurat原生绘图
MEs_all <- GetMEs(seurat_obj, harmonized = FALSE)
seurat_obj@meta.data <- cbind(seurat_obj@meta.data, MEs_all)

# 获取所有非灰色模块的名称
mods <- levels(GetModules(seurat_obj)$module)
mods <- mods[mods != "grey"]

p <- DotPlot(
  seurat_obj,
  features = mods,
  group.by = "cell_type"
) +
  RotatedAxis() +
  scale_color_gradient2(low = "blue", mid = "grey95", high = "red")
ggsave("Dotplot_Modules_across_CellTypes.pdf", p, width = 8, height = 6)

四、结果解读与生物学意义

4.1 看模块与性状的关联

雷达图（ModuleRadarPlot）和热图能直观告诉你哪些模块在特定细胞类型中特异性高表达。比如某个模块在“小胶质细胞”中富集，说明该模块可能与小胶质细胞的功能密切相关。

4.2 看Hub基因

导出 hub_genes_top20.csv 文件，聚焦每个模块排名前5的基因。这些基因往往是模块功能的核心执行者。如果一个模块富集了“炎症反应”通路，那么它的Hub基因可能就是炎症的关键调控因子。

4.3 看模块的生物学功能

提取感兴趣模块的基因列表，结合我们之前学的 GO/KEGG富集分析（第八篇推文），对模块基因进行功能注释。例如：

五、常见问题与避坑指南

5.1 细胞数量太少怎么办？

hdWGCNA要求每个分组至少有50-100个细胞。如果某类细胞数量不足：

合并相似的细胞类型（如OPC + Oligodendrocyte）。
降低k值（元细胞大小），从25降到15，获得更多元细胞。
如果实在太少，建议跳过WGCNA，改用其他方法（如SCENIC）。

5.2 如何确定软阈值β？

运行TestSoftPowers后，查看生成的WGCNA_soft_power.pdf。

重点关注 Scale-free R² 曲线，选择 R² 首次达到 0.8 以上 对应的β值。
如果所有β值的R²都很低（<0.8），说明数据可能不适合WGCNA，可能是细胞异质性太高或基因筛选太宽松。

5.3 模块太多或太少怎么办？

模块太多（>50个）：提高软阈值β，或增大ConstructNetwork中的minModuleSize参数（默认30）。
模块太少（<5个）：降低软阈值β，或降低minModuleSize。

5.4 运行报错“无法分配足够内存”？

单细胞WGCNA确实比较吃内存。对于几万个细胞，建议在服务器上运行（64G+内存）。

六、总结

至此，我们完成了从Seurat对象到hdWGCNA共表达网络分析的完整流程，代码可直接复现。通过这一套流程，你能够：

将数万个基因降维成几十个有生物学意义的“模块”；
识别每个模块的核心驱动基因（Hub基因）；
建立模块与细胞类型/样本分组/临床性状的关联；
将模块与功能富集、细胞通讯、分化轨迹串联，讲出完整的生物学故事。

WGCNA帮助你从“单个基因差异”的视角上升到“基因模块协同”的层面——不是问“哪个基因变了”，而是问“哪些基因一起变了、为什么变、由谁驱动”。这种全局视角尤其适合解析复杂的生物学过程，如发育、疾病进展和细胞状态转换。

这里是两栖生物手册，中科院生物医学博士，持续记录生物医学实验和生物信息学笔记，干湿结合两不误～正努力成为最贴心负责的生信数据分析者和实验技术分享者，欢迎大家关注～

2026版单细胞测序scRNA-seq入门实战（十三）：hdWGCNA共表达网络分析——解码基因的“社交网络”