WGCNA原理详解：从共表达网络到基因模块识别的完整思路

我做了快十年的生信分析WGCNA这个工具几乎是每个做转录组、表达谱芯片、单细胞下游分析的人都会撞上的坎。第一次听到“加权基因共表达网络分析”这个名字时我脑子里只有两个字——劝退。但真把这个工具吃透之后你会发现它其实是一条特别清晰的思路从一堆基因里找到协作小组再把小组和表型挂钩。这篇文章先把最核心的“简介、原理”讲明白后续再聊具体怎么跑、参数怎么调。我会尽量避开教科书腔用做项目时积累下来的理解和踩坑经验把这套方法背后的逻辑一层层剥开。不管你是刚进入生信领域的新人还是已经跑过几轮WGCNA但总觉得哪里没理顺的老手这篇文章应该都能帮你把整条链路打通。1. WGCNA到底解决什么问题——从“单个基因”到“基因团队”思维的转变1.1 你手里那堆差异基因真的够用吗先说说我自己的经历。早些年做转录组分析拿到一批差异表达基因之后常规操作就是往GO、KEGG里一丢看看富集到哪些通路挑几条画个通路图然后开始编故事。这事本身没问题但慢慢你会发现有几个特别头疼的场景差异基因数量动辄上千个富集结果里全是“代谢通路”“免疫应答”这类大而空的条目根本看不出重点。不同批次实验或者不同组织之间差来差去总是那一批经典通路的基因有价值的新发现被淹没掉了。有些基因单独看表达量差异不大p值也确实不显著但它在整个基因调控网络中可能是承上启下的关键节点这种基因用差异表达的思路永远捞不到。WGCNA处理问题的角度完全不同。它不关心基因单独变化了多少倍它关心的是在整个样本群体里哪些基因的表达模式总是同步起伏。如果一组基因在样本A里集体上调在样本B里集体下调步伐高度一致那么这群基因很可能属于同一个生物学程序——比如同一个通路、同一个细胞类型的功能模块、同一类调控事件的下游响应。这就是“共表达模块”的基本概念。你说这不就是以基因为变量算相关性吗对思路就这么朴素但落地的时候数学细节就多了后面慢慢展开。1.2 为什么“加权”是这套方法最得意的设计标准的相关性网络很简单两个基因表达量的Pearson相关系数超过某个阈值就连一条边没过阈值就不连这是个二值化操作要么有边要么没边。WGCNA的“加权”就体现在它不这么简单粗暴。它对相关性做了幂次变换公式大概是这样的连接强度 |cor(gene_i, gene_j)|^ββ是软阈值不是固定为1或者2而是在分析时专门选择的。这种处理方式的影响我给你打个比方把基因想象成社交网络里的人。硬阈值就像是“只要你们两个没有说过话就不算朋友”而软阈值则是“就算你们只是偶尔见过面也会有一点交情只是亲疏程度不同”。这样一来每条边的强度是连续的既有强联系也有弱联系整个网络的拓扑结构就丰富多了。后面讲原理的时候我会重点解释为什么这个β值通常在5到20之间调来调去以及“无标度网络”又是怎么回事。1.3 WGCNA不是万能的——它适合什么不适合什么先泼一盆冷水WGCNA不是所有场景都能用。它最典型的使用场景是你手头有一批样本的表达谱数据样本数不能太少我个人建议至少15个以上越多越稳而且你对样本的某个属性特别关心。比如肿瘤不同分期、不同恶性程度的表达谱想找与分级相关、可能驱动恶性转化的基因模块。不同发育阶段的时序转录组想找在特定时间点协同变化的基因程序。不同药物处理浓度或不同处理时间点想找和药物剂量/时间相关的响应模块。它不太适合的场景也有几个样本数太少比如每组3个生物学重复一共也就9个样本这种数据算出来的相关性噪声太大。只关心少数几个明星基因不关心基因间的协同关系。数据本身是单细胞测序产出直接拿原始表达矩阵跑。单细胞数据零膨胀严重、稀疏性高不预处理跑WGCNA会有大量假模块通常需要先做伪bulk或者做基因模块分析的其他替代方案。理解边界比理解工具本身更重要。这就像你不能用螺丝刀钉钉子工具本身再好也要用对地方。2. 核心原理拆解一——共表达网络与无标度拓扑2.1 共表达关系是怎么算出来的咱们从最基础的步骤开始。假设你有一份表达矩阵行是基因列是样本。第一步通常是过滤掉低表达基因比如在所有样本里表达量都很低、基本检测不到的基因然后对表达量做标准化常见的做法是取log2让数据分布更接近正态。接下来任意两个基因i和j在整个样本集合里算Spearman或Pearson相关系数。我们得到的是一个基因×基因的相关矩阵SS[i][j]就是基因i和基因j的共表达强度取值范围是-1到1。到这一步你可能会问负相关怎么办WGCNA里有两种设计选择取绝对值|S[i][j]|这种做法把正负相关统一当成“关联强度”因为负相关也可能是生物学上有意义的关系比如一个基因激活另一个被抑制。不取绝对值用signed类型把相关系数转换为0到1之间的数正相关接近1负相关接近0。实际项目里我一般用signed网络因为对生物学解释更友好。比如两群基因一个正调控一个负调控它们之间的强烈负相关关系在unsigned网络里会被当成强关联连在一起模块会糊成一大片。2.2 无标度网络——为什么“少数基因掌握话语权”反而成了判断标准这是WGCNA里让很多人头疼的概念scale-free topology无标度拓扑。自然界的真实网络经常呈现出一种所谓“无标度”特征。通俗地讲就是网络里绝大多数节点只有很少的连接而极少数节点连接了几乎整个网络。一个典型例子是机场网络——大量小机场只通一两条航线而少数枢纽机场连接了几十个上百个城市。生物学里的代谢网络、蛋白互作网络也有类似特征。WGCNA的创始人Barabási等人的思路是如果真实的基因调控网络具有无标度特征那我们构造的共表达网络也应该尽量贴近这种特征。就像测量物理量要校准仪器一样WGCNA要选择合适的β让网络拓扑结构尽量近似这种无标度形态。具体操作上分析过程是这样的对于候选的β值从1到20左右每次选定一个β值计算任意两个基因的邻接矩阵a[i][j] |cor(i,j)|^β。把每个基因节点的连接度算出来即该基因与其他所有基因的连接强度之和。统计不同连接度分布的情况然后用线性回归拟合log(p(k))与log(k)之间的关系。看看回归的R²有多大也就是网络的“无标度拟合程度”。如果R²接近0.8或以上说明这批候选β值构造的网络已经比较接近无标度。换个角度说这时网络已经拉开了层级少数hub基因连接度很高大部分基因只有少量连接信息传递有明确的主干道。我见过不少人看不懂这个步骤总觉得R²越高越好。其实标准没那么死板对很多真实生物学数据来说R²到不了0.9能稳定在0.8以上就差不多可以接受了。另外还有一个实用经验尽量选最小的、能让R²达到0.8以上的β值。原因很简单β越高对相关性的压缩越狠动态范围越小容易丢失比较弱的共表达关系导致模块太少。这里额外说一点很多人不知道WGCNA包里有自动选择软阈值的函数也就是pickSoftThreshold。它会把你需要的1到20的候选β全部跑一遍输出一个结果图让你自己判断选哪个。在实际项目里我通常是根据它的推荐值再结合自己数据的模块数量和稳定性做微调标定之后遇到RNA-seq数据基本落在7到16之间。2.3 邻接矩阵与拓扑重叠矩阵——为什么不能直接用相关系数聚类选完β之后我们有了邻居矩阵。这时有个问题用邻接矩阵直接聚类行不行可以但不是最优。比如基因A和基因B表达相关性很高基因B和基因C相关性也很高但A和C之间相关性一般。这时候如果直接对邻接矩阵做聚类A和B、B和C可能被分到不同模块而实际上它们共享的调控机制可能是一致的。这就是所谓“间接关系”的问题。WGCNA的解决办法是引入拓扑重叠矩阵。TO思想的本质是两个基因之间有多相似不只取决于它们自己的直接联系还取决于它们是否共享了很多共同邻居。类比一下社交网络你和某人是朋友关系这是直接联系。你们两人都在同一个百人群里而且有80个共同的熟人说明你们的关系圈高度重叠就算你俩直接交流不多也大概率是同一个圈子的。拓扑重叠度高就意味着两个基因在整个网络的局部拓扑语境中扮演相似的角色属于同一个“圈子”的成员。TO的计算公式在WGCNA文档里有详细描述代码层面不需要你手写包里的TOMsimilarity函数直接搞定。这里重点理解它的生物学含义就行。我在实际项目中的体会是TO矩阵做聚类比直接对相关性矩阵聚类稳定得多模块很少因为个别基因的表达噪声而碎掉这是WGCNA设计里很核心的一环。3. 核心原理拆解二——模块识别、模块特征基因与模块-性状关联3.1 模块识别从TO距离到基因树状图再到动态剪枝有了TO矩阵之后WGCNA把TO矩阵转换为距离矩阵然后做层次聚类得到一棵系统进化树风格的聚类树。树上的每根大枝、小枝就对应一组表达模式相近的基因。但这里有个经典问题聚类树该切在哪个高度如果切太高模块太少、太粗切太低模块碎成一地很难解释。WGCNA默认用的动态树切割算法不像传统的固定高度切割那么死板它会根据树状图的拓扑结构自动找出稳定的分支作为模块。实际过程中有几个参数需要你注意minModuleSize模块最小基因数默认是30。如果你的数据量不大这个值可以调到20甚至10否则模块数量会非常少。但如果盲目调低得到一堆只有三五个基因的小模块后面生物学解释也不容易。deepSplit取值0到4数值越大切得越细模块越多。我一般从2开始试然后看模块数量是否合理再适当调整。mergeCutHeight模块合并阈值。有时候初始划分会出现两个模块它们的特征基因相关性极高比如高达0.9这个时候要思考一下它们其实是同一个程序的两个侧面。合并阈值默认0.25意思就是模块特征基因相关度高于0.75的两个模块会被合并。经验值给到这里样本量20到30个、基因数1万到2万时通常能得到15到40个模块不等每个模块几十到几千个基因。如果你得到了上百个模块大概率是参数选得太激进或者数据噪声太大。3.2 模块特征基因——用“一个值”代表“一群基因”的聪明设计模块划分完成之后每个模块里可能有几百上千个基因。你不可能把每个基因都拿出来和性状去统计那样又回到了“单个基因”的思维。WGCNA的处理方式是为主成分分析对模块内所有基因在样本维度上做主成分分析取第一主成分称为模块特征基因module eigengene简称ME。这可以理解成这个模块在某个样本里的“平均表达水平”的一种加权表达但它比简单平均更能代表模块内基因的整体表达变化趋势。ME是一个向量维度等于样本数代表该模块在每个样本里的综合表达高低。这样就把一群基因压缩成了一个“模块级别的表达量”。从这一步开始分析的语言就变了不再说“基因X在样本组A里上调”而是说“模块M在样本组A里整体激活”。从基因维度到模块维度的提升让后续统计量大大简化而且更容易和临床表型挂上钩。3.3 模块-性状关联找到你关心的那一个或那几个模块有了每个样本的模块特征基因表达值你就可以计算模块与性状之间的关联了。这里的性状可以是连续型变量如血压、肿瘤大小、生存时间也可以是二分组变量如患病组/对照组。WGCNA的标准输出里有模块特征基因与性状的相关性热图每一格就是模块和性状的相关系数和p值。格子越红越正相关越蓝越负相关。这里要特别提醒自己一个问题不要只盯着p值还得看相关系数的绝对值大小。一个模块与肿瘤分级相关系数是0.3p值0.01虽然统计上显著但生物学上这个相关强度只能算中等。0.6以上才算比较强的关联。当然这也取决于你实验设计本身的信号强度。另外一个常见操作是把模块特征基因的显著性gene significance即单个基因与性状的相关性绝对值和模块内连通性module membership即基因与模块特征基因的相关性结合起来。如果模块里的基因大多与性状相关性高而且与模块本身的相关性也高那这个模块就值得重点挖掘——里面的hub基因可能就是核心调控因子。4. 实操过程中踩过的坑——WGCNA使用的几个核心教训4.1 样本量真的不能太少关于样本量的经验国内国外论坛上帖子很多但最根本一点WGCNA是建立在一堆两两相关系数之上的每一对相关性估计都需要足够的样本数来支撑。样本数少于15时相关系数估计的置信区间很宽算出来的网络非常不稳前后两次抽样可能结果差异极大。我建议探索性分析至少15到20个样本。做严谨结论最好30个样本以上。如果样本量很少比如只有9个可以试着用limma或者variance stabilization之类的办法做表达量标准化但依旧要警惕结果的可靠性。如果样本实在不够还有个替代策略如果你有足够多的基因比如全转录组测序先用方差筛选出变化较明显的基因把维度降下来再跑WGCNA。但要注意这个做法改变了问题的定义你不再探索全基因组的共表达结构而是只针对高变基因找模块解读时要自我保护。4.2 千万别用“先筛差异基因再跑WGCNA”的思路这是新手最容易犯的错误我当年也在这上面栽过跟头。如果先做差异表达筛选再只挑差异基因跑WGCNA会有一个数学和逻辑上的双重问题逻辑上你丢掉了大量“表达量变化不大但在网络中发挥连接者作用”的基因模块结构会被严重扭曲。数学上差异基因往往是表达变化幅度大的基因它们之间的相关性容易虚高导致模块异常紧凑缺少层次。正确做法是把经过质量控制和标准化后的全基因组表达矩阵送入WGCNA让算法自己决定模块划分再回头去看哪些模块里的基因与性状相关。先定模块再看表型而不是先看表型筛基因再拿剩余基因定模块。4.3 输入数据的分布问题——别忽视标准化这一步RNA-seq的read counts如果不做任何标准化就拿去算相关性结果会受测序深度影响极大。测序深度高的样本整体表达量都高样本之间的假阳性相关性就会很明显很容易产生一个包含上千个基因的大模块纯粹就是技术噪声。常见的处理方法是用DESeq2的variance stabilizing transformation或rlog转换。或更直接的用log2(counts 1)只保留在较多样本中表达量达到一定阈值的基因。我自己偏好用variance stabilizing transformation它对低表达基因的方差压缩效果比log2(counts1)更平滑模块结构也更清晰。这个选择在跑WGCNA之前就应该定下来而不是跑完发现模块糊了再回头调。4.4 模块太多或者太少的时候怎么办模块数量是WGCNA分析中特别敏感的输出之一。模块太少比如五六个每个模块里的基因过于杂乱解释起来没有边界模块太多几十上百个大部分模块基因数太少不够稳健。处理策略如下按优先级排序先调整deepSplit参数1到3之间试影响动态剪切树的细碎程度。再调整minModuleSize影响小模块是否保留。实在不行再调整mergeCutHeight也就是模块合并阈值把过于相似的模块合并。我一般先固定minModuleSize 30和deepSplit 2跑一轮看模块数量和大小分布如果模块太少就改deepSplit 3或把minModuleSize降到20如果模块太多就反过来调。千万别一次性把所有参数都调掉那样你根本不知道哪个参数把结果带偏了。4.5 模块稳定性验证模块稳定是从结果可信度角度说必不可少的环节。WGCNA官方包不直接提供置换检验或分组重抽样的代码但实际项目里可以用一个简单的办法做自助抽样随机保留80%的样本重新跑一遍全流程看模块划分是否和全样本一致。如果你关心的关键模块在多次自助抽样里都能稳定出现结论的可信度就高很多。有些时候你会在不同批次数据里发现同一个模块始终存在但具体成员基因有一些变化这种时候我建议用模块特征基因做跨数据集的匹配而不是逐基因对比。5. 从原理到项目的落桥——一个典型的WGCNA实战流程框架5.1 全流程步骤梳理到了这一步把前面讲的所有原理串起来一个标准WGCNA项目大致长这样表达矩阵清洗过滤低表达基因去除离群样本检查样本聚类情况。数据标准化选择合适的标准化方式确保不同测序深度或批次效应不会主导相关性。计算软阈值β用pickSoftThreshold结合无标度拟合R²和模块数量的稳定性选β。构建邻接矩阵和TO矩阵调用adjacency函数和TOMsimilarity函数。层次聚类把1-TO矩阵转成距离用hclust聚类。动态树切割识别模块用cutreeDynamic确定初始模块。合并相似模块按模块特征基因的相关性合并高度相似的模块。模块-性状关联计算每个模块特征基因与性状的相关性和p值。模块内部挖掘找hub基因、模块内功能富集、模块成员与性状的关联。结果验证与可视化模块热图、特征基因聚类图、性状关联热图、hub基因网络图。这个流程不是死板的有些步骤可以根据研究目的调整。比如你只是想看看某个已知通路基因在哪些模块里聚集前面的模块识别流程可以简化。5.2 输出结果怎么解读——从一个实际项目实例说起给你举个例子我以前处理过一组肿瘤样本的转录组数据样本分早、中、晚期一共36例。跑完WGCNA之后我得到一个跟临床分级高度相关的模块相关性0.82p值远小于0.001模块里有240个基因。我对这个模块做GO富集发现显著富集到细胞外基质组织、胶原纤维组织、细胞迁移和血管生成等条目与肿瘤恶性进展的经典特征高度吻合。接着我找到模块内连通性最高的几个hub基因它们与模块特征基因的相关度都在0.9以上并且其中有两个基因在文献里已经有报道说与肿瘤侵袭转移有关。整个过程就是从“一群没有边界的差异基因”收敛到“一个具体的生物学程序”再到“少数几个候选核心基因”层层落实逻辑紧凑。这个例子并不是说WGCNA能告诉你哪个基因是“机制”但它能很有力地告诉你在那个与表型密切相关的基因圈子里哪些基因处于结构性的核心位置。后续验证工作比如细胞功能实验、组织验证、动物模型就都有了明确的下手目标。6. 个人使用心得与这个系列后续的安排6.1 工具选型为什么不直接手写而用WGCNA包WGCNA是一个R包运行环境要求不算高普通的16G内存工作站就能跑常见规模的数据。唯一的痛点是TO矩阵的构建对内存和CPU压力大。基因数超过2万、样本数接近100的时候构建TO矩阵会比较慢这时候可以设置参数把计算过程并行化或者在服务器上跑。目前WGCNA包在R 4.x版本下依然维护良好。官方的教程文档做得挺细致建议每个刚接触的人先去把官方Tutorial完整跑一遍再上手自己的数据。网上也有很多中文教程但整体水平良莠不齐尤其是对原理的解释建议以官方文档为准。6.2 这部分内容之后还会继续写什么这篇文章只聊了简介和原理对应这个系列的题目1。后续我会继续写WGCNA的实操篇包括从表达矩阵到模块划分的完整R代码流程每条代码都有注释。软阈值选择的可视化怎么看什么情况该选哪个值。模块与性状关联分析的实操细节连续型变量和二分组变量各自怎么处理。hub基因筛选和网络可视化的几种常用方案。以及我踩过的各种坑包括反复出现的“模块全是噪声”“结果在不同参数下完全不一样”这类问题都是怎么定位和处理的。把原理彻底理解了后面的实操才不会变成“抄代码”。这个道理我觉得是所有生信分析都通用的工具可以迭代理解才是真正能迁移的东西。