WGCNA 实跑复盘，强负相关的基因为什么会进入同一个模块

看到两个基因的相关系数接近负一，却被 WGCNA 分到同一模块，很多人会怀疑程序分错了。实际上，是否保留相关方向，是建网时必须明确的选择。本次让两组人工基因分别跟随同一个变化趋势及其相反方向，实际运行两种网络，得到截然不同的归组结果。

这篇文章围绕这一组对照解释相关系数、连接权重、模块成员和模块特征向量之间的关系。同时运行输入质量检查、一个异常点相关性练习和一次样本子集对照。所有数据均为模拟表达值，不是某种组织的真实发现，也没有把构造出来的高相关当作生物学验证。

输入矩阵的方向先检查清楚

本次安装 WGCNA 1.74，输入按样本为行、特征为列组织。有效数据有四十个样本、一百四十个特征。甲组四十个特征跟随第一个潜在趋势，乙组四十个特征跟随它的相反方向，丙组四十个特征跟随另一个趋势，剩余二十个是独立噪声特征。

模拟值是连续的表达型数据，可以出现负值。它们不是原始测序计数，也没有经过真实数据的标准化或方差稳定变换。因此，本实验只检验建网和结果解释，不验证从原始读段到表达矩阵的准备过程。将自己的原始计数直接代入这里的示例，并不自动获得相同的输入条件。

为了检验入口检查，我们额外加入一列常数、一列全缺失，以及一个全部缺失的样本。实际调用质量检查后，这三项被排除，输入从四十一行、一百四十二列变成四十行、一百四十列。检查结果与预先安排的异常一致，说明流程确实处理了这些无法提供有效相关信息的项目。

这类检查不等于样本质量已经全面合格。一个样本即使没有缺失值，也可能存在批次影响、错标或异常分布。质量检查给出的通过，只对应其检查项目。正式分析仍应查看样本关系和实验记录，而不是用一个全部通过的布尔值替代数据理解。

qc <- goodSamplesGenes(dirty, verbose = 0)
clean <- dirty[qc$goodSamples, qc$goodGenes]
stopifnot(identical(dim(clean), c(40L, 140L)))

相关系数和连接权重不是同一个数

在有效矩阵中，甲组第一个特征与乙组第一个特征的相关系数约为负零点九六八，符合我们安排的相反变化趋势。随后分别计算无符号网络和有符号网络的连接权重，幂都固定为六。这个幂值只是本次控制变量的设置，没有通过软阈值选择过程优化，也不作为通用推荐值。

无符号网络使用相关强度的绝对值，因此强负相关仍然对应强连接。有符号网络保留方向信息，同样的强负相关被转换为接近零的连接。两种定义服务于不同问题，不能仅根据哪一种得到的模块更多就评价优劣。相关网络和加权连接的基本思想见 WGCNA 方法论文。

第一行相同，后两行不同，是因为变换规则改变了，并非重新观察到了另一份表达数据。尤其不能把连接权重零点八二解释为两者有百分之八十二的相互作用概率。它是当前定义和参数下的网络权重，含义需要连同计算规则一起报告。

相同幂值也不保证两种网络中的权重可以逐项当成同一尺度解释。表中的正相关连接在两种定义下也有差别。比较网络时，应首先明确希望保留哪类关系，再在该定义内评估参数，而不是把某个数字较大直接当作证据更充分。

两种定义实际得到了怎样的模块

我们在同一矩阵上运行模块识别，最小模块大小设为二十，合并阈值设为零点二五，相关方法均使用 Pearson。两次分别指定网络类型及对应的拓扑重叠类型。其余主要条件保持一致，以便把结果差异与相关方向的处理联系起来。

无符号网络将甲组和乙组合并为一个八十特征的模块，丙组形成四十特征的模块，二十个噪声特征未被分配。有符号网络将甲、乙、丙分成三个各四十特征的模块，同样留下二十个未分配特征。人工构造的相反变化模式，在这两套定义下得到的归组正好不同。

未分配类别在本次数字标签中记为零，常被显示成灰色。它不应和识别出的模块一起，直接计入具有同等含义的模块数量。若报告里写有三个颜色类别和三个模块，读者可能误以为灰色也是一个同样紧密的共表达群体，实际含义并非如此。

模块编号和颜色也不是跨运行稳定的身份标识。本次甲组在两次运行中的编号不同，并不说明基因名称改变。比较两份结果，应根据成员关系或特定映射核对，而不是把第一份结果中的蓝色，与第二份结果中的蓝色直接视作同一生物学对象。

模块与性状正相关，不代表所有成员同向变化

在无符号网络的八十特征模块中，模块特征向量与第一个人工趋势几乎完全正相关。但是甲组第一个特征与该向量的相关约为零点九九零，乙组第一个特征与该向量的相关约为负零点九八七。两者都在同一个模块，却处在相反方向。

这给出一个很具体的反例。若仅看到模块与某个性状正相关，就把模块内所有基因都写成随性状升高而升高，会错误解释其中一半人工特征。模块的摘要关系不能替代逐基因方向检查。特别是使用无符号网络时，成员可能围绕同一变化轴呈相反趋势。

模块特征向量本身还有方向约定的问题，不能脱离当前输出将符号当作固定生物学标签。更稳妥的做法是查看代表特征、成员与摘要向量的相关，以及表达模式图，并把比较方向写清楚。数值标签只是整理工具，实际解释应回到数据中的变化关系。

本例中模块与人工趋势的相关接近一，是因为表达值本来就按该趋势生成。它证明代码找回了设计中的结构，不是独立发现了一个强生物学关联。真实项目里的性状若与批次、组织比例或其他因素相关，还需要进一步分析这些可能共同驱动表达的来源。

一个异常点能把相关方向翻过来

我们另做了一个独立的双变量练习，未将它混入上述模块数据。先让四十对数值沿同一方向变化，再把最后一对改成一个变量特别大、另一个特别小。实际计算的 Pearson 相关约为负零点八四三，稳健相关约为正零点九九六。

这说明一个高影响观测可以主导普通相关，甚至使方向与其余大多数观测相反。但不能据此把所有与总体模式不同的点都删除。它可能来自技术异常，也可能来自真实的亚群或双峰分布，是否应当降低影响需要结合研究问题判断。

本次稳健相关只用于这一对变量，没有重跑完整的稳健网络。因此，文章不会声称稳健方法已经修复了所有模块，也没有据此比较两套完整网络的优劣。实现和方法背景可参考 稳健相关函数论文。

stats::cor(x, y)
bicor(x, y, maxPOutliers = 0.1)

作者的 常见问题说明 对稳健相关在双峰数据中的行为有专门提醒，也讨论了限制异常比例的设置。该说明还强调小样本相关的不稳定性。选择相关方法前，先看数据分布及样本结构，比仅把方法名称从普通改成稳健更重要。

子集运行说明了什么，又没有说明什么

为了查看样本变化，我们预先指定取前面三十个样本，重新运行有符号网络。结果仍有三个主要模块，但其中一个模块变成四十一个特征，未分配数量从二十降到十九。具体变化是噪声特征 N16 从未分配类别进入了一个模块。

这个例子说明，整体模块结构看起来相近时，边界成员仍可能改变。对于准备挑选后续实验候选的研究，不能只比较模块总数，还应查看候选是否在不同合理设置下保持同样关系。若候选恰好处于归组边缘，报告中应保留这种不确定性。

但我们只运行了一个固定子集，没有进行反复抽样、置换检验、独立队列复现或模块保存性检验。结果记录中的配对归类一致率还包含未分配标签，不能当作标准模块保存性统计量。因此本文不把一次子集结果写成网络已经稳定验证，而只把它作为检查入口。

需要正式评价稳定性时，应预先确定评价目标。是关注模块整体是否保留，还是关注某个基因的成员关系，或者关注模块与性状的关联是否重现。不同问题需要不同证据，不能用一个很高的一致率替代所有层面的验证。

参数控制与可复现运行

完整脚本见 本组实跑资料。驱动程序在固定环境中安装软件，生成输入，运行质量检查、两种网络和子集对照，并保存表达矩阵、人工趋势、模块成员及特征向量。进入 WGCNA 对应目录后执行以下命令即可重建本实验。

python3 run.py run160

下面是有符号网络对照中的主要设置。它用于说明本次实际运行条件，不是可以不看数据就使用的最佳参数组合。全部参数及随机种子保存在驱动脚本中，便于检查每次修改究竟改变了什么。

net <- blockwiseModules(
  clean, power = 6,
  networkType = "signed", TOMType = "signed",
  corType = "pearson", minModuleSize = 20,
  mergeCutHeight = 0.25, numericLabels = TRUE,
  reassignThreshold = 0, pamRespectsDendro = FALSE,
  maxBlockSize = 200, randomSeed = 160, verbose = 0
)

本次特征数量小于分块上限，实际没有用这个实验评价大规模分块建网。真实分析中，稠密关系矩阵的大小会随特征数平方增长。二万个特征的一份双精度方阵，仅数值部分约需三点二吉字节，还没有计入其他矩阵和对象副本。这是维度推算，不是本次测得的内存峰值。

因此，扩大分析前应先明确特征过滤与分块策略，避免把所有中间矩阵都无条件保留。另一方面，也不宜为了得到更明显的模块，先按同一个性状筛出显著差异基因再建网。那会改变输入的相关结构，使后续与该性状的关系更难作为独立证据解释。

优点、局限和实际选用

WGCNA 的实用优势是能够把大量特征的共同变化组织成模块，再以模块层面的摘要连接到研究问题。本次两种网络的对照说明，它允许使用者明确选择是否将相反方向的变化视作相近关系。可检查的成员表和特征向量，也方便把概括性结论追溯到具体特征。

它的局限在于相关结构容易受样本组成、异常观测和参数影响。模块并不自动等于通路，连接也不自动等于调控关系。一个模块中挑出的枢纽特征，仍需要根据表达质量、稳定性和独立证据决定是否值得后续验证，不能仅凭连接度较高就称其为驱动因素。

运行成功之后，最应补上的不是更多颜色，而是对网络定义和成员方向的解释。强负相关为什么被分到一起，模块关联是否被少量观测驱动，候选是否处于归组边界。把这些问题说清楚，模块结果才真正能够支持下一步研究判断。

当明确的建网方案需要扩展到更多队列或重复分析时，可结合特征规模了解 云生信生物信息学分析服务器服务。