BioF3 组学数据分析

04 WGCNA 与共表达模块

导出日期：2026年6月27日

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 模块在整合分析中有两个用途：

  1. 特征降维：用 ME 代替上千个基因，减少后续整合模型的输入维度。比如把 20 个模块的 ME 和蛋白组数据一起送进 MOFA2。
  2. 跨层验证：在 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 上跑

下一步

接着深入：

横向延伸：

参考资源