WGCNA共表达网络分析全流程：从表达矩阵到基因模块挖掘

拿到转录组表达矩阵以后，很多人第一反应是跑差异表达，筛出一堆显著基因，然后再去补一个GO/KEGG富集分析。但真到解释生物学意义的时候，往往会卡住：差异基因太多太散，彼此之间到底怎么协同工作？哪些基因模块跟临床性状、表型数据最相关？哪些基因是模块里的“核心枢纽”？这些问题如果不借助共表达网络，很难回答得让人信服。WGCNA（Weighted Gene Co-expression Network Analysis，加权基因共表达网络分析）解决的就是这个环节。

这是一套在R语言里完成的系统性分析流程，从表达矩阵清洗、样本质检，到软阈值选择、网络构建、模块识别，再到模块和性状的关联分析，一步扣一步。作为“全步骤”系列的第一篇，我会把从零开始跑通WGCNA上游核心流程的完整思路写清楚，所有代码都附上逐行解读，不只是让你能复制粘贴，更要让你知道每一行在干什么、为什么这么干。如果你正拿着转录组表达数据不知道下一步怎么挖，或者刚接触WGCNA看着官方教程一头雾水，这篇应该能帮你省掉不少折腾的时间。

1. 先搞清楚WGCNA到底在解决什么问题

1.1 共表达网络的基本逻辑：基因不是孤立工作的

我们可以把一个样本里成千上万个基因的表达量想象成一群人的工作状态。某些基因的表达量在所有样本里总是同涨同跌，说明它们大概率在同一个生物学通路里协同工作，或者受同一个上游转录因子调控。WGCNA把这种“同步变化”的关系抽象成网络：每个基因是一个节点，两个基因之间有没有边、边的权重有多大，取决于它们的表达量在所有样本里的相关性强不强。

传统做法是直接计算两两基因间的Pearson相关系数，然后设一个阈值，比如相关系数大于0.8就认为有连接。这样做的最大问题是：阈值是拍脑袋定的，而且硬性切割会把很多弱但真实存在的共表达关系丢掉。WGCNA的不同之处在于，它不设硬阈值，而是用一个软阈值（soft-thresholding power）对相关系数做幂指数加权，让网络尽量符合无标度拓扑特征——简单说，就是让网络里少数基因拥有大量连接，大多数基因只有少量连接，这是很多真实生物网络具备的特征。

1.2 WGCNA和普通相关分析的核心差别

普通相关分析是在“基因对”层面做检验，几万个基因会产生上亿个基因对，多重检验校正会让结果变得非常保守，而且你很难从海量相关关系里看出整体结构。WGCNA的处理方式是先构建全基因网络，再用层次聚类的方法把网络划分成若干个模块（module），每个模块里的基因表达模式高度相似。这样分析单元就从“单个基因”升级成了“基因模块”，后续做模块与性状关联、筛选枢纽基因、富集分析，逻辑上都更清晰。

1.3 一篇WGCNA能回答哪些问题

在实际项目中，WGCNA最常见的用途有三个：

• 找到与某个性状（比如疾病/正常、药物响应、发育阶段、产量高低）显著相关的基因模块。
• 在每个模块内部筛选高连接度的枢纽基因（hub gene），作为后续实验验证的候选靶点。
• 把模块基因拿去做功能富集，解释模块代表的生物学通路。

这套流程不是万能的，它对样本量有一定要求，对输入数据质量也很敏感，但这些细节放在后面实操部分再讲。先把逻辑框架搭好，后面每一步代码才不会变成“黑盒操作”。

2. 数据准备：两个文件的格式和你最容易踩的坑

2.1 表达矩阵的格式：行是样本，列是基因

WGCNA的输入数据格式非常明确：表达矩阵必须是“行=样本，列=基因”的二维数据框。我们经常拿到的原始表达矩阵是“行=基因，列=样本”，所以读入后通常要做一次转置。这也是新手第一处容易出错的地方，很多报错都源自行列搞反了。

library(WGCNA) # 关闭字符串自动转因子，避免后面矩阵运算出幺蛾子 options(stringsAsFactors = FALSE) # 读入表达矩阵，行名是基因ID，列名是样本名 datExpr_raw <- read.csv("expression_data.csv", row.names = 1, check.names = FALSE) # 转置成 行=样本，列=基因 datExpr0 <- as.data.frame(t(datExpr_raw))

这里有个细节：check.names = FALSE是为了防止R把样本名里的横杠、空格之类改写成点号，否则后面和性状文件匹配样本名时经常会莫名其妙对不上。

2.2 表达量数据到底该用什么

很多人问，WGCNA能不能直接输入DESeq2得到的标准化counts？理论上可以，但我更推荐用FPKM、TPM这类已经校正过基因长度和测序深度的表达量，而且要注意数据的分布形态。标准的RNA-seq表达矩阵数值跨度很大，直接拿来做相关性计算会被高表达基因主导，一般都需要做log2(x + 1)变换，让数据更接近正态分布。

# 如果数据还没有log变换，做一个log2(x+1)变换 datExpr0 <- log2(datExpr0 + 1)

如果你是拿芯片数据或者qPCR数据来跑，也要先确认数值范围合理。表达量差距在几个数量级、又没有做变换的数据跑出来，模块结构通常会非常碎，而且很难解释。

2.3 基因过滤：挑表达量稳定且非零的基因

WGCNA官网教程里的示例数据是已经预处理好的，但我们自己处理的数据里经常有大量低表达基因。这类基因的“表达变化”很多是测序噪声，计算相关矩阵时会把网络搅乱，还会显著拖慢计算速度。

我一般会分两步过滤：

第一步，用goodSamplesGenes检查缺失值和标准差为零的基因，把它们剔除。

# 检查样本和基因的基本质量 gsg <- goodSamplesGenes(datExpr0, verbose = 3) # 如果存在不合格的基因或样本，直接过滤掉 if (!gsg$allOK) { datExpr0 <- datExpr0[gsg$goodSamples, gsg$goodGenes] }

第二步，过滤表达量太低的基因，以及所有样本里表达量基本不变的基因，并用中位绝对偏差（MAD）筛出表达量变化最明显的那些基因。

# 计算每个基因在所有样本中的中位数和MAD gene_median <- apply(datExpr0, 2, median) gene_mad <- apply(datExpr0, 2, mad) # 保留中位数表达量较高、且变异较大的基因 keep <- gene_median > 1 & gene_mad > 0.5 datExpr <- datExpr0[, keep]

是不是一定要做这一步？如果基因数量在2万左右，机器配置也足够，网络构建其实也能跑完。但过滤之后模块会更稳定，后续模块-性状关联的显著性也更容易出现。这里需要把握一个度：如果过滤太狠，可能会丢掉重要基因；如果不过滤，噪声又会影响结果。比较稳妥的做法是先做常规低表达过滤，再用MAD筛掉后20%的基因，保留约8000到15000个基因进行网络构建。

2.4 样本聚类：肉眼剔除异常样本

这一步很多人会跳过，但我觉得它比后面的参数调优更关键。样本聚类树可以直观地暴露问题：如果一个样本和其他样本离得非常远，说明它的整体表达模式异常，可能是实验批次差异、样品污染或者数据预处理出了问题。

# 样本聚类 sampleTree <- hclust(dist(datExpr), method = "average") # 画图观察 pdf("sample_cluster.pdf", width = 12, height = 6) plot(sampleTree, main = "Sample clustering to detect outliers", sub = "", xlab = "", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2) abline(h = 100, col = "red") dev.off()

如果图里出现明显的离群样本，直接用下标把它剔除：

# 假设样本聚类图中第3个样本明显离群 datExpr <- datExpr[-3, ]

cutHeight阈值不是固定不变的，要看聚类树的高度分布来定。我的习惯是先用目测选一个把大多数样本聚在一起、只把极少数离群样本切出去的高度，不用刻意追求统一标准。

2.5 性状数据文件怎么整理

性状数据是WGCNA做模块关联分析的核心输入。格式要求是“行=样本，列=性状”，样本名要和表达矩阵的行名完全一致。性状可以是连续变量，比如年龄、血压、药物浓度；也可以是分组变量，但需要转成数值型，比如疾病组=1、对照组=0，或者多分组设计使用0/1哑变量编码。

# 读取性状数据 datTraits <- read.csv("sample_traits.csv", row.names = 1, check.names = FALSE) # 查看两个数据集的样本名交集 common_samples <- intersect(rownames(datExpr), rownames(datTraits)) datExpr <- datExpr[common_samples, ] datTraits <- datTraits[common_samples, , drop = FALSE]

这里要特别注意，必须保证表达矩阵和性状文件的样本一一对应。我一开始跑的时候就因为表达矩阵和性状文件样本顺序不一致，直接栽过跟头。用intersect统一样本后，后续分析就不会出现张冠李戴的问题。

3. 软阈值选择：别只会默认选9，要学会看两张图

3.1 无标度拓扑准则到底是什么

WGCNA里最核心的一个概念是“软阈值”power。简单理解，它就是一个加权系数，把基因间的相关系数取绝对值的power次方，得到基因之间的邻接权重。power越大，弱相关被压制得越厉害，网络的稀疏程度也越高。

但power并不是越大越好，WGCNA选择power的依据是让网络尽可能地符合无标度拓扑：我们希望网络中存在少数连接度极高的“枢纽节点”，而大部分节点连接度较低。衡量网络是否符合这个特征的指标是SFT.R.sq，也就是无标度拓扑拟合指数。通常来说，我们希望这个值尽量高，尤其是要超过0.85。

3.2 pickSoftThreshold的代码和输出

选择软阈值不需要自己瞎试，WGCNA包提供了pickSoftThreshold函数：

# 选择一系列候选power值 powers <- c(1:10, seq(from = 12, to = 30, by = 2)) # 计算不同power下的无标度拟合指数和平均连接度 sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 5) # 把结果整理成一个数据框查看 fit <- sft$fitIndices print(fit[, c("Power", "SFT.R.sq", "mean.k.", "median.k.")])

输出结果会看到一张类似下面的表格：

眼睛不要只盯着哪个power大，要同时看SFT.R.sq和mean.k.。判断标准是：选择最小的、让SFT.R.sq首次进入0.85以上区间的power，同时平均连接度不能太低，否则网络会过度稀疏，模块识别失去意义。

3.3 组合图怎么看

通常还会画一张双面板的组合图，左边是power和拟合指数SFT.R.sq的折线，右边是power和平均连接度的趋势线：

pdf("soft_threshold.pdf", width = 10, height = 5) par(mfrow = c(1, 2)) cex1 <- 0.9 # 左图：SFT.R.sq plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], xlab = "Soft Threshold (power)", ylab = "Scale Free Topology Model Fit (R^2)", type = "n", main = "Scale independence") text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], labels = powers, col = "red", cex = cex1) abline(h = 0.85, col = "red", lty = 2) # 右图：平均连接度 plot(sft$fitIndices[, 1], sft$fitIndices[, 5], xlab = "Soft Threshold (power)", ylab = "Mean Connectivity", type = "n", main = "Mean connectivity") text(sft$fitIndices[, 1], sft$fitIndices[, 5], labels = powers, col = "red", cex = cex1) dev.off()

这里有一个WGCNA代码里比较反直觉的地方：左图纵轴为什么是-sign(x[,3]) * x[,2]？因为当power比较低时，无标度拟合指数可能是负斜率，乘上负号能让所有点都朝上显示，方便观察。看的时候只要看绝对值就可以了。

3.4 如果R平方一直上不了0.85怎么办

这是实际分析中最常见的问题之一。遇到这种情况，先不要急着怀疑数据有问题。有几种处理思路：

• 检查是否没有过滤低表达基因，噪声太大影响了网络结构，先回去重新做数据清洗。
• 尝试更广的power范围，比如seq(1, 40, by = 2)。
• 对表达矩阵再做一次更严格的标准差筛选。
• 样本量太小（比如少于15个样本）时，SFT.R.sq可能很难达到0.85，此时可以考虑退而求其次，选择曲线转折点附近的power，同时结合网络生物学可解释性来选定。

我个人在临床小样本数据上遇到过几次R平方只能到0.7的情况。我的做法是选取一个中间偏上的power，比如根据曲线趋势选10或者12，然后看后续模块是否稳定、模块与性状的关联是否有意义。如果模块结果乱七八糟，再回头调整power。

4. 核心代码块：网络构建与模块识别

4.1 blockwiseModules函数参数逐个拆解

选定power之后，就到了整个WGCNA分析的重头戏——构建网络并识别模块。标准代码并不复杂，核心就是blockwiseModules：

# 用选定的power构建网络 net <- blockwiseModules( datExpr, power = 9, TOMType = "unsigned", minModuleSize = 30, reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE, pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = "TOM-blockwise", verbose = 3 )

这个函数到底做了什么？它可以拆成三层理解：先计算基因间的相关性矩阵，再换算成邻接矩阵，然后计算TOM相似度（拓扑重叠），最后基于TOM相异度做层次聚类和动态剪枝，得到模块。下面逐个参数说明，每个参数我都会结合自己踩过的坑来讲。

power：就是软阈值。选9还是选12，取决于第3节pickSoftThreshold的结果，不要无脑用默认值。

TOMType：推荐用"unsigned"。它的意思是把负相关也当作共表达关系，因为基因调控网络里存在大量负调控，一个转录因子抑制下游基因时，它们的表达模式是负相关的。如果改用"signed"，则只保留正相关基因之间的共表达关系，模块会更保守。

minModuleSize：最小模块基因数。默认值是30，但是基因数少的平台数据可以降到20甚至10。设置太小会导致模块数量爆炸，很多只有一两个基因的“迷你模块”没有生物学意义。

mergeCutHeight：相似模块合并的阈值。默认0.25，意思是模块间特征基因相关性高于0.75就需要考虑合并。如果后续发现模块多而碎，可以适当调低mergeCutHeight到0.2，或者在后期手动做一次模块合并。

numericLabels：如果为TRUE，模块会以数字编号命名，输出的是0、1、2…… 如果为FALSE，会用颜色命名，比如turquoise、blue、brown。我个人建议代码阶段先用数字，导出绘图时再映射成颜色，这样在后续匹配基因时更方便。

pamRespectsDendro：是否在动态剪枝后继续用PAM（Partitioning Around Medoids）对模块边界做进一步细化。如果设为TRUE，生成的模块更紧凑，但有时会过度切割；如果设为FALSE，会更容易得到比较大的、稳健的模块。默认是FALSE，我一般也保持FALSE。

saveTOMs和saveTOMFileBase：是否把TOM矩阵保存到本地。TOM矩阵是WGCNA里最大的中间产物，会占很大磁盘空间，但保存下来可以避免后续重复计算，尤其是在做不同power对比时能节省大量时间。

4.2 从net对象里能拿到什么

blockwiseModules跑完之后，net对象里包含了做后续分析所需的几乎所有关键结果：

# 查看模块划分结果 table(net$colors) # 查看模块特征基因 net$MEs # 查看模块数量 length(unique(net$colors))

net$colors是一个和基因列表长度相同的向量，记录了每个基因被分配到了哪个模块，0表示未进入任何模块的灰色基因。net$MEs是每个样本在每个模块上的特征基因表达值，后面做模块-性状关联时主要靠它。

4.3 把聚类树画出来看看整体结构

模块识别之后，最直观的可视化方式是画聚类树和模块颜色条：

pdf("dendrogram_modules.pdf", width = 12, height = 8) plotDendroAndColors( net$dendrograms[[1]], moduleColors = net$colors, groupLabels = "Module colors", main = "Gene dendrogram and module colors", dendroLabels = FALSE, addGuide = TRUE, hang = 0.03 ) dev.off()

这张图是整个WGCNA结果里我最先看的一张。理想情况下，树状图会分成几个明显的大分支，每个分支对应一个颜色模块。如果看到颜色条像斑马线一样反复横跳，或者模块数量特别多，就要考虑调整minModuleSize、mergeCutHeight或者回去检查数据预处理。

刚接触WGCNA的时候，我拿到模块结果第一反应是赶紧去看模块和性状的关联，结果画出来的聚类树乱得没法看。后来才意识到，聚类树的稳定程度本身就能反映数据质量，如果这一层就有问题，后面所有统计都会跟着出问题。

4.4 基因数量与模块大小的合理范围

一个合理的网络通常会有10到30个模块，其中最大的模块基因数量占10%到30%，同时会有几个中等模块和少量小模块。如果所有基因都堆在同一个巨大模块里，说明power选得太低，模块没有分开；如果模块数量超过50个，要么过滤不充分，要么minModuleSize设得太小。这些标准不是绝对死线，但能帮助快速判断结果是否靠谱。

5. 模块合并与模块-性状关联分析

5.1 为什么要做模块合并

blockwiseModules里的mergeCutHeight参数会自动完成相似模块的合并。但不同版本代码或者手动调整后，可能还需要自己再检查一遍。相似模块指的是两个模块的特征基因（ME，Module Eigengene）高度相关，它们本质上是同一群基因的不同变体，合并后更容易解释，也减少后续多重检验的次数。

如果需要手动合并，代码是这样的：

# 计算模块特征基因 MEs <- moduleEigengenes(datExpr, net$colors)$eigengenes # 对模块进行合并 merge <- mergeCloseModules(datExpr, net$colors, cutHeight = 0.25, verbose = 3) # 合并后的模块颜色 mergedColors <- merge$colors # 合并后的模块特征基因 mergedMEs <- merge$newMEs

合并后一定要重新把模块颜色画一遍确认，不要直接跳过可视化。

5.2 模块-性状关联的计算方式

每个模块的ME是这个模块所有基因表达模式的第一主成分，代表整个模块在样本间的主要变化趋势。把ME和每个性状做相关分析，就能得到模块与性状的关联矩阵。

# 计算模块特征基因与性状的相关性 nSamples <- nrow(datExpr) moduleTraitCor <- cor(mergedMEs, datTraits, use = "p") moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples)

corPvalueStudent是WGCNA为这个场景专门封装好的函数，它会基于相关系数和样本量给出p值，比手算方便很多。

模块与性状关联的结果通常画成热图，横轴是性状，纵轴是模块，每个格子里的数字是相关系数和括号里的p值：

pdf("module_trait_heatmap.pdf", width = 8, height = 10) textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = "") dim(textMatrix) <- dim(moduleTraitCor) labeledHeatmap( Matrix = moduleTraitCor, xLabels = colnames(datTraits), yLabels = colnames(mergedMEs), ySymbols = colnames(mergedMEs), colorLabels = FALSE, colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.6, main = "Module-trait relationships" ) dev.off()

重点看两个指标：相关系数的绝对值大小，以及p值的显著性水平。比如某个模块与“疾病状态”的相关系数是0.72，p值小于0.001，那这个模块就是你后续要重点挖掘的候选模块。

5.3 从模块里挑基因：基因显著性与模块成员度

找到一个与性状显著相关的模块后，还需要知道模块里哪些基因起主导作用。WGCNA提供了两个经典指标：

• 基因显著性（GS, Gene Significance）：基因表达量与性状之间的相关绝对值，表示该基因与性状的关联强度。
• 模块成员度（MM, Module Membership）：基因与该模块特征基因的相关性，表示该基因在模块内的核心程度。

# 选择一个性状，比如第1列 trait_column <- 1 # 计算每个基因与性状的相关性 geneTraitSignificance <- as.data.frame(cor(datExpr, datTraits[, trait_column], use = "p")) GS <- as.numeric(geneTraitSignificance$V1) # 计算每个基因与模块特征基因的相关性 MM <- as.data.frame(cor(datExpr, mergedMEs, use = "p"))

然后可以拿GS和MM做散点图，筛选兼具模块核心地位且和性状显著相关的基因。比如模块内MM > 0.8且GS > 0.2的基因，往往就是值得后续验证的候选枢纽基因。

5.4 灰色模块怎么处理

grey模块是WGCNA里一个特殊的存在，它是所有没能被划分到任何模块的基因集合。灰色模块通常不会和性状显著相关，如果有也一样要关注一下，可能意味着数据里还存在另一种独立的表达模式，值得单独做一次亚聚类分析。多数情况下，灰色模块基因不参与后续分析。

6. 实操中的零散问题与排查建议

6.1 样本量太小怎么办

WGCNA对样本量的最低要求保守说至少要有15到20个样本，如果少于这个数，基因相关性估计会非常不稳定，模块结果也很难重复。真遇到小样本数据，我有几个折中经验：把minModuleSize调低到10左右，power选择稍微偏大一点，让网络更稀疏，同时使用signed TOM而不是unsigned，往往会稳一些。但也要明确，小样本条件下的WGCNA结果只能作为探索性分析，不适合直接下强结论。

6.2 计算太慢、内存爆掉

基因数量接近两万时，blockwiseModules默认会做分块计算，因为一次性计算全部基因的TOM矩阵对内存压力很大。如果仍然卡死，建议先确认R是64位版本，然后适当降低保留基因数，比如从15000降到10000。还可以开启多线程：

# 启用多线程，靠CPU核心数决定 enableWGCNAThreads()

注意，这条命令要在加载WGCNA后、构建网络之前运行，而且Windows系统下多线程支持不如Linux/macOS稳定。

6.3 结果可重复性问题

WGCNA里有一些步骤依赖随机性，尤其是样本量不特别大的时候，模块识别结果可能会有轻微波动。想保证结果可重复，建议在脚本开头设置随机种子：

set.seed(2024)

另外，分析过程中所有中间结果都要及时保存：

save(datExpr, datTraits, net, mergedColors, mergedMEs, file = "wgcna_step1.rda")

这样即使后面改参数，也不需要从最开始的读取文件重新跑一遍。

6.4 包版本差异

WGCNA包这些年更新不算频繁，但不同小版本的默认参数可能略有差异，比如blockwiseModules里的某些参数在新版本提示deprecated。如果碰到函数调用报错，先看包自带的NEWS文档，再看官方教程。网上很多老教程用的代码在新版本下可能会跑不通，这不是你写错了，需要根据报错信息微调函数名或参数。

最后再分享一个偷懒技巧

整套流程里最耗时的往往是不同power下的网络构建对比。我会在跑正式分析前先拿一小部分基因，比如随机抽3000个基因做一个快速测试，看看模块数量是否合理、聚类树是否稳定。参数基本满意后，再用全部基因跑正式版本。这个小技巧能省下大量反复调试的时间。

另外，blockwiseModules生成的TOM文件占空间非常大，分析做完如果没有特殊需要，记得及时清理，不然一个项目下来几百GB一点也不夸张。把模块结果、基因颜色、特征基因这些核心结果保留好就足够了。

这一篇把WGCNA从数据准备到模块-性状关联的完整前半程梳理完了，代码基本可以直接照着改路径和数据跑通。下一篇我会继续写模块内部的可视化细节，包括基因网络导出到Cytoscape、hub gene筛选、以及怎么把模块基因批量提交给富集分析工具。跑代码的过程中遇到具体的报错和诡异结果，欢迎照着这篇的排查思路先自己试一圈，多数问题都出在数据格式和样本匹配上。