企业官网建设流程全解析

1. 认识WGCNA模块与性状展示图

1.1 这张图到底在展示什么

WGCNA 的“模块与性状展示图”，在转录组分析里可以算是出场率最高的可视化之一。做植物育种的用它看产量相关模块，做肿瘤研究用它找生存相关基因群，做微生物组分析的也常借鉴这套思路。这张图之所以通用，是因为它把“基因-基因”的共表达关系和“基因群-性状”的表型关联压在同一张图里，一眼就能锁定需要优先关注的模块。

从实际操作来看，一张标准的模块-性状展示图分好几个信息层。最上方或左侧是模块颜色标签，每种颜色代表一组表达模式高度相似的基因；中间部分是模块与性状的相关性矩阵，格子里通常写相关系数和 p 值，格子背景颜色从蓝色到红色渐变，表示负相关到正相关。如果你关心的是某个体性状，比如疾病分期、药物响应或作物株高，那么你只需要横向扫过对应列，找到颜色最深而且 p 值低于阈值的小格，再横向看对应的模块颜色，那就是候选模块。

那模块是怎么来的呢？WGCNA 首先计算基因之间表达谱的相关性，再通过层次聚类，把相关性高的基因聚成一个模块。每个模块被赋予一个颜色名称，比如 turquoise、blue、brown。模块与性状之间没有直接的计算接口，需要借助“模块特征基因”来搭桥。模块特征基因（module eigengene，ME）是对模块内全部基因表达矩阵做主成分分析后的第一主成分，它相当于整个模块表达趋势的“虚拟代表基因”。有了每个模块的代表表达谱，再跟每个性状做相关分析，于是就有了这张热图。

这个思路在实际项目中价值很大。比如在某个作物的抗逆项目里，基因数有两万多，逐个去找响应逆境的基因非常低效。先通过 WGCNA 把两万多个基因收敛成二十多个模块，再从中锁定与抗逆性状显著相关的两三个模块，范围一下子就缩小到几百个基因。之后再结合模块内的基因功能注释、蛋白互作网络，筛选出候选基因做验证。整套流程里，模块与性状展示图是最关键的“过滤闸门”。

1.2 图上每一处信息该怎么读

很多人第一次看图时，会纠结格子里先写相关系数还是先写 p 值。最常见的默认格式是相关系数在上、p 值在下，但这个可以自己定义。我更建议在论文方法学部分写清楚格式定义，免得审稿人误会。格子背景颜色表示相关性的方向和强度：红色调代表正相关、蓝色调代表负相关，颜色越深绝对值越大。默认的 blue-white-red 渐变里，白色区域表示相关性趋近于零。

需要特别提醒的是，不要只盯着 p 值，也要看相关系数的绝对值。有时候某个模块与性状相关的 p 值小于 0.001，但相关系数只有 0.3，说明样本量较大、相关关系稳定，但效应量其实一般。反过来，某个模块相关系数达到 0.7，但 p 值不显著，这通常是样本量太少造成的。做候选模块筛选时，我更倾向于综合看两个指标，优先选择“相关性强且 p 值显著”的一批，而不是只挑 p 值最小的那个。

还有一点，这张图是模块级别的分析结果，不是基因级别。也就是说，它告诉你“哪个基因群”值得关注，至于具体是群里哪个基因起主导作用，还需要看模块内成员的连接度、kME 值等指标。很多初学者把模块与性状图当成了结论图，以为画完就万事大吉，其实它只是筛选漏斗的开始。

2. 画图之前必须先确认的三件事

2.1 表达矩阵和性状表的格式规范

先说数据布局。表达矩阵应该行是基因、列是样本。性状矩阵应该行是样本、列是性状，样本顺序和数量要与表达矩阵的列完全一一对应。WGCNA 的很多函数默认这个布局，数据放反了轻则报错，重则画出一张完全不可信的网络图。

表达量数据不建议用原始 Counts 直接输入，因为测序深度差异会产生虚假的相关。芯片数据一般取对数，转录组数据我习惯先做 TPM 或 CPM 的对数转换，也可以先用 DESeq2 的 vst 函数做方差稳定变换。处理完之后，还要检查一下有没有缺失值。WGCNA 包里有一个好用的检查函数goodSamplesGenes()，它会同时检查基因和样本的缺失情况以及零方差基因，跑出来有 warning 就要先处理。

样本量方面，我的经验是少于 15 个样本的 WGCNA 需要非常谨慎。不是说绝对不能用，而是模块划分结果很容易不稳定。以前我接过一个只有 8 个样本的数据，跑出来的模块图看起来还挺整齐，但后续候选基因验证基本全军覆没。后来复盘才发现，模块的划分过度依赖了样本之间的噪声相关性。如果你样本少，建议补样本或者换一种思路。

性状的编码也需要动脑筋。连续型性状比如体重、株高、激素浓度，直接用原始数值就行。二分类性状比如患病/健康，我建议编码为 0/1，这样正相关和负相关的含义很直观。多分类无序性状最好转成哑变量，比如三个分组，可以生成两列 0/1。有些小伙伴喜欢把分组标成 1、2、3 当成连续变量塞进去，这也不是绝对不行，但前提是你得确认 1 到 3 的线性变化有生物学意义。如果分组本身没有梯度，这种编码很容易产生误导性结果。

2.2 软阈值选择与方法背后的逻辑

软阈值是新手最容易困惑的环节。WGCNA 构建的不是普通的相关性阈值网络，而是加权网络，基因间相关性会被取 β 次幂。选软阈值时，最常见的参考是看“无标度拓扑拟合指数图”和“平均连接度图”。默认准则是，选择使 R² 首次达到 0.8 以上、同时平均连接度不至于太小的最小 β 值。

为什么要在意无标度拓扑？生物学网络通常呈现“少数枢纽基因、大量外围基因”的分布，这和无标度网络的特征一致。如果 β 选得太小，网络过于稠密，模块会被拆成很多小碎块；β 太大，网络过于稀疏，大部分基因变成孤点，灰色模块占比升高。两种极端都不是想要的。

我自己的习惯是，把 1 到 20 的 β 值对应的模块数量、灰色基因占比整理成一个小表格，结合标准图来做决定。有一个常见情况是，R² 一直达不到 0.8，这种真实数据集里很常见，可能只有 0.7 左右。这时不用死磕，选一个候选 β，看最终模块划分是否合理再作微调。判断标准可以看模块数量：如果模块数量在十几到三十几个之间，基因大多能分配到有颜色的模块，灰色模块占比不高，那这个网络基本是可用的。

这个步骤要养成记录习惯，因为不同数据、不同标准化方式，甚至不同批次样本，选出来的 β 都会变。它不是固定的 6 或者 8，而是受数据分布影响很大的参数。做项目时把每个候选 β 的评估结果记下来，比只盯着两张标准图更容易做后续的比较和复现。

3. 模块与性状关联的核心计算细节

3.1 模块特征基因是怎么提炼出来的

模块特征基因的完整计算，是对模块内所有基因的表达矩阵做主成分分析，返回第一主成分。这里有一个质量指标值得留意：第一主成分对模块方差的解释率。解释率越高，说明模块中的基因表达越同步，模块信息越凝聚。如果解释率过低，比如低于 0.3，那这模块可能是一堆表达谱相似度有限的基因被勉强聚到一起，后续结果的可信度就得打问号。

在 R 中，moduleEigengenes()函数计算完 ME 之后，还会得到每个基因与该模块 ME 的相关性，官方称为 kME 或基因显著性。我们在做模块内候选基因筛选时，会把 kME 值较高的基因视为模块的核心基因。这个指标用途很广，后续看 hub 基因、做富集分析、挑实验验证对象，都会用到它。

为什么要用第一主成分而不用模块内所有基因的平均表达值？因为后者容易被少数高表达基因拉偏。第一主成分会根据每个基因表达方差自动确定权重，相当于“加权平均”，比“均权平均”稳健得多。转录组数据本身高维度、高噪声，加权操作通常能带来更稳定的模块代表值。理解了这一层，你就能明白为什么模块-性状图上的相关性，往往比单个基因与性状的相关性更可靠。

3.2 显著性的多重检验问题

模块与性状相关分析中，模块数量通常在十几个到几十个之间，这确实存在多重检验的问题。WGCNA 的标准函数默认不做事后校正，你会看到很多 p < 0.05 的格子。我的做法是两步走。第一步用 p < 0.05 做初筛，第二步对筛选出来的显著模块再判断是否符合生物学逻辑。如果某个模块与性状显著相关，但它的基因功能注释与目标性状完全不沾边，那这种相关很可能只是统计噪声。

也有朋友问，要不要直接把 FDR 校正后的 q 值画上去。我的看法是：如果模块数量少于 10 个，直接展示原始 p 值问题不大；如果模块数量明显多于 20 个，最好标记 BH 校正后的 q 值，或者在方法学部分做说明。发表论文时如果审稿人提出严格统计控制要求，再补一个 FDR 也不复杂。

另外还要注意，Pearson 相关对异常值非常敏感。如果某个样本在某个性状上明显异常，比如记录错误或者极端值，它会单独拉动相关系数。所以我一般会在正式出图前，先跑一遍模块-性状相关矩阵，瞄一眼有没有奇怪的极端相关。发现问题先核对原始记录，再去决定是否剔除样本，而不是直接把问题丢给统计算法处理。

4. 实操：从数据到成品图的完整流程

4.1 一份可以直接套用的代码

以下是一段我在实际项目中经常使用的核心流程，整理成了最小可运行框架。你可以根据自己的数据进行替换。

library(WGCNA) # 允许 WGCNA 线程设置，避免后续计算警告 allowWGCNAThreads() # 读取表达矩阵：行基因、列样本，然后转置为行样本、列基因 datExpr <- read.csv("expr.csv", header = TRUE, row.names = 1) datExpr <- as.data.frame(t(datExpr)) # 检查样本和基因的缺失情况 gsg <- goodSamplesGenes(datExpr, verbose = 3) if (!gsg$allOK) { datExpr <- datExpr[, gsg$goodGenes] } # 读取性状数据，行名与表达矩阵样本名一致 datTraits <- read.csv("trait.csv", header = TRUE, row.names = 1) # 挑选软阈值 powers <- c(1:20) sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 5) # 以选定的 power 构建网络和模块 net <- blockwiseModules( datExpr, power = 6, TOMType = "unsigned", minModuleSize = 30, reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE ) # 计算模块特征基因 MEs <- moduleEigengenes(datExpr, colors = net$colors)$eigengenes MEs <- orderMEs(MEs) # 模块与性状的相关性和 p 值 nSamples <- nrow(datExpr) moduleTraitCor <- cor(MEs, datTraits, use = "p") moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples) # 生成图上显示的文字 textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = "") dim(textMatrix) <- dim(moduleTraitCor) # 输出 PDF 格式的模块与性状展示图 pdf("module_trait.pdf", width = 8, height = 10) labeledHeatmap( Matrix = moduleTraitCor, xLabels = colnames(datTraits), yLabels = colnames(MEs), ySymbols = colnames(MEs), colorLabels = FALSE, colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.5 ) dev.off()

代码里值得关注的几个地方。TOMType = "unsigned"表示不区分正负相关，只要表达趋势一致就算连接，这个设置适合大多数转录组研究。如果你明确只关注正相关共表达，可以改成"signed"，模块切割结果会有差异。mergeCutHeight = 0.25是模块合并的阈值，两个模块的特征基因相关性高于这个阈值时会被合并，这能有效避免产生过多相似模块。minModuleSize = 30设的是模块最小基因数，样本量小或者想看到更细粒度模块时可以适当降低。

4.2 出图参数和图形观感调整

实际画图时，我一般会先输出一个 PDF，而不是直接在屏幕上看，因为热图的信息密度太高，屏幕显示容易失真。width和height要根据模块数和性状数调整。模块多、性状多时，图幅太小会导致格子挤在一起，文字重叠。

cex.text控制格子里文字的大小，模块数多的时候建议调低到 0.4 到 0.5，模块数少时可以调到 0.7 以上。blueWhiteRed(50)里的 50 表示颜色等级数，等级越多颜色过渡越细腻。如果你希望高相关区域更醒目，可以自定义颜色向量，把蓝色和红色设置得更深。

标签也是一个容易忽视的坑。默认的模块标签会是“MEturquoise”这种很长的名字，在热图的行标签里占地方。我通常会先做一轮模块重命名，把最关注的模块改成简洁的名字，比如“ME_wt”或“ME_groupA”，列名太长时整体图会显得乱。此外，setStdMargins = FALSE配合手动调整页面留白，可以避免标签被切掉。

5. 常见问题与排查手记

5.1 模块几乎全变成灰色怎么办

灰色是 WGCNA 里“未分配基因”的集合，如果灰色模块占基因总数的 30% 以上，说明网络模块化程度不高。常见原因有几个：软阈值选得太高导致网络太稀疏，minModuleSize设置过大，或者数据本身噪声偏高。

排查思路我一般固定走三步。第一步，降低软阈值往回调两三个档位，重新看模块数量和灰色占比；第二步，把minModuleSize从 30 降到 20 甚至 15，给基因更多分配机会；第三步，检查表达矩阵的标准化是否合理，原始 Counts 直接跑最容易出现大片灰色。还有一种容易被忽略的情况，是个别离群样本把整体相关性拉低了，这时候要先看样本聚类树，找出离群样本，分析原因后再决定是否剔除。

5.2 所有模块与性状都不显著怎么办

这是很多第一次跑 WGCNA 的项目会撞上的问题。一种可能是模块划分不合理，另一种可能是你对“显著”的预期太高。基因表达与复杂性状之间的相关，本来就不会像临床化验指标那样高，相关系数 0.5 以上已经算相当强了。我见过不少真实数据，生物学相关就落在 0.3 到 0.5，p 值能到 0.05 以下已经不错。

如果全不显著，先看样本量，再看性状本身的变异程度。如果某个性状在你所有样本里几乎不变，比如所有样本的株高都差不多，那它跟任何模块都不会有显著相关。还有一个排查点是性状编码，分类变量胡乱编码成连续数字，也可能把信号掩盖掉。把性状的分布画出来看一眼，往往能发现问题。

5.3 热图文字重叠和布局混乱

labeledHeatmap的默认布局有时不能满足复杂数据的需求。模块多时行标签太长，容易跟图右侧的内容冲突；性状名太长时会压迫格子区域。解决的思路有两个：一是把 PDF 图幅调大，增加留白；二是缩短标签文本，模块名和性状名都改为简短别名，图下加注释说明全称。

还有一个小技巧，格子里的文字内容可以调整。如果 p 值位数太多显得拥挤，可以只保留相关系数和显著性星号。虽然我建议论文里放完整数值，但在探索阶段，用星号标注可以更清爽。等确定好要展示的重点再补全细节。

5.4 模块数量异常多或异常少的判断

做过几轮 WGCNA 之后，你会对模块数量有个直觉：一般十来个个到三十几个比较常见。如果一次跑出 50 多个模块，说明切割过细，后续做富集分析会发现每个模块基因太少、功能不集中。反之只有三五个模块，基因群过大，筛选候选基因的精度不够。

模块数量异常时，优先调整两个参数：minModuleSize和mergeCutHeight。minModuleSize决定模块最小基因数，调大可以减少小碎块；mergeCutHeight决定模块间合并阈值，调低会让模块更难合并，模块数增加。两个参数可以组合起来试，每次改一个，记录结果，直到模块数量落在合理区间。

6. 让这张图真正服务于课题结论

6.1 从显著模块到候选基因的常规路线

模块与性状展示图完成了“模块筛选”这一步之后，接下来的工作才是重头戏。拿到显著模块后，第一件事就是把模块内的基因提取出来，做一个 GO 和 KEGG 富集分析。这样可以快速判断模块的功能主题与目标性状是否吻合。比如在植物抗逆项目里，与耐旱性状正相关的模块富集出脱落酸信号通路，那这个模块的生物学可信度就很高。

接下来是筛选模块内的中心基因。常用指标包括 kME 值、模块内连接度（intramodular connectivity）以及与其他模块基因的共表达关系。挑选 kME 较高且蛋白互作网络中处于核心位置的基因作为候选，这类基因更可能在模块中发挥调控作用。最后再用独立数据集或者 qPCR 实验做一轮验证，整个流程就闭环了。

6.2 一些值得坚持的工作习惯

我做了不少 WGCNA 相关项目后，最大的体会是：这个分析流程本身不难，难的是每一步都留下可复现的记录。建议你从一开始就建立一个分析日志文件，记下表达矩阵的标准化方式、软阈值的选择依据、模块参数设置、合并阈值，以及每次调整后模块数量变化。这样即使过了一个月回头看，也能快速还原分析过程。

版本管理同样重要。同一个项目，我不会直接覆盖原来的图，而是按日期或者参数命名，比如module_trait_power6_min30.pdf、module_trait_power8_min20.pdf。比较不同参数下的结果时，这种命名方式能省很多时间。还有一点，出图时顺手把 ME、模块颜色、基因分配表都导出成 CSV，这样后续做分析、写文章的时候不需要重新跑一遍代码。

模块与性状展示图只是 WGCNA 分析中的一个节点，但把这张图做扎实了，后面的路会顺畅很多。最后分享一个小技巧：在最终提交论文前，把这张图和样本聚类树放在同一个补充图里，审稿人可以直观地看到样本关系与模块的对应情况，这一细节虽然不起眼，却经常能减少一轮返修中的质疑。