企业官网建设流程全解析

WGCNA这名字在转录组文章里出现得太频繁了，以至于刚接触生信的人大概率都会问一句：它到底解决了什么问题？我当年的理解是，差异分析告诉我们哪些基因变了，而WGCNA告诉我们哪些基因是“一起变”的。后来真正跑完一遍才发现，这句话只对了一半。WGCNA（Weighted Gene Co-expression Network Analysis，加权基因共表达网络分析）核心价值是把成百上千个基因归类成少数模块，再把模块和样本性状关联起来，从中锁定期值得往下走的候选基因。这篇文章就从零开始，把数据格式、软阈值选择、网络构建、模块识别、性状关联到候选基因筛选的整个流程掰开揉碎讲透，配合可直接照抄的R代码和我在真实项目里踩过的坑，适合刚入门转录组分析、手上有一批表达谱数据却不知道下一步该怎么做的同学。

1. WGCNA到底在干什么——从“差异基因”到“基因社交圈”

1.1 一个直观的类比：基因也有自己的“朋友圈”

想象一下你把一个班级的人按作息习惯聚类：有人习惯早起，有人通宵赶due，喜欢熬夜的人经常同一时间出现在自习室，早起的人总在食堂同一桌碰面。WGCNA做的事情和这个非常相似——它会扫描每个基因在不同样本里的表达曲线，把表达趋势高度相似的基因归进同一个模块，这些模块在文献里也被称为基因共表达模块。

和传统的差异分析相比，WGCNA解决的是完全不同的问题。差异分析回答“谁发生了变化”，比如疾病组vs对照组里有158个基因上调、97个下调。但到了这一步，你得到的仍然是一长串基因列表，几千个基因在通路富集里堆在一起，根本不知道谁在上游、谁跟谁协同工作。WGCNA则把问题变成了“哪些基因像一个班级里的小团体一样总是一起行动”，以及“这些小团体里谁是话语权最大的核心人物”。

在实际项目中，这两步通常是串联关系：先用DESeq2或limma筛出一批显著差异基因，再用这批基因构建WGCNA网络，把富集分析和模块分析结合成一条证据链。这样做的逻辑是，模块本身可能包含很多差异基因，也可能包含差异不显著但和性状相关性很强的基因——后者往往是新手最容易忽略、但实验验证价值同样很高的部分。

1.2 核心假设：生物网络接近“无尺度网络”

WGCNA这个名字里的“Weighted”是整个方法论的灵魂。很多人的第一反应是：既然要研究基因之间的关系，直接计算两两基因之间的相关系数，然后把相关系数大于0.8的基因对连接起来不就行了？

这就是硬阈值网络的思路，但现实中它有两个致命问题：一是阈值怎么定，0.7和0.8画出来的网络可能天差地别，结果对参数极其敏感；二是硬截断会丢掉相关系数0.79、0.76这类边缘信息，在转录组这么大的噪声背景下，这种取舍非常危险。

WGCNA的思路是加权：基因间相关系数r的绝对值经过一个幂指数变换后直接作为边的权重，a_ij = |cor(i,j)|^β。β越大，权重之间的差异越明显，弱相关被压得更低，强相关则保留下来，整个网络会逐渐呈现出一种特殊的拓扑结构——少数基因连接大量基因，大多数基因只连接少数基因，这就是无尺度网络（scale-free network）。

听起来很抽象，但你可以把它类比成航空网络：少数几个大枢纽机场串联全球航线，大量小机场只飞往就近的枢纽。生物体系里的调控关系同样如此，少量转录因子负责切换大范围的基因表达程序，大量功能基因是“乘客”。WGCNA选择β的核心标准，就是让基因网络尽量接近这种无尺度结构。理解了这一点，你就理解了软阈值为什么那么重要，也理解了WGCNA为什么能比普通相关热图更接近真实的生物学调控网络。

2. 数据准备：别让第一步成为最大的坑

2.1 表达矩阵应该长什么样

WGCNA的输入是一张表达谱矩阵，通常行是基因、列是样本。但在R里跑blockwiseModules时，它内部默认的数据格式是样本在行、基因在列，所以读入原始数据后第一件事往往是转置。

library(WGCNA) options(stringsAsFactors = FALSE) allowWGCNAThreads() # 多核运行，省时间的关键 expr <- read.csv("expression_matrix.csv", row.names = 1) # 假设原始数据是基因行、样本列，转换成WGCNA要求的样本行、基因列 datExpr <- as.data.frame(t(expr)) datExpr[] <- lapply(datExpr, as.numeric) # 检查缺失和零方差 gsg <- goodSamplesGenes(datExpr, verbose = 3) if (!gsg$allOK) { datExpr <- datExpr[gsg$goodSamples, gsg$goodGenes] }

这里有个容易被忽略的细节：表达量需要先做好标准化。如果是RNA-seq的read count，请先转成log2 CPM或VST表达量，直接拿原始整数count跑相关系数会被文库大小差异干扰得面目全非。芯片数据则通常用log2后的表达值，一般问题不大。goodSamplesGenes这个函数会帮你识别缺失值和方差为0的基因，但它只是兜底检查，真正的数据清理工作应该在读入之前完成。

2.2 缺失值、离群样本与批次效应

WGCNA官方文档里明确说，缺失值会让相关系数计算产生偏差，因为cor函数默认按pairwise complete来处理缺失，这导致不同的基因对用了不同的样本子集，模块结构会被扭曲。处理策略我给出两条经验：缺失比例超过5%的基因直接删除；缺失比例低的基因，可以用impute.knn这类方法填补，但填补后的结果一定要比较一下聚类树是否有异常。

样本离群比缺失值更隐蔽。一个简单实用的检查方法是先对所有样本做层次聚类，画一张树状图：

sampleTree <- hclust(dist(datExpr), method = "average") pdf("sample_tree.pdf", width = 12, height = 6) plot(sampleTree) dev.off()

如果某个样本孤零零挂在树的远端，高度明显高于其他样本之间的连接高度，那它很可能是质量问题，比如RNA降解严重、测序量过低、标签弄错。对待离群样本要冷静：如果样本数本来就少（20个左右），删掉一个影响就很大，先检查是否可以通过批次校正救回来；如果样本数足够多，剔除明显离群样本会让后续软阈值达标率明显提升。

批次效应的处理是另一个大坑。多批次测序、不同时间提取RNA、不同实验员操作都会在表达谱里留下系统性信号，不校正的话，WGCNA很容易把“批次”当成“模块”，因为同一批次的样本表达谱天然相似。常用的做法是先用limma包里的removeBatchEffect或者sva包的ComBat_seq做校正，然后再进入WGCNA流程。校正之后务必重新做一次样本聚类，确认批次不再是树上最大的分组因素。

2.3 选哪些基因进入网络：全基因组还是差异基因

很多人第一次跑WGCNA，兴奋地把两万个基因全部丢进去，结果跑了两小时，模块切出来稀碎，后期富集分析更是无从谈起。问题的根源在于：表达变化极小的基因有大量随机噪声，它们会让相关系数下降，还会干扰模块边界。

基于我在项目中的经验，更推荐两种基因筛选策略。第一种：先用DESeq2或limma筛差异基因（比如padj<0.05且|log2FC|>1），然后用这些基因建网络。这种策略下模块更容易富集到有生物学意义的功能簇，因为你和表型关联的搜索范围已经高度聚焦。第二种：如果你不想预设差异阈值，可以按方差排序，取表达标准差前5000到10000个基因，这些基因提供了全基因组层面的共表达信息，适合探索性分析。

权衡标准很现实：全基因组网络样本量需求大、计算量大、模块往往太多；差异基因网络模块少、解释容易、和性状关系更直接。我的习惯是先拿差异基因跑一遍，把模块性状关联看出来之后，再考虑要不要扩大到全基因组。你还要留个心眼：同一项目里如果已知存在某些调控核心基因，即使它们差异不显著，也建议手工保留进基因集，否则模块里的hub gene可能并不完整。

3. 软阈值是WGCNA的灵魂，理解了它你就懂了一半

3.1 为什么必须用加权网络而不是硬阈值

前面已经简单提过硬阈值的问题，这里再展开说透。硬阈值方式大概是这样的逻辑：先算所有基因两两之间的Pearson相关系数，再把相关系数超过0.85的基因对保留为“网络连接”，其余一律视为无连接。看起来直观，但相关系数本身的分布受样本量、表达动态范围、批次噪声的影响很大，0.85在不同数据集里的意义完全不同。

更麻烦的是，硬阈值会破坏无尺度拓扑。想象一个基因表达波动比较平缓，它和大量基因的相关性都徘徊在0.75到0.85之间，硬阈值一设，这些潜在连接要么全部保留、要么全部删除，网络结构极其脆弱。而WGCNA的软阈值通过幂指数加权，让相关系数0.9的基因对获得接近1的权重，相关系数0.5的基因对权重变得极小，既不会误连低相关对，又能保留弱相关的梯度信息。在实际数据里，加权网络的模块重复性和生物学可解释性都比硬阈值网络强得多，这是它在2010年前后迅速成为转录组分析标配的重要原因。

另一个要点是networkType的选择。WGCNA提供signed和unsigned两种网络类型，unsigned取相关系数的绝对值，正相关和负相关都算连接；signed只看正相关，负相关基因对被看作不连接。大多数入门教程默认用unsigned，但我建议你在正式跑之前想清楚自己的生物学问题：如果关注的是功能协同单元，signed往往更干净，因为负相关的基因很可能处于相互抑制的关系，放在同一个模块里解释起来比较别扭；如果关注的是全局调控关系，unsigned能保留更多信息。模块数量也会受这个选择影响，切换后要重新评估所有结果。

3.2 pickSoftThreshold：R²到底该怎么读

选定软阈值的核心函数是pickSoftThreshold，它会遍历一组候选的power值，评估每个power下构建的网络有多接近无尺度拓扑，并返回平均连接度等指标：

powers <- c(1:20) sft <- pickSoftThreshold(datExpr, powerVector = powers, networkType = "unsigned", verbose = 5)

输出结果里最关键的是两列：SFT.R.sq和mean.k。SFT.R.sq代表无尺度拟合指数，反映当前网络和理想无尺度网络的接近程度，官方经验值是把它推到0.8到0.9以上。mean.k是所有基因的平均连接度，代表网络整体的连通水平。

选择power的基本原则是：找到能让SFT.R.sq首次达到0.8以上的最小power，并且观察此时mean.k是否已经开始断崖式下跌。如果选一个特别大的power，虽然无尺度拟合很好，但平均连接度会变得非常低，网络里只剩下少数枢纽基因，模块数量也会少得可怜；power选小了，网络会偏向随机关联，模块碎片化。实际操作时，我会把fitIndices打印出来，逐行看，不要只交给默认逻辑。

print(sft$fitIndices[, c("power", "SFT.R.sq", "mean.k", "median.k")])

常见的坑是：数据质量一般时，SFT.R.sq死活到不了0.8。这时候不要强行往下跑，先回头检查样本聚类，剔除离群样本，或者从unsigned切成signed、把Pearson相关改成Spearman相关试试。实在不行，可以选择SFT.R.sq曲线开始变平拐点处的power作为手动选择值——这是最常用的补救方式，但要清楚这是权衡之举，结果解读要更谨慎。

4. 网络构建与模块识别：从邻接矩阵到基因模块

4.1 TOM矩阵：为什么不能直接用相关系数矩阵聚类

选定软阈值之后，WGCNA会做三件环环相扣的事：先构建邻接矩阵，再把它转换成拓扑重叠矩阵（TOM），最后基于TOM距离做层次聚类和动态模块切割。

邻接矩阵本身只是一组加权后的相关系数，但它有一个先天不足：只描述两个基因之间的直接关系，不考虑它们在网络中的共同邻居。一个生活化的例子：A和B之间互动多、B和C之间互动多，如果只看两两关系，A和C可能看起来不熟，但在真实的社交网络里，A和C之所以认识，正是因为有B这个共同好友，他们其实处于同一个紧密社区。TOM矩阵做的正是这一件事：它把“两个基因共享多少个网络邻居”纳入考量，如果A和C共同连接了大量基因，那么即使它们的直接相关系数不高，TOM值同样会很高。

这就是为什么WGCNA要用TOM矩阵而不是相关系数矩阵做聚类。TOM让模块边界更平滑，显著减少了对单个噪声基因对（某个异常的相关系数）的敏感性。实际计算中，blockwiseModules会自动完成从邻接矩阵、TOM矩阵到层次聚类的全过程，你也可以手动分步跑，方便在中间步骤检查矩阵。

这里有个小提示：TOM矩阵的体积和基因数的平方成正比，一万个基因的TOM矩阵就是上亿个元素，内存占用非常可观。所以前一步的基因筛选不只是为了生物学解释，也是在为这一步的算力减压。

4.2 动态树切割和模块合并：参数怎么调才有意义

聚类完成后，面临的第一个问题是：树状图上怎么切分模块？传统方法是设定一个固定的高度阈值一刀切，但基因表达聚类树的各个分支高度差异很大，一刀切既可能把同一功能簇劈成几块，又可能把不相关的分支错误合并。WGCNA官方推荐使用dynamicTreeCut包里的动态切割算法，它根据树状图的局部结构特征决定切分位置，能识别出嵌套在深层的精细模块。

net <- blockwiseModules( datExpr, power = 12, # 这里填入你上面选的软阈值 networkType = "unsigned", TOMType = "unsigned", minModuleSize = 30, # 模块最少基因数 reassignThreshold = 0, mergeCutHeight = 0.25, # 模块合并阈值 numericLabels = TRUE, saveTOMs = TRUE, saveTOMFileBase = "geneTOM", verbose = 3 )

blockwiseModules把整个流程打包了，比手动一步步做省心很多，初学者建议直接用。关键参数有三个：minModuleSize控制模块最小基因规模，基因太少无法稳定估计共表达结构，默认30一般够用；deepSplit控制切割灵敏度，取值0到4，越大切出的模块越多越细，默认2，如果模块太粗可以调到3或4；mergeCutHeight控制模块合并，默认0.25对应的含义是：当两个模块的特征基因相关性超过0.75时，它们会被合并成一个模块。

模块合并这个步骤经常被轻视，但它对结果影响很大。聚类算法有时会把同一个生物学通路劈成两个勉强分开的分支，合并就是用来消除这种低水平噪声的。如果你在结果里看到几个模块的ME相关性明显很高、热图颜色又非常接近，那可能是mergeCutHeight设得太大或太小。判断标准只有一个：合并后的模块是否能更简洁地解释下游富集和性状关联。保存聚类树图时，模块颜色会标在树的底部：

moduleColors <- labels2colors(net$colors) table(moduleColors) # 看看每个模块的基因数

灰色（grey）模块永远是特殊的存在——它代表所有没有被分进任何模块的基因。分析模块时默认要排除灰色，因为它是“无组织”基因的垃圾桶，不具有共表达结构上的稳定性。

5. 模块与性状关联：从热图到候选基因

5.1 module eigengene、GS与MM：三个核心概念一次讲清

模块是基因集合，但和性状做关联时不能把模块里几百个基因逐个拿去算相关，那样会陷入多重检验的地狱。WGCNA的解决办法是提取模块的第一主成分，称为模块特征基因（module eigengene，ME），它浓缩了这个模块整体表达模式。从数学上讲，主成分就是样本维度上最能代表模块表达趋势的“虚拟基因”，你可以把它理解成一个集体发言人。

一旦得到ME，后面的计算就顺理成章了：模块与性状的关系（module-trait relationship）就是ME和性状向量的相关系数，相关性检验给出p值；基因显著性GS（gene significance）是单个基因表达谱和性状的相关系数；模块成员MM（module membership）则是单个基因表达谱和所在模块ME的相关系数。

MEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes nSamples <- nrow(datExpr) datTraits <- read.csv("traits.csv", row.names = 1) # 样本行、性状列 moduleTraitCor <- cor(MEs, datTraits, use = "p") moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples)

MM在这里特别有用，因为它是判断hub gene的核心依据。一个模块里MM的绝对值越接近1，说明这个基因的表达趋势和模块整体越一致，越可能是模块的核心执行者。实际操作中，我通常把|MM| > 0.8作为候选hub gene筛选条件，结合模块内连接度kWithin排序，两者都靠前的基因基本就是你需要重点盯住的对象。要注意的是，MM是一个有符号的值，正负号代表该基因与模块ME是同向还是反向变化，两者方向相反时生物学含义完全不同。

5.2 模块-性状关联热图的解读技巧

热图是WGCNA结果里最直观的一张图，它展示每个模块ME和每个性状间的相关系数与p值。绘制代码是：

textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = "") labeledHeatmap( Matrix = moduleTraitCor, xLabels = names(datTraits), yLabels = names(MEs), colorLabels = FALSE, colors = blueWhiteRed(50), textMatrix = textMatrix, cex.lab = 0.8, main = "Module-trait relationships" )

看懂这张热图有两个层次。第一层是找“显著模块”：看哪个模块的p值小于0.05、相关系数绝对值较大，那个模块就是和该性状关联最密切的候选模块。第二层是把相关系数的符号和生物学知识对照起来——模块ME和患病状态强烈正相关，意味着模块整体在患病组上调；如果有人把模块当成“疾病上调模块”去解读，却忘了检查具体是哪个方向，很容易得出相反的结论。

对于性状类型，还要做一点区分。连续型性状（药物浓度、年龄、BMI）用Pearson相关检验没有任何问题；二分类性状（患病vs对照、治疗vs未治疗）虽然也可以用相关分析近似，但更严谨的做法是把分组作为因子，对模块ME做t检验或线性模型。两种做法的p值在大多数情况下结论一致，但遇到边界情况时，t检验给出的p值更可靠。

选择候选模块时还应该引入另一个指标：模块显著性（module significance），它是模块内所有基因GS的均值。如果一个模块内基因普遍和性状相关，说明这个模块不是靠少数极端基因撑起来的，而是整体性地响应性状变化，这样的模块更值得投入下游验证。

6. 常见问题与实操心得：这些坑我都替你踩过了

6.1 问题速查表

把我在真实项目里遇到过的问题整理成一张表，遇到类似情况可以直接对号入座。

这里想多强调一句：如果你看到“全不显著”，先别急着怀疑算法。回顾一下样本数——WGCNA在样本少于15个时模块结构会非常不稳定，少于10个基本没有可信度；回顾一下性状分组——有些表型是连续变量但被误存成字符型，读入R后变成了因子，相关性计算直接报错或给出无意义结果。这些都是我帮人排查过很多次的低级但高频问题。

6.2 实操心得：为什么你的结果和别人不一样

第一，多线程务必打开。allowWGCNAThreads()可以让你在大型数据集上省下几个小时，默认单线程跑一万个基因会让你怀疑人生。

第二，基因名重复要提前处理。同一个基因对应多个探针或isoform时，不处理会出现完全相同的行，看起来相关是1，实际上是冗余信息。建议先取表达均值或方差最大的探针作为代表。

第三，保存工作空间。save.image("wgcna_step1.RData")这种习惯能救你的命。软阈值计算、TOM矩阵构建都是耗时步骤，一旦R崩溃，没有进度存档就只能从头再来。

第四，模块特征基因的符号问题。在筛选hub gene时，别只取MM的绝对值。一个MM为-0.9的基因和模块ME呈强负相关，它可能是在模块核心表达程序里扮演抑制作用的关键角色。如果你把它和MM为0.9的基因混在一个列表里，下游富集和实验验证的逻辑都会乱掉。

第五，可视化时导出Cytoscape格式做网络图。exportNetworkToCytoscape可以把指定模块的边权重导出，筛选权重前5%的连接，导入Cytoscape画出的网络图比静态热图更能展示hub gene的地位。这一步在论文补充材料里非常加分。

第六，灰色模块永远不要纳入hub gene分析。灰色是未分配基因的集合，不具备共表达模块的生物学意义，有些人图省事把灰色里相关性好的基因拉出来当hub gene，答辩时被问住很尴尬。

第七，结合差异表达和富集分析一起讲。WGCNA自己不能做功能注释，模块鉴定出来后马上用clusterProfiler做GO/KEGG富集，把模块里的基因集合和通路对应起来，才算走完从数据到机制假设的闭环。

6.3 后续还能往哪个方向走

这篇文章的定位是从零基础开始跑通主流程，但WGCNA的应用远不止于此。跑熟单数据集后，你可以试试多数据集整合的共识网络分析blockwiseConsensusModules，它解决了跨批次、跨平台数据不一致的问题；也可以对模块内基因做蛋白-蛋白互作网络叠加分析，把共表达和物理互作两个维度结合起来找更加可靠的调控核心；还可以把hub gene的表达趋势放进时序数据里做动态轨迹分析，看它是在早期、中期还是后期起主导作用。

这些方向的前提都一样：把单数据集WGCNA的基本盘打扎实，明白每一步在做什么，而不是只知道调包跑完出图。

我个人在实际操作中体会最深的一点是，WGCNA不是一个“一键出结果”的黑箱工具，它更像一个需要你和它反复对齐认知的分析框架。软阈值选多少、模块合并到什么程度、hub gene的阈值定多严，都和你手头数据的质量、样本量、生物学背景密切相关。套默认参数跑完一次很容易，但要让模块和性状的关联经得起推敲，要在参数选择时多问自己一句“为什么是它”。最后再分享一个操作习惯：分析过程中的聚类树图、软阈值曲线、模块颜色分布、ME热图，每一张图都存好并记录当时用的参数组合。我吃过亏，跑完一堆分析回头看时，已经不记得某个模块是在哪个参数组合下得到的，重跑一遍的成本高得离谱。数据时代，记录和结果本身同样珍贵。