scRNA-hdWGCNA共表达网络：模块特征差异与相关性精选分析

1 概述

此前两篇教程已完整梳理了hdWGCNA的核心分析流程：从共表达网络构建与基因模块识别，到富集分析及大规模球状网络图的生成。网络框架搭建完毕，模块也已划定，但研究者最关心的往往是——究竟哪些模块在特定条件下呈现出显著变化？这些变化与实验设计的生物学变量之间存在怎样的关联？

本文直接回应这两个核心问题。我们将分别执行差异模块分析（DME）和模块-性状相关性分析，将模块的生物学解读落到实处，使分析结果具备可解释性。

首先加载所需数据与依赖包：

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复现/')
seurat_obj <- readRDS('./Agrp_hdWGCNA_obj.rds')

2 差异模块分析：Differential module eigengene (DME) analysis

前期分析已揭示数据中存在的基因模块及其大致功能。接下来需要回答的关键问题是：这些模块在不同条件下的表达模式是否存在实质性差异？以本演示数据集为例，实验设计包含性别（Sex）和营养状态（Nutr_State）两个变量。例如，雌性禁食组相较于雄性禁食组，是否存在某些模块出现特异性的上调或下调？这正是DME分析要解决的核心问题。

hdWGCNA提供FindDMEs函数用于组间差异模块检测，其语法与Seurat中的FindMarkers高度一致，学习成本极低。

group1 <- seurat_obj@meta.data %>% 
  subset(sexXnutr == 'F_Fast' & cell_type3 == 'Agrp') %>% rownames
group2 <- seurat_obj@meta.data %>% 
  subset(sexXnutr == 'M_Fast' & cell_type3 == 'Agrp') %>% rownames

# 结果为 group1 vs group2
DMEs <- FindDMEs(
  seurat_obj,
  barcodes1 = group1,
  barcodes2 = group2,
  harmonized = TRUE,  # 默认T，前面对ME做了harmony矫正
  test.use = 'wilcox',
  wgcna_name = 'ARH'
)

head(DMEs)

##               p_val  a vg_log2FC pct.1 pct.2   p_val_adj   module
## Agrp-M10 0.0003634394  0.98248499 0.855 0.676 0.005088152 Agrp-M10
## Agrp-M11 0.0004068578 -0.86945801 0.809 0.919 0.005696009 Agrp-M11
## Agrp-M6  0.0006566490 -0.35091115 0.800 0.635 0.009193085  Agrp-M6
## Agrp-M13 0.0016294910  0.70616560 0.818 0.608 0.022812875 Agrp-M13
## Agrp-M4  0.0020561082  0.01432549 0.955 0.865 0.028785515  Agrp-M4
## Agrp-M12 0.0031978977  1.30940087 0.891 0.784 0.044770567 Agrp-M12

获取结果后，使用PlotDMEsLollipop进行可视化展示。但这里有一个关键细节：该函数要求DME结果中的模块名称必须与seurat_obj中存储的模块名称完全一致，否则图表渲染会失败。实际运行中发现，FindDMEs执行后，模块名称中的下划线被自动替换为短横线——这是FindMarkers的经典问题。

验证seurat_obj中的模块名称：

mod_color_df <- GetModules(seurat_obj)
unique(mod_color_df$module)

## [1] grey     Agrp_M1  Agrp_M2  Agrp_M3  Agrp_M4  Agrp_M5  Agrp_M6  Agrp_M7 
## [9] Agrp_M8  Agrp_M9  Agrp_M10 Agrp_M11 Agrp_M12 Agrp_M13 Agrp_M14
## 15 Levels: grey Agrp_M1 Agrp_M2 Agrp_M3 Agrp_M4 Agrp_M5 Agrp_M6 ... Agrp_M14

对比DME结果中的模块名称：

unique(DMEs$module)

## [1] "Agrp-M10" "Agrp-M11" "Agrp-M6"  "Agrp-M13" "Agrp-M4"  "Agrp-M12"
## [7] "Agrp-M5"  "Agrp-M3"  "Agrp-M9"  "Agrp-M1"  "Agrp-M2"  "Agrp-M8" 
## [13] "Agrp-M7"  "Agrp-M14"

确实，下划线已全部变为短横线。修复方法很简单：手动将短横线替换回下划线即可。

DMEs$module <- gsub("-", "_", DMEs$module)
rownames(DMEs) <- DMEs$module

名称修正后即可顺利绘图。两种常规可视化方式：棒棒糖图适合快速目视模块的差异方向与显著性水平。

PlotDMEsLollipop(seurat_obj, DMEs, wgcna_name = 'ARH', 
                 pvalue = "p_val_adj") +
  theme_classic() +
  labs(title = "F_Fast vs M_Fast") +
  theme(plot.title = element_text(hjust = 0.5))

火山图则能更密集地展示信息，直观呈现哪些模块显著上调和下调：

PlotDMEsVolcano(seurat_obj, DMEs, wgcna_name = 'ARH') +
  labs(title = "F_Fast vs M_Fast") +
  theme(plot.title = element_text(hjust = 0.5))

3 模块-性状相关性分析：Module Trait Correlation

DME分析解答的是“哪些模块在组间存在差异”，而模块-性状相关性分析则将模块与具体的生物学变量或技术变量直接关联。这一逻辑与经典WGCNA完全一致：通过计算模块特征基因（ME）与性状的相关性，可以快速锁定那些与研究变量（如性别、营养状态、测序深度等）高度相关的模块。

操作上，使用ModuleTraitCorrelation函数。需注意该函数会基于指定的细胞类群进行计算——因为某些性状可能仅在特定细胞类型的特定模块中呈现相关性，而在其他细胞类型中不显著。因此先提取Agrp细胞子集，再进行分析。

# 提取Agrp细胞子集
agrp_obj <- subset(seurat_obj, cell_type3 == "Agrp")

# 确保因子顺序正确
agrp_obj$Sex <- as.factor(agrp_obj$Sex)
agrp_obj$Sex <- factor(agrp_obj$Sex, levels = c("F", "M"))
agrp_obj$Nutr_State <- as.factor(agrp_obj$Nutr_State)
agrp_obj$Nutr_State <- factor(agrp_obj$Nutr_State, levels = c("Fast", "Fed"))

# 指定要分析的性状
cur_traits <- c('Sex', 'Nutr_State', 'nCount_RNA', 'nFeature_RNA')

agrp_obj <- ModuleTraitCorrelation(
  agrp_obj,
  traits = cur_traits,
  group.by = 'cell_type3',
  wgcna_name = 'ARH'
)

分析结果存储于对象内，可通过GetModuleTraitCorrelation提取。返回的列表包含三个数据框：相关系数、p值及FDR校正后的p值。

Agrp_module_cor <- GetModuleTraitCorrelation(agrp_obj)
Agrp_module_cor

可视化采用热图形式，直观展示每个模块与各性状的相关性强度及显著性（以星号标注FDR）：

PlotModuleTraitCorrelation(
  agrp_obj,
  label = 'fdr',
  label_symbol = 'stars',
  text_size = 3,
  text_digits = 3,
  text_color = 'black',
  high_color = '#F97B72',
  mid_color = 'white',
  low_color = '#1798E5',
  plot_max = 0.2,
  combine = TRUE
)

该热图可快速识别出在性别或营养状态上具有显著偏好性的模块，为后续深入解析提供清晰而直接的切入点。

来源：互联网