做转录组数据分析的朋友对WGCNA这个名字应该都不陌生。这是一套基于R语言的加权基因共表达网络分析方法全称Weighted Gene Co-expression Network Analysis。简单说它把你手上成百上千个基因的表达量数据按照表达模式的相关性组织成网络再从中找到表达行为相似的基因模块进一步和样本的性状、临床信息、处理条件等关联起来最终锁定那些可能起关键作用的基因。这套方法在医学、农学、植物学等领域都应用得非常广泛尤其在寻找生物标志物、候选基因、关键调控因子的场景里基本属于必跑的分析流程。这篇文章是系列的第一篇我会从拿到表达矩阵之后的第一步开始把数据预处理、样本聚类检查和软阈值选择这三个最基础也最容易卡住的环节用R代码加逐段解读的方式完整过一遍。内容适合已经有R语言基础、能看懂数据框操作、但对WGCNA流程还不熟的朋友。如果你是完全零基础建议先把R的基本数据结构、tidyverse常用的几个函数跑熟再回来看这篇会更顺手。第二篇再讲模块识别、模块与性状关联、hub基因筛选这些下游分析。1. WGCNA到底在做什么核心思路与应用场景1.1 从相关矩阵到共表达网络的底层逻辑WGCNA的中心思想听起来不复杂如果两个基因在不同样本里的表达量变化趋势高度一致那它们很可能参与同一个生物学过程或被同一个调控机制控制。传统做法是算皮尔逊相关系数设置一个阈值相关系数超过阈值的基因对就算“有关系”这是硬阈值。问题在于阈值怎么定都有点拍脑袋而且信息丢失严重——相关系数0.79和0.81的差别可能本身没有生物学意义却因为一条线被划成了完全不同的结果。WGCNA改用了软阈值的策略不把相关关系硬切掉而是把相关系数做幂指数运算权重按照相关系数的大小连续变化。公式是连接强度 (a_{ij} |cor(i,j)|^\beta)这个(\beta)就是软阈值。这样做的好处是保留了基因间关系的强度信息网络结构更加稳定也更符合生物系统连续渐变的特性。后面代码里pickSoftThreshold就是在帮你自动挑选合适的(\beta)值。另外一个关键点是WGCNA关注的不是单个基因对而是模块module。模块是一组高度互联的基因集合你可以把它理解为转录组里的“功能单元”。分析的核心产出之一就是把几千个基因归并成几十个以内可解释的模块再用模块特征基因module eigengene简称ME来代表每个模块的表达模式跟表型数据做关联分析。这一步能大幅降低分析的维度把“几千个基因”变成“十来个模块”生物学解释的可行性一下就上来了。1.2 什么时候该用WGCNA适用场景与前置条件WGCNA不是万能工具它有自己的适用边界。最理想的应用场景是样本量在15到50个之间每个样本都有对应的表型数据比如疾病组和对照组、不同发育时期、不同处理浓度等你有全转录组或全基因组的表达谱数据。样本量太小相关性估计不稳定样本量太大计算时间会非常感人不过WGCNA包对大数据集做了分块处理几百个样本也能跑只是后面我会提到内存和时间的代价。还有一类场景特别适合WGCNA那就是你不满足于只做差异表达分析想知道哪些基因“协同变化”想从系统层面理解转录调控。差异表达分析回答的是“哪些基因变了”WGCNA回答的是“哪些基因一起变这些基因和什么性状有关”两者互补不冲突。实操中我一般两个都做先跑差异表达拿到候选基因列表再用WGCNA看这些基因在共表达网络里的位置和归属。前置条件里最容易被忽略的是数据质量。WGCNA对缺失值敏感对极端离群样本敏感对批次效应也很敏感。正式跑网络构建前花点时间做数据清洗和样本聚类检查后面会省下大量排查问题的时间。这就是接下来第二部分要详聊的内容。2. 环境准备与数据预处理从表达矩阵到干净数据2.1 R包安装与版本选择WGCNA包在CRAN上可以直接安装但声明一下从CRAN安装的版本足够日常使用。安装命令如下# 安装BiocManager如果还没有 if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) # 安装WGCNA BiocManager::install(WGCNA) # 加载 library(WGCNA)需要注意的一点是WGCNA对R版本有要求建议安装前先更新R到当前的最新稳定版而不是用系统自带的老版本。装完之后建议跑一下:# 开启多线程后面计算会快很多 allowWGCNAThreads()这条命令会调用机器上可用的CPU线程对后续的网络构建加速非常明显。如果你在Windows上遇到线程相关的报错可以先运行disableWGCNAThreads()确认问题再决定是否启用。实际分析中多线程带来的提速在大样本量时候特别明显我曾经用24线程跑过一个200个样本的分析耗时从单线程的40分钟缩短到了不到10分钟所以这个细节建议不要忽略。2.2 输入数据格式与标准化的正确姿势WGCNA标准输入是一个表达矩阵行是基因列是样本值是表达量。至于这个表达量是RPKM、FPKM还是TPM理论上有争议但实操中大部分人直接采用经过标准化的表达值。如果做的是芯片数据用RMA或MASS归一化后的结果如果是RNA-seq建议用DESeq2的varianceStabilizingTransformation或EdgeR的CPM/logCPM转换把数据拉到一个近似正态分布的尺度上。数据导入之后的第一个正经步骤是过滤掉低表达基因。低表达基因的技术噪音大相关性不可靠留着只会干扰网络构建。常见的过滤策略有两种一是保留在所有样本中表达量大于某个阈值的基因二是保留方差排前一定比例的基因。下面是我常用的代码# 读取表达矩阵行为基因列为样本 expr - read.csv(expression_matrix.csv, row.names 1, check.names FALSE) # 过滤低表达基因至少在80%样本中表达量大于1 keep - rowSums(expr 1) 0.8 * ncol(expr) expr_filtered - expr[keep, ] # 过滤低变异基因保留MAD绝对中位差前75%的基因 mads - apply(expr_filtered, 1, mad) expr_filtered - expr_filtered[order(mads, decreasing TRUE)[1:ceiling(nrow(expr_filtered) * 0.75)], ] dim(expr_filtered)最后一步用mad排序再截取前75%是为了在保留生物学信号的同时尽量去掉噪音基因。你也可以保留更多比如前5000个高变基因这在芯片时代很常见。我个人更推荐从“所有表达达标的基因”出发而不是预先砍到几千个因为WGCNA本身就能通过模块化把基因分组你把输入砍得太狠反而可能丢掉与表型相关的基因。对于基因注释信息建议在表达矩阵里保留一列基因Symbol或Entrez ID后面做富集分析和数据库注释要用。分析过程中可能会遇到基因名重复的问题记得先去重否则相关性计算会出错# 假设第二列是基因名 merged - aggregate(expr_filtered[, -1], by list(gene expr_filtered$gene), FUN mean) rownames(merged) - merged$gene merged - merged[, -1]这种按基因名求平均的去重做法比简单保留第一个出现的基因更稳因为它把相同基因名下的多条探针或转录本表达量合并了。2.3 样本聚类与离群样本处理别急着建网络样本质量检查这一步很多人会跳过但我强烈建议不要跳过。样本聚类能直观地显示出是否有离群样本尤其是技术重复、批次差异导致的异常样本。如果再配合表型信息你还能看出聚类结果和分组是否一致。不一致才是正常的转录组数据的样本聚类本身就能反映样本的生物学特征。但如果同一个样本明显偏离所有其他样本通常意味着这个样本的数据质量有问题或者它收集时混入了一些非目标组织类型的细胞。操作方法很直接对样本做层级聚类画个树状图看看# 转置矩阵让行变成样本列变成基因 sampleTree - hclust(dist(expr_filtered), method average) pdf(sample_clustering.pdf, width 12, height 6) plot(sampleTree, main Sample clustering to detect outliers, sub , xlab , cex.lab 1.5, cex.axis 1.5) abline(h 阈值, col red) dev.off()这里dist()默认计算欧氏距离把每个样本的基因表达向量当成一个高维空间里的点计算点与点之间的距离。两个样本距离越近表示它们的基因表达模式越相似。method average表示聚类时类间距离取平均距离这种连接方式相对稳定是WGCNA官方教程里推荐的选择。聚类图里如果出现一个样本单独挂在一根很长的分支下面距离其他样本都很远就需要考虑是否剔除。判断标准一般看距离阈值比如把聚类高度大于某个值的样本视为离群。我在实际项目中遇到过不止一次因为某个样本的RNA质量差导致全基因组表达量异常偏低的情况这种样本不做剔除后续模块检测结果会被明显带偏。判断的时候要结合样本的QC指标综合判断不要仅凭聚类图就做决定。如果剔除样本后分组样本量不够可以考虑补充样本或者谨慎保留但一定要在文章里如实说明样本筛选过程。剔除离群样本后的矩阵建议重新保存一份用于后续分析# 假设sampleTree识别出的离群样本是 SAMPLE_017 keep_samples - !colnames(expr_filtered) %in% SAMPLE_017 expr_clean - expr_filtered[, keep_samples]3. 软阈值选择网络构建的关键一步3.1 无标度网络与软阈值的来龙去脉进入网络构建之前先解决一个绕不开的概念问题为什么WGCNA要用软阈值。WGCNA假设基因调控网络具有无标度网络的特征。无标度网络的特点是少数节点拥有大量连接hub大多数节点只有少数连接生物学网络比如蛋白质互作网络、代谢网络都近似符合这种特性。如果用硬阈值来构建网络相关性刚好超过阈值的基因对全部保留、刚好低于阈值的全部丢弃网络结构容易变得碎片化不稳定。更麻烦的是硬阈值选择的“拍脑袋”属性太强阈值从0.8改成0.85网络结构可能面目全非。软阈值的作用是给每对基因的连接强度做一个幂函数变换在保留网络连续性的同时让网络结构尽量逼近无标度特性。你的任务是找一个合适的(\beta)值使得加权后的网络连接度分布尽可能拟合无标度分布。WGCNA包里的pickSoftThreshold函数专门做这件事。它会计算一系列候选(\beta)值对应的拟合指数scale-free fit index就是R方同时统计网络的平均连接度然后帮你画一张图你根据这张图来选定最终的(\beta)。通常标准是选择第一个R方达到0.85以上的(\beta)值。但注意这个0.85不是绝对的。有些数据集比如来自复杂组织的转录组数据天然就很难达到高R方这时候综合考虑平均连接度的下降趋势选一个折中的(\beta)值就能继续分析。另外(\beta)也不需要设得过高(\beta)太高意味着对强相关的依赖度更大网络会变得稀疏且过于集中下游模块检测的稳定性反而受影响。3.2 用pickSoftThreshold挑选参数并解读输出下面这段代码是软阈值选择的模板# 允许多线程加速 allowWGCNAThreads() # 候选软阈值向量 powers - c(1:10, seq(from 12, to 30, by 2)) # 计算软阈值 sft - pickSoftThreshold( datExpr t(expr_clean), # WGCNA要求样本为行、基因为列 powerVector powers, verbose 5 ) # 查看软阈值评估结果 print(sft$fitIndices)这里有一个非常容易踩坑的地方datExpr参数需要的是行是样本、列是基因的矩阵和正常的表达矩阵行是基因、列是样本正好相反。所以传参时一定要转置。很多人第一步就在这里报错要么提示“Error: datExpr must be a data frame or matrix”要么后续结果乱套原因就是忘了转置。pickSoftThreshold的输出里有两列最关键SFT.R.sq和mean.k.。前者是候选阈值下网络拟合无标度分布的R方后者是平均连接度。画图可以用WGCNA自带的绘图函数# 设置画布分两栏 pdf(soft_threshold.pdf, width 12, height 6) par(mfrow c(1, 2)) # 左图无标度拟合指数 plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], xlab Soft Threshold (power), ylab Scale Free Topology Model Fit (Signed R^2), type n, main Scale independence) text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], labels powers, col red) abline(h 0.85, col red) # 右图平均连接度 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) dev.off()这段绘图代码里有个细节-sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2]是在处理无标度拟合的判定系数。WGCNA包输出里第三列是斜率斜率为负代表网络满足无标度特性即高连接度节点数量少所以取负号把R方修正为正值。如果某一行的斜率是正的说明网络结构发散不满足无标度假设这个候选阈值就不适合。实际判断的时候我会先看左图找R方首次超过0.85的(\beta)然后看右图确认平均连接度没有跌到特别低。比如之前分析一个植物转录组数据R方在(\beta6)时达到了0.91但(\beta6)时平均连接度还有70多网络不会太稀疏我就选了6。另一个动物组织数据R方在(\beta12)时是0.87但平均连接度只有6模块后续非常碎这时候我会下调到(\beta10)来平衡保证模块的规模能用生物学语言解释。4. 网络构建与模块识别实操4.1 一步法构建网络与参数含义选定软阈值之后就可以构建网络并识别模块了。有两种路径一种是完全手动分布执行先算邻接矩阵、再算TOM矩阵、再做层级聚类、最后动态剪枝另一种是直接用blockwiseModules函数一步完成。日常分析推荐直接使用后者它把全套流程封装好了参数也足够灵活。# 网络构建与模块识别 net - blockwiseModules( datExpr t(expr_clean), power 6, # 上一步选定的软阈值 TOMType unsigned, # TOM的计算类型 minModuleSize 30, # 最小模块基因数 reassignThreshold 0, # 重分配阈值 mergeCutHeight 0.25, # 模块合并的聚类高度阈值 numericLabels TRUE, # 模块用数字标识 pamRespectsDendro FALSE, saveTOMs TRUE, saveTOMFileBase TOM_data, # TOM矩阵保存的文件名 verbose 3 )这段代码里最容易让人犹豫的是TOMType参数。它有三个取值unsigned、signed和signed hybrid。unsigned只考虑相关性的绝对值基因正相关和负相关都会算作连接signed只把正相关算作连接负相关不计入。转录组分析中共调控的基因通常是正相关关系signed更符合生物学预期但负相关的基因也有可能属于同一个调控通路上的抑制关系。我个人做转录组默认用signed如果做的是芯片数据且表达谱分布比较宽可以考虑unsigned。关键是前后一致不要一会儿signed一会儿unsigned导致结果对不上。mergeCutHeight是模块合并的阈值。因为直接聚类出来的模块可能非常碎把聚类树上距离接近的模块合并掉可以减少下游分析中模块太多、过度碎片化的问题。值越小合并越少常见范围0.15到0.3。我一般先用默认0.25跑一遍看模块数量和合并结果如果模块数量太多、每个模块基因太少就调大一点如果模块数太少导致很多表达模式差异明显的模块被合并吞掉就调小一点。这一点需要反复尝试结合你后续模块与性状关联的结果来最终确定。4.2 识别模块结果与模块特征基因运行结束后net对象里包含多个结果组件。最常用的是net$colors它给每个基因打了一个模块编号的标签。数字0表示没有被分配到任何模块的基因这些基因会被后面的分析忽略。可以用代码快速看下模块的基因数分布# 查看模块颜色与基因数 moduleColors - labels2colors(net$colors) table(moduleColors)labels2colors函数会把数字标签转成颜色名方便画图时使用。模块基因数如果出现很多个低于30的小模块说明minModuleSize设小了或者mergeCutHeight设小了。对于下游分析模块大小在50到2000之间比较理想。太小的模块生物学重复性差太大的模块内部异质性高难以解读。模块识别完之后WGCNA会为每个模块计算模块特征基因module eigengene用来代表这个模块的整体表达模式# 计算模块特征基因 MEs - moduleEigengenes(t(expr_clean), moduleColors)$eigengenes # 查看模块特征基因的样本×模块矩阵 head(MEs)特征基因本质上是对模块内所有基因的表达矩阵做PCA取第一主成分。它用一组数值概括了每个样本在该模块上的整体表达水平。有了这个值后面跟表型做相关性分析就变得非常直接。一个模块的特征基因如果与某个细胞类型、临床指标高度相关那这个模块很可能反映了该表型相关的生物学过程。查看模块特征基因之间的相关性以及样本间的模块表达模式通常用moduleEigengenes之后接着做heatmap或者plotEigengeneNetworks。这一步能帮你快速发现哪些模块之间表达模式相似可能共同参与某个通路。我在实践中经常发现两个模块的特征基因相关系数超过0.8这就是你后续做模块合并或者解读功能时要特别注意的信号。5. 常见报错与排查技巧实操中踩过的坑5.1 数据格式与数据处理类报错从一开始就把格式处理好能避免很多不必要的折腾。但我还是把实际中遇到过的典型报错整理出来方便你排查对照。报错一Error: datExpr must be a data frame or matrix几乎都是因为没转置矩阵。WGCNA整个流程里datExpr要求的是行是样本、列是基因。如果你读入的表格是行基因列表记得在传给pickSoftThreshold或blockwiseModules前运行t()转置。还有个容易忽略的点t()之后原本的data.frame会变成matrix有些后续操作需要再转回data.frame否则某些函数会报类型错误。报错二Error: y must be numeric或聚类时报相似性矩阵找不到这通常是因为表达矩阵里有非数值列比如基因注释列没有去掉。dist()函数只能对数值矩阵计算距离如果你的矩阵里混入了字符型数据就会报错。解决方法是在读取数据后把基因注释列拆出来单独存表达量部分的矩阵确保每个列都是numeric类型。报错三运行blockwiseModules时内存不足或卡死当基因数量超过2万、样本数量超过100时TOM矩阵的计算量迅速增长内存消耗会非常惊人。一种方案是设置maxBlockSize参数让WGCNA自动分块计算。另一种是在服务器上跑增加内存。如果只能在本机跑我建议把输入基因先做个过滤比如只保留MAD前10000的基因能大幅降低计算量。实际分析中影响不大但内存问题能很快缓解。5.2 软阈值选择与模块识别的思路调整软阈值选定后模块识别结果不理想的情况也经常出现。比较常见的情形有两种一是没有达到0.85的R方二是模块太多太碎、难以解读。这两种情况处理思路不太一样。先看第一种没有达到0.85。这时需要回到软阈值图找一个拐点也就是R方上升变缓的位置。比如候选阈值从6到8R方从0.82涨到0.86但平均连接度从80跌到30我一般会选6而不是8。因为更高的(\beta)让网络过于稀疏模块会变得碎片化后续结果稳定性差。0.85是参考值而不是绝对标准判断时看到的是网络拓扑特性整体是否合理。再看第二种模块太多太碎。最常见的原因是mergeCutHeight太小模块合并不够。我的习惯是先输出不同mergeCutHeight值下的模块数量画个对照表再从中挑一个模块数量适中、模块大小分布合理的值。另一个原因是minModuleSize设太小如果设置15很容易产生大量几十个基因的小模块。建议至少设到30如果模块数量还是多可以试试50。当然如果你的研究目标就是找特定的小模块那可以适当调低阈值但要做好结果稳定性检查。操作中还有一个小细节值得注意WGCNA的运行结果和随机种子有关。blockwiseModules内部用到的一些随机过程不是完全确定的如果你的分析文章需要可重复性建议在代码开头设置随机种子set.seed(2024)这样别人按同一份代码跑结果能完全复现。不同版本R或不同平台下分块方式不同会导致结果细微差异但设置种子能保证绝大多数情况下的一致性写文章和做报告时这个细节会显得你很专业。实操总结与两点补充心得到这里WGCNA分析流程的第一阶段就跑通了。梳理一下第一篇的内容数据清洗与低表达基因过滤、样本聚类与离群样本剔除、软阈值选择、一步法网络构建与模块识别、模块特征基因计算。跑完这套流程你已经能看到基因模块的概貌了。第二篇会继续往前走做模块与性状的关联分析、模块内基因的功能富集解读、hub基因筛选以及最后那些漂亮的网络可视化图是怎么画出来的。最后分享两个实际操作中积累的小经验。第一个是关于数据分析记录的WGCNA分析涉及大量参数选择和反复尝试建议每次调参都记录对应的结果比如模块数量、模块大小分布等方便事后追溯和复现也方便写文章时说明参数选择的依据。第二个是关于结果的解释WGCNA跑出来只是第一步模块里富集到哪些通路、与哪些性状显著相关才决定模块能否转化为一个有故事可讲的生物学发现。这两点在第二篇里都会有更详细的实际操作演示到时候结合代码一起看会更有感觉。