WGCNA 加权共表达网络入门

跟完这篇，您能从一份表达矩阵出发，走通 WGCNA 的标准链条：摆正数据、定软阈值 β、把上万个基因划成若干「共表达模块」、把模块和您关心的性状关联起来、再从关键模块里挑出 hub 基因；同时您会认清样本量和预处理这两处最容易翻车的地方。

先想清楚 WGCNA 在给您什么

WGCNA（加权基因共表达网络分析）是一种无监督方法。它不看您的分组，而是问一个更基础的问题：在所有样本里，哪些基因的表达一起高、一起低，总是同进同退？把这些同进同退的基因聚成一团，就是一个模块（module）。

它最终给您三样东西：

模块：一批批协同表达的基因，各用一种颜色命名（grey 是「没归进任何模块」的散兵）。
模块特征基因（module eigengene，ME）：每个模块的「代表值」，可以理解为该模块所有基因表达的第一主成分，用一列数就概括整个模块在各样本上的高低。
模块—性状关联：用 ME 去和临床性状（分期、分型、生存、某个表型打分）算相关，找出「和这个性状最相关的模块」，再从中挑 hub 基因。

一句话：WGCNA 把「上万个基因」压成「十几个模块」，再让模块去对接您的生物学问题。

先摆正数据：方向、单位、样本量

这一步偷懒，后面全错。三件事必须先对：

1. 方向要转过来。 平时的表达矩阵是「行基因、列样本」，而 WGCNA 要的是行样本、列基因，正好是转置。这是新手最常见的错，务必先转再算。

library(WGCNA)
options(stringsAsFactors = FALSE)

# expr 是常规的「行基因、列样本」矩阵，转置成「行样本、列基因」
datExpr <- as.data.frame(t(expr))
dim(datExpr)   # 应为 样本数 × 基因数

# 检查有没有缺失过多或方差为零的基因/样本
gsg <- goodSamplesGenes(datExpr, verbose = 3)
gsg$allOK      # 为 TRUE 才算干净；否则按它标记的把坏基因/坏样本删掉

2. 单位要「归一化且方差平稳」。 WGCNA 的整套计算建立在相关系数上，直接喂原始 counts 会被少数高表达基因主导，相关结构失真。请喂经过归一化、且做过方差稳定的值：

数据类型	推荐喂给 WGCNA 的形式
RNA-seq（counts）	DESeq2 的 `vst()` / `rlog()`，或 limma-voom 的 log2-CPM
芯片	已归一化的 log2 表达值
不要喂	原始 counts、TPM/FPKM 未取对数、行列没转置的矩阵

VST/rlog 怎么得到，见 DESeq2 入门；矩阵整理的通用步骤见表达矩阵怎么整理。

3. 样本量要够。 WGCNA 靠样本间的相关来估共表达，样本太少相关不稳、模块不可复现。WGCNA 作者的建议是不要在少于 15 个样本上跑，20 个以上更稳妥。这是被广泛接受的经验下限，不是绝对红线，但样本个位数时得到的模块基本没有解释力，别硬做。

过滤基因、剔除离群样本

过滤基因只做一件安全的事：删掉在绝大多数样本里几乎不表达的基因（低表达噪声大，对网络没贡献）。

有一件事千万别做：不要先按「组间差异表达」筛一批基因再拿去做 WGCNA。WGCNA 的无标度假设建立在基因整体的相关结构上，先按差异挑基因会破坏这个结构、让模块划分产生偏倚。想缩减计算量时，人们常保留方差最大的一批基因，这是常见做法，但请注意按方差截断也会一定程度上改变网络拓扑——按表达水平过滤比按方差过滤更稳妥，能算得动就尽量少截。

# 惯例做法：保留在足够多样本里达到一定表达量的基因（阈值按数据调，不是铁律）
keep    <- colMeans(datExpr > 0) > 0.5        # 例:一半以上样本有表达
datExpr <- datExpr[, keep]

剔除离群样本：对样本做层次聚类，明显「吊在外面」的样本多半是技术异常，去掉再往下走。

sampleTree <- hclust(dist(datExpr), method = "average")
plot(sampleTree)   # 目视找离群样本;定一条高度线切开

# 按高度切树,保留主群(cutHeight 按图上情况定)
clust   <- cutreeStatic(sampleTree, cutHeight = 15000, minSize = 10)
datExpr <- datExpr[clust == 1, ]

cutHeight 没有通用数值，必须看着您自己的聚类图定——它取决于数据的尺度，照抄别人的数字没有意义。

定软阈值 β：让网络接近无标度

WGCNA 不做「相关就连、不相关就断」的硬切，而是把相关系数取 β 次幂作为连接强度——这就是「加权」。β 越大，弱相关被压得越狠、强相关越突出。定 β 的准则是让网络尽量满足无标度拓扑（scale-free topology）：少数基因高度连接、多数基因连接很少，这是真实生物网络的普遍特征。

用 pickSoftThreshold 扫一遍候选幂次，看两条线：

powers <- c(1:10, seq(12, 20, 2))

# networkType 选 "signed" 会区分正负相关(如今较常推荐);"unsigned" 只看相关强弱
sft <- pickSoftThreshold(datExpr, powerVector = powers,
                         networkType = "signed", verbose = 5)

sft$powerEstimate     # 函数给的建议 β(可能为 NA,需自己看图定)
head(sft$fitIndices)  # 含各幂次的无标度拟合 R² 与平均连接度

判读准则（常用惯例，按数据调，不是铁律）：

无标度拟合 R²：选第一个使 R² 达到约 0.8～0.9 的 β。R² 越高越接近无标度，但一味求高会把 β 推得过大、网络过稀。
平均连接度（mean connectivity）：随 β 增大而下降，选点时它应已降到相对较低、趋于平稳的位置。
两者权衡：在「R² 够高」的前提下取尽量小的 β。

如果 R² 怎么都上不去（样本偏少时常见），用 WGCNA 作者给的经验表按样本数兜底：

样本数	signed / signed hybrid 网络	unsigned 网络
< 20	18	9
20–30	16	8
30–40	14	7
> 40	12	6

注意这是拟合失败时的兜底默认，能从拟合曲线清楚选出 β 时，优先用曲线。

建网络、划模块

定好 β，就可以建网络并把基因切成模块。最省事的是一步到位的 blockwiseModules，它把「算邻接矩阵 → 算拓扑重叠 TOM → 层次聚类 → 动态剪枝分模块 → 合并相似模块」一次做完。

net <- blockwiseModules(
  datExpr,
  power         = 6,          # 用上一步定的 β
  networkType   = "signed",
  TOMType       = "signed",
  minModuleSize = 30,         # 模块最小基因数,默认 30(太小的模块不稳)
  mergeCutHeight = 0.25,      # ME 相关高于 0.75 的模块合并,默认 0.25
  numericLabels = TRUE,
  saveTOMs      = TRUE,
  verbose       = 3
)

# 把数字标签转成颜色;0 号 = grey = 未归入任何模块
moduleColors <- labels2colors(net$colors)
table(moduleColors)          # 看每个模块多少基因

几个参数的含义值得记住：

参数	作用	常用默认
`minModuleSize`	一个模块至少多少基因，低于此不成模块	30
`mergeCutHeight`	两模块 ME 相关高到一定程度就合并，避免碎片化	0.25（即相关 > 0.75 合并）
`TOMType` / `networkType`	是否区分正负相关；一般与软阈值步保持一致	`"signed"`

想更透明地控制每一步（自定义 TOM、手动剪枝），可以用「分步版」：adjacency() → TOMsimilarity() → hclust() → cutreeDynamic() → mergeCloseModules()。结果与一步法一致，只是每一步都摊开可调。

画出「基因树 + 模块颜色带」可以直观看模块划分：

plotDendroAndColors(
  net$dendrograms[[1]],
  moduleColors[net$blockGenes[[1]]],
  "Module colors",
  dendroLabels = FALSE, addGuide = TRUE
)

把模块和性状关联起来

这一步才把无监督的模块接回您的生物学问题。做法是：先算每个模块的 ME，再让 ME 和性状矩阵求相关。

nSamples <- nrow(datExpr)

# 每个模块一列 ME;orderMEs 让相似模块相邻,看热图更顺
MEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs <- orderMEs(MEs)

# traits: 行样本、列性状的数值矩阵(分类性状先转成 0/1 等数值)
moduleTraitCor <- cor(MEs, traits, use = "p")
moduleTraitP   <- corPvalueStudent(moduleTraitCor, nSamples)

结果是一张「模块 × 性状」的表：每格一个相关系数和一个 p 值。用带标注的热图看最直观：

textMatrix <- paste0(signif(moduleTraitCor, 2), "\n(",
                     signif(moduleTraitP, 1), ")")

labeledHeatmap(
  Matrix = moduleTraitCor,
  xLabels = colnames(traits),
  yLabels = names(MEs),
  textMatrix = textMatrix,
  colors = blueWhiteRed(50)
)

怎么读这张图：

先看颜色和数值：格子越红（正相关）或越蓝（负相关）、绝对值越大，说明该模块整体表达和该性状越同向或反向。
再看 p 值：相关要显著（p < 0.05 是常用惯例）才值得往下追。性状和模块成百上千对，别只挑最好看的一格讲故事。
锁定目标模块：与您关心的性状相关最强、又显著的那个（或几个）模块，就是下一步挖 hub 基因的对象。

挑 hub 基因：MM 与 GS 两把尺

选定目标模块后，从里面挑 hub 基因靠两个指标，缺一不可：

模块成员度（Module Membership，MM，也叫 kME）：某基因表达和模块 ME的相关。MM 高 = 这个基因很「像」整个模块的代表，处在模块核心。
基因显著性（Gene Significance，GS）：某基因表达和性状的相关（绝对值）。GS 高 = 这个基因本身和性状关系紧。

# MM:每个基因对每个模块 ME 的相关
geneMM <- as.data.frame(cor(datExpr, MEs, use = "p"))

# GS:每个基因对某个性状(取 traits 里一列)的相关
geneGS <- as.data.frame(cor(datExpr, traits$stage, use = "p"))

hub 基因 = 在目标模块里 MM 高、同时 GS 也高的基因。 把目标模块内所有基因画成「MM（横轴）对 GS（纵轴）」的散点，落在右上角的就是重点候选——它们既是模块核心，又和性状强相关。这类基因在模块内连接度（intramodular connectivity）通常也最高。

# 模块内连接度:直接从表达算,也是衡量 hub 的常用指标
kIM <- intramodularConnectivity.fromExpr(datExpr, moduleColors, power = 6)

# 若只想快速拿「每个模块一个最核心基因」
hubs <- chooseTopHubInEachModule(datExpr, colorh = moduleColors,
                                 power = 6, type = "signed")

一句要紧的提醒：这里的 hub 是「统计意义上的网络枢纽」，不等于「生物学关键基因」。 MM/GS 高只说明它在共表达结构里居中、和性状相关，不代表因果或可成药。hub 基因是候选、是假设，真正的作用要靠后续实验或独立数据验证。挑出后常见的下一步是把目标模块的基因整体拿去做富集，看它落在哪些通路，见富集分析。

几个最容易翻车的坑

WGCNA 出错往往不在代码，而在数据和解读。逐条核对：

样本太少还硬跑。少于 15 个样本时相关不稳、模块换个子集就变，结论不可复现。样本不够，先别做 WGCNA。
批次/技术结构没清干净。WGCNA 对批次极其敏感——技术差异也会「同进同退」，被当成模块。有明显批次时，先用 ComBat 或 limma::removeBatchEffect 校正，再进网络。
先按差异表达筛基因。这会破坏无标度假设、让模块偏倚。过滤只按表达水平做，别按组间差异做。
矩阵没转置 / 喂了原始 counts。样本没放到行、或喂了未做方差稳定的值，是最高频的两个低级错，前面「先摆正数据」一节务必落实。
把模块颜色当固定含义。颜色只是标签，turquoise、blue 这些名字不带任何生物学意义，也不跨数据集通用，别用颜色去附会。
过度解读 hub。统计枢纽 ≠ 关键基因，见上一节。

在平台上跑一遍

不想配 WGCNA 环境、也不想跟软阈值和 TOM 计算缠斗，可以在百沐一下不写代码做这套流程：上传做过方差稳定的表达矩阵和一张性状表，用一句话说清目的，例如「按无标度拟合定软阈值、划模块，并把模块和分期、生存关联，给出目标模块的 hub 基因」，组学智能体会替您完成软阈值选择、模块划分、模块—性状热图和 hub 候选，并附可复现的 R 代码，您只需核对数据方向、样本量和阈值是否合理。用法见组学智能体使用指南。