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）两个变量。那么，有没有某些模块在雌性禁食组 vs 雄性禁食组中表现出了特异性上调或下调？这正是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
)
从这个热图可以快速锁定哪些模块在性别或营养状态上有显著偏好，为后续的深入挖掘提供了非常清晰的切入点。
