04 WGCNA 与共表达模块
WGCNA（Weighted Gene Co-expression Network Analysis）通过构建基因共表达网络，把上万个基因压缩成几十个模块，每个模块代表一组协同变化的基因。在多组学整合中，WGCNA 的模块可以作为"特征压缩"的手段，也可以用来检验某个模块在不同数据集或不同组学层之间是否保守。
基本流程
library(WGCNA)
allowWGCNAThreads()
# 输入：行是样本，列是基因（注意转置）
datExpr <- t(rna_final)
# 1. 选择软阈值（soft threshold）
powers <- c(1:20)
sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 3)
# 画 scale-free topology fit 图
plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
xlab = "Soft Threshold (power)",
ylab = "Scale Free Topology Model Fit (signed R^2)",
type = "n", main = "Scale independence")
text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2],
labels = powers, col = "red")
abline(h = 0.85, col = "red")
选择 R^2 首次超过 0.85 的 power 值。通常 RNA-seq 数据在 6-12 之间。
构建网络与模块识别
# 2. 一步法构建网络
net <- blockwiseModules(datExpr,
power = 8,
TOMType = "unsigned",
minModuleSize = 30,
reassignThreshold = 0,
mergeCutHeight = 0.25,
numericLabels = TRUE,
pamRespectsDendro = FALSE,
verbose = 3)
# 模块颜色
moduleColors <- labels2colors(net$colors)
table(moduleColors)
# 画聚类树 + 模块颜色
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
"Module colors", dendroLabels = FALSE,
hang = 0.03, addGuide = TRUE, guideHang = 0.05)
模块与表型的关联
每个模块用 module eigengene（ME，即模块内基因表达的第一主成分）来代表。然后计算 ME 与临床表型的相关：
MEs <- net$MEs
# 计算 ME 与表型的相关
trait_data <- data.frame(
subtype = as.numeric(factor(col_data$subtype)),
stage = as.numeric(factor(col_data$stage))
)
cor_ME_trait <- cor(MEs, trait_data, use = "p")
pval_ME_trait <- corPvalueStudent(cor_ME_trait, nrow(datExpr))
# 热图展示
library(ComplexHeatmap)
Heatmap(cor_ME_trait, name = "Correlation",
cell_fun = function(j, i, x, y, w, h, fill) {
if (pval_ME_trait[i, j] < 0.05)
grid.text("*", x, y)
})
模块保守性分析
如果你有两个独立数据集（比如 TCGA 和自己的队列），可以检验某个模块在第二个数据集中是否保守。这用 modulePreservation 函数：
# 准备第二个数据集
datExpr2 <- t(rna_validation)
# 设置多集数据
multiExpr <- list(
Discovery = list(data = datExpr),
Validation = list(data = datExpr2)
)
multiColor <- list(Discovery = moduleColors)
# 计算保守性
mp <- modulePreservation(multiExpr, multiColor,
referenceNetworks = 1,
nPermutations = 200,
randomSeed = 1,
verbose = 3)
# Zsummary > 10 表示高度保守，2-10 中等，< 2 不保守
stats <- mp$preservation$Z$ref.Discovery$inColumnsAlsoPresentIn.Validation
print(stats[, c("moduleSize", "Zsummary.pres")])
在多组学整合中的角色
WGCNA 模块在整合分析中有两个用途：
- 特征降维：用 ME 代替上千个基因，减少后续整合模型的输入维度。比如把 20 个模块的 ME 和蛋白组数据一起送进 MOFA2。
- 跨层验证：在 RNA 层发现的模块，检查对应基因在蛋白层或甲基化层是否也有一致的模式。
常见坑
坑 1：soft threshold 选得太低或太高
R² 没达到 0.85 就随便选了个 power = 4，网络不满足 scale-free 假设，模块划分没有生物学意义。反之选 power = 20 会让网络太稀疏，大多数基因被归入灰色模块。一定要看 pickSoftThreshold 的拟合曲线，选 R² 首次超过 0.85 的最小 power。
坑 2：模块稳定性差却不做验证
WGCNA 的模块划分对参数（minModuleSize、mergeCutHeight）敏感。换个 seed 或调一下 mergeCutHeight 模块就完全不同，说明结果不稳定。应该用 modulePreservation 在独立数据集上验证，或者用 bootstrap 检验模块的可重复性。Zsummary < 2 的模块不应该报告。
坑 3：样本量不足导致相关矩阵不可靠
WGCNA 依赖基因间的相关矩阵。样本 < 20 时，相关系数的估计噪声很大，构建出的网络是假信号。官方建议至少 15 个样本，实际经验 30 以上效果才稳定。样本太少时考虑用其他降维方法（PCA、MOFA）替代。
坑 4：module eigengene 与表型的相关当成因果
ME 和 subtype 显著相关不代表"这个模块驱动了亚型形成"。ME 可能只是 confounding（比如该模块反映的是增殖速率，而增殖本身就和分期高度相关）。需要结合富集分析和文献验证，给出合理的生物学解释。
不想本地装环境？在 BioF3 上跑
下一步
接着深入：
- 05 MOFA2 因子分析实战 — 把 WGCNA 模块的 eigengene 作为 MOFA 输入做多组学整合
- 03 跨层相关性探索 — 模块间的跨层相关检验
横向延伸：
- bulk RNA-seq 06 高级建模 — 转录组层面的共表达和批次校正
- 蛋白质组 07 STRING 网络 — 蛋白互作网络的构建和模块分析
参考资源
离线资料下载
手册 HTML / PDF 已在后台预生成，点击后直接下载网站静态资源。