拿到转录组表达矩阵以后很多人第一反应是跑差异表达筛出一堆显著基因然后再去补一个GO/KEGG富集分析。但真到解释生物学意义的时候往往会卡住差异基因太多太散彼此之间到底怎么协同工作哪些基因模块跟临床性状、表型数据最相关哪些基因是模块里的“核心枢纽”这些问题如果不借助共表达网络很难回答得让人信服。WGCNAWeighted Gene Co-expression Network Analysis加权基因共表达网络分析解决的就是这个环节。这是一套在R语言里完成的系统性分析流程从表达矩阵清洗、样本质检到软阈值选择、网络构建、模块识别再到模块和性状的关联分析一步扣一步。作为“全步骤”系列的第一篇我会把从零开始跑通WGCNA上游核心流程的完整思路写清楚所有代码都附上逐行解读不只是让你能复制粘贴更要让你知道每一行在干什么、为什么这么干。如果你正拿着转录组表达数据不知道下一步怎么挖或者刚接触WGCNA看着官方教程一头雾水这篇应该能帮你省掉不少折腾的时间。1. 先搞清楚WGCNA到底在解决什么问题1.1 共表达网络的基本逻辑基因不是孤立工作的我们可以把一个样本里成千上万个基因的表达量想象成一群人的工作状态。某些基因的表达量在所有样本里总是同涨同跌说明它们大概率在同一个生物学通路里协同工作或者受同一个上游转录因子调控。WGCNA把这种“同步变化”的关系抽象成网络每个基因是一个节点两个基因之间有没有边、边的权重有多大取决于它们的表达量在所有样本里的相关性强不强。传统做法是直接计算两两基因间的Pearson相关系数然后设一个阈值比如相关系数大于0.8就认为有连接。这样做的最大问题是阈值是拍脑袋定的而且硬性切割会把很多弱但真实存在的共表达关系丢掉。WGCNA的不同之处在于它不设硬阈值而是用一个软阈值soft-thresholding power对相关系数做幂指数加权让网络尽量符合无标度拓扑特征——简单说就是让网络里少数基因拥有大量连接大多数基因只有少量连接这是很多真实生物网络具备的特征。1.2 WGCNA和普通相关分析的核心差别普通相关分析是在“基因对”层面做检验几万个基因会产生上亿个基因对多重检验校正会让结果变得非常保守而且你很难从海量相关关系里看出整体结构。WGCNA的处理方式是先构建全基因网络再用层次聚类的方法把网络划分成若干个模块module每个模块里的基因表达模式高度相似。这样分析单元就从“单个基因”升级成了“基因模块”后续做模块与性状关联、筛选枢纽基因、富集分析逻辑上都更清晰。1.3 一篇WGCNA能回答哪些问题在实际项目中WGCNA最常见的用途有三个找到与某个性状比如疾病/正常、药物响应、发育阶段、产量高低显著相关的基因模块。在每个模块内部筛选高连接度的枢纽基因hub gene作为后续实验验证的候选靶点。把模块基因拿去做功能富集解释模块代表的生物学通路。这套流程不是万能的它对样本量有一定要求对输入数据质量也很敏感但这些细节放在后面实操部分再讲。先把逻辑框架搭好后面每一步代码才不会变成“黑盒操作”。2. 数据准备两个文件的格式和你最容易踩的坑2.1 表达矩阵的格式行是样本列是基因WGCNA的输入数据格式非常明确表达矩阵必须是“行样本列基因”的二维数据框。我们经常拿到的原始表达矩阵是“行基因列样本”所以读入后通常要做一次转置。这也是新手第一处容易出错的地方很多报错都源自行列搞反了。library(WGCNA) # 关闭字符串自动转因子避免后面矩阵运算出幺蛾子 options(stringsAsFactors FALSE) # 读入表达矩阵行名是基因ID列名是样本名 datExpr_raw - read.csv(expression_data.csv, row.names 1, check.names FALSE) # 转置成 行样本列基因 datExpr0 - as.data.frame(t(datExpr_raw))这里有个细节check.names FALSE是为了防止R把样本名里的横杠、空格之类改写成点号否则后面和性状文件匹配样本名时经常会莫名其妙对不上。2.2 表达量数据到底该用什么很多人问WGCNA能不能直接输入DESeq2得到的标准化counts理论上可以但我更推荐用FPKM、TPM这类已经校正过基因长度和测序深度的表达量而且要注意数据的分布形态。标准的RNA-seq表达矩阵数值跨度很大直接拿来做相关性计算会被高表达基因主导一般都需要做log2(x 1)变换让数据更接近正态分布。# 如果数据还没有log变换做一个log2(x1)变换 datExpr0 - log2(datExpr0 1)如果你是拿芯片数据或者qPCR数据来跑也要先确认数值范围合理。表达量差距在几个数量级、又没有做变换的数据跑出来模块结构通常会非常碎而且很难解释。2.3 基因过滤挑表达量稳定且非零的基因WGCNA官网教程里的示例数据是已经预处理好的但我们自己处理的数据里经常有大量低表达基因。这类基因的“表达变化”很多是测序噪声计算相关矩阵时会把网络搅乱还会显著拖慢计算速度。我一般会分两步过滤第一步用goodSamplesGenes检查缺失值和标准差为零的基因把它们剔除。# 检查样本和基因的基本质量 gsg - goodSamplesGenes(datExpr0, verbose 3) # 如果存在不合格的基因或样本直接过滤掉 if (!gsg$allOK) { datExpr0 - datExpr0[gsg$goodSamples, gsg$goodGenes] }第二步过滤表达量太低的基因以及所有样本里表达量基本不变的基因并用中位绝对偏差MAD筛出表达量变化最明显的那些基因。# 计算每个基因在所有样本中的中位数和MAD gene_median - apply(datExpr0, 2, median) gene_mad - apply(datExpr0, 2, mad) # 保留中位数表达量较高、且变异较大的基因 keep - gene_median 1 gene_mad 0.5 datExpr - datExpr0[, keep]是不是一定要做这一步如果基因数量在2万左右机器配置也足够网络构建其实也能跑完。但过滤之后模块会更稳定后续模块-性状关联的显著性也更容易出现。这里需要把握一个度如果过滤太狠可能会丢掉重要基因如果不过滤噪声又会影响结果。比较稳妥的做法是先做常规低表达过滤再用MAD筛掉后20%的基因保留约8000到15000个基因进行网络构建。2.4 样本聚类肉眼剔除异常样本这一步很多人会跳过但我觉得它比后面的参数调优更关键。样本聚类树可以直观地暴露问题如果一个样本和其他样本离得非常远说明它的整体表达模式异常可能是实验批次差异、样品污染或者数据预处理出了问题。# 样本聚类 sampleTree - hclust(dist(datExpr), method average) # 画图观察 pdf(sample_cluster.pdf, width 12, height 6) plot(sampleTree, main Sample clustering to detect outliers, sub , xlab , cex.lab 1.5, cex.axis 1.5, cex.main 2) abline(h 100, col red) dev.off()如果图里出现明显的离群样本直接用下标把它剔除# 假设样本聚类图中第3个样本明显离群 datExpr - datExpr[-3, ]cutHeight阈值不是固定不变的要看聚类树的高度分布来定。我的习惯是先用目测选一个把大多数样本聚在一起、只把极少数离群样本切出去的高度不用刻意追求统一标准。2.5 性状数据文件怎么整理性状数据是WGCNA做模块关联分析的核心输入。格式要求是“行样本列性状”样本名要和表达矩阵的行名完全一致。性状可以是连续变量比如年龄、血压、药物浓度也可以是分组变量但需要转成数值型比如疾病组1、对照组0或者多分组设计使用0/1哑变量编码。# 读取性状数据 datTraits - read.csv(sample_traits.csv, row.names 1, check.names FALSE) # 查看两个数据集的样本名交集 common_samples - intersect(rownames(datExpr), rownames(datTraits)) datExpr - datExpr[common_samples, ] datTraits - datTraits[common_samples, , drop FALSE]这里要特别注意必须保证表达矩阵和性状文件的样本一一对应。我一开始跑的时候就因为表达矩阵和性状文件样本顺序不一致直接栽过跟头。用intersect统一样本后后续分析就不会出现张冠李戴的问题。3. 软阈值选择别只会默认选9要学会看两张图3.1 无标度拓扑准则到底是什么WGCNA里最核心的一个概念是“软阈值”power。简单理解它就是一个加权系数把基因间的相关系数取绝对值的power次方得到基因之间的邻接权重。power越大弱相关被压制得越厉害网络的稀疏程度也越高。但power并不是越大越好WGCNA选择power的依据是让网络尽可能地符合无标度拓扑我们希望网络中存在少数连接度极高的“枢纽节点”而大部分节点连接度较低。衡量网络是否符合这个特征的指标是SFT.R.sq也就是无标度拓扑拟合指数。通常来说我们希望这个值尽量高尤其是要超过0.85。3.2 pickSoftThreshold的代码和输出选择软阈值不需要自己瞎试WGCNA包提供了pickSoftThreshold函数# 选择一系列候选power值 powers - c(1:10, seq(from 12, to 30, by 2)) # 计算不同power下的无标度拟合指数和平均连接度 sft - pickSoftThreshold(datExpr, powerVector powers, verbose 5) # 把结果整理成一个数据框查看 fit - sft$fitIndices print(fit[, c(Power, SFT.R.sq, mean.k., median.k.)])输出结果会看到一张类似下面的表格PowerSFT.R.sqmean.k.median.k.10.0812145.31942.120.2121023.7856.440.553366.8245.360.742178.496.280.853103.647.3100.91266.224.2120.94845.113.6眼睛不要只盯着哪个power大要同时看SFT.R.sq和mean.k.。判断标准是选择最小的、让SFT.R.sq首次进入0.85以上区间的power同时平均连接度不能太低否则网络会过度稀疏模块识别失去意义。3.3 组合图怎么看通常还会画一张双面板的组合图左边是power和拟合指数SFT.R.sq的折线右边是power和平均连接度的趋势线pdf(soft_threshold.pdf, width 10, height 5) par(mfrow c(1, 2)) cex1 - 0.9 # 左图SFT.R.sq plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], xlab Soft Threshold (power), ylab Scale Free Topology Model Fit (R^2), type n, main Scale independence) text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], labels powers, col red, cex cex1) abline(h 0.85, col red, lty 2) # 右图平均连接度 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, cex cex1) dev.off()这里有一个WGCNA代码里比较反直觉的地方左图纵轴为什么是-sign(x[,3]) * x[,2]因为当power比较低时无标度拟合指数可能是负斜率乘上负号能让所有点都朝上显示方便观察。看的时候只要看绝对值就可以了。3.4 如果R平方一直上不了0.85怎么办这是实际分析中最常见的问题之一。遇到这种情况先不要急着怀疑数据有问题。有几种处理思路检查是否没有过滤低表达基因噪声太大影响了网络结构先回去重新做数据清洗。尝试更广的power范围比如seq(1, 40, by 2)。对表达矩阵再做一次更严格的标准差筛选。样本量太小比如少于15个样本时SFT.R.sq可能很难达到0.85此时可以考虑退而求其次选择曲线转折点附近的power同时结合网络生物学可解释性来选定。我个人在临床小样本数据上遇到过几次R平方只能到0.7的情况。我的做法是选取一个中间偏上的power比如根据曲线趋势选10或者12然后看后续模块是否稳定、模块与性状的关联是否有意义。如果模块结果乱七八糟再回头调整power。4. 核心代码块网络构建与模块识别4.1 blockwiseModules函数参数逐个拆解选定power之后就到了整个WGCNA分析的重头戏——构建网络并识别模块。标准代码并不复杂核心就是blockwiseModules# 用选定的power构建网络 net - blockwiseModules( datExpr, power 9, TOMType unsigned, minModuleSize 30, reassignThreshold 0, mergeCutHeight 0.25, numericLabels TRUE, pamRespectsDendro FALSE, saveTOMs TRUE, saveTOMFileBase TOM-blockwise, verbose 3 )这个函数到底做了什么它可以拆成三层理解先计算基因间的相关性矩阵再换算成邻接矩阵然后计算TOM相似度拓扑重叠最后基于TOM相异度做层次聚类和动态剪枝得到模块。下面逐个参数说明每个参数我都会结合自己踩过的坑来讲。power就是软阈值。选9还是选12取决于第3节pickSoftThreshold的结果不要无脑用默认值。TOMType推荐用unsigned。它的意思是把负相关也当作共表达关系因为基因调控网络里存在大量负调控一个转录因子抑制下游基因时它们的表达模式是负相关的。如果改用signed则只保留正相关基因之间的共表达关系模块会更保守。minModuleSize最小模块基因数。默认值是30但是基因数少的平台数据可以降到20甚至10。设置太小会导致模块数量爆炸很多只有一两个基因的“迷你模块”没有生物学意义。mergeCutHeight相似模块合并的阈值。默认0.25意思是模块间特征基因相关性高于0.75就需要考虑合并。如果后续发现模块多而碎可以适当调低mergeCutHeight到0.2或者在后期手动做一次模块合并。numericLabels如果为TRUE模块会以数字编号命名输出的是0、1、2…… 如果为FALSE会用颜色命名比如turquoise、blue、brown。我个人建议代码阶段先用数字导出绘图时再映射成颜色这样在后续匹配基因时更方便。pamRespectsDendro是否在动态剪枝后继续用PAMPartitioning Around Medoids对模块边界做进一步细化。如果设为TRUE生成的模块更紧凑但有时会过度切割如果设为FALSE会更容易得到比较大的、稳健的模块。默认是FALSE我一般也保持FALSE。saveTOMs和saveTOMFileBase是否把TOM矩阵保存到本地。TOM矩阵是WGCNA里最大的中间产物会占很大磁盘空间但保存下来可以避免后续重复计算尤其是在做不同power对比时能节省大量时间。4.2 从net对象里能拿到什么blockwiseModules跑完之后net对象里包含了做后续分析所需的几乎所有关键结果# 查看模块划分结果 table(net$colors) # 查看模块特征基因 net$MEs # 查看模块数量 length(unique(net$colors))net$colors是一个和基因列表长度相同的向量记录了每个基因被分配到了哪个模块0表示未进入任何模块的灰色基因。net$MEs是每个样本在每个模块上的特征基因表达值后面做模块-性状关联时主要靠它。4.3 把聚类树画出来看看整体结构模块识别之后最直观的可视化方式是画聚类树和模块颜色条pdf(dendrogram_modules.pdf, width 12, height 8) plotDendroAndColors( net$dendrograms[[1]], moduleColors net$colors, groupLabels Module colors, main Gene dendrogram and module colors, dendroLabels FALSE, addGuide TRUE, hang 0.03 ) dev.off()这张图是整个WGCNA结果里我最先看的一张。理想情况下树状图会分成几个明显的大分支每个分支对应一个颜色模块。如果看到颜色条像斑马线一样反复横跳或者模块数量特别多就要考虑调整minModuleSize、mergeCutHeight或者回去检查数据预处理。刚接触WGCNA的时候我拿到模块结果第一反应是赶紧去看模块和性状的关联结果画出来的聚类树乱得没法看。后来才意识到聚类树的稳定程度本身就能反映数据质量如果这一层就有问题后面所有统计都会跟着出问题。4.4 基因数量与模块大小的合理范围一个合理的网络通常会有10到30个模块其中最大的模块基因数量占10%到30%同时会有几个中等模块和少量小模块。如果所有基因都堆在同一个巨大模块里说明power选得太低模块没有分开如果模块数量超过50个要么过滤不充分要么minModuleSize设得太小。这些标准不是绝对死线但能帮助快速判断结果是否靠谱。5. 模块合并与模块-性状关联分析5.1 为什么要做模块合并blockwiseModules里的mergeCutHeight参数会自动完成相似模块的合并。但不同版本代码或者手动调整后可能还需要自己再检查一遍。相似模块指的是两个模块的特征基因MEModule Eigengene高度相关它们本质上是同一群基因的不同变体合并后更容易解释也减少后续多重检验的次数。如果需要手动合并代码是这样的# 计算模块特征基因 MEs - moduleEigengenes(datExpr, net$colors)$eigengenes # 对模块进行合并 merge - mergeCloseModules(datExpr, net$colors, cutHeight 0.25, verbose 3) # 合并后的模块颜色 mergedColors - merge$colors # 合并后的模块特征基因 mergedMEs - merge$newMEs合并后一定要重新把模块颜色画一遍确认不要直接跳过可视化。5.2 模块-性状关联的计算方式每个模块的ME是这个模块所有基因表达模式的第一主成分代表整个模块在样本间的主要变化趋势。把ME和每个性状做相关分析就能得到模块与性状的关联矩阵。# 计算模块特征基因与性状的相关性 nSamples - nrow(datExpr) moduleTraitCor - cor(mergedMEs, datTraits, use p) moduleTraitPvalue - corPvalueStudent(moduleTraitCor, nSamples)corPvalueStudent是WGCNA为这个场景专门封装好的函数它会基于相关系数和样本量给出p值比手算方便很多。模块与性状关联的结果通常画成热图横轴是性状纵轴是模块每个格子里的数字是相关系数和括号里的p值pdf(module_trait_heatmap.pdf, width 8, height 10) textMatrix - paste(signif(moduleTraitCor, 2), \n(, signif(moduleTraitPvalue, 1), ), sep ) dim(textMatrix) - dim(moduleTraitCor) labeledHeatmap( Matrix moduleTraitCor, xLabels colnames(datTraits), yLabels colnames(mergedMEs), ySymbols colnames(mergedMEs), colorLabels FALSE, colors blueWhiteRed(50), textMatrix textMatrix, setStdMargins FALSE, cex.text 0.6, main Module-trait relationships ) dev.off()重点看两个指标相关系数的绝对值大小以及p值的显著性水平。比如某个模块与“疾病状态”的相关系数是0.72p值小于0.001那这个模块就是你后续要重点挖掘的候选模块。5.3 从模块里挑基因基因显著性与模块成员度找到一个与性状显著相关的模块后还需要知道模块里哪些基因起主导作用。WGCNA提供了两个经典指标基因显著性GS, Gene Significance基因表达量与性状之间的相关绝对值表示该基因与性状的关联强度。模块成员度MM, Module Membership基因与该模块特征基因的相关性表示该基因在模块内的核心程度。# 选择一个性状比如第1列 trait_column - 1 # 计算每个基因与性状的相关性 geneTraitSignificance - as.data.frame(cor(datExpr, datTraits[, trait_column], use p)) GS - as.numeric(geneTraitSignificance$V1) # 计算每个基因与模块特征基因的相关性 MM - as.data.frame(cor(datExpr, mergedMEs, use p))然后可以拿GS和MM做散点图筛选兼具模块核心地位且和性状显著相关的基因。比如模块内MM 0.8且GS 0.2的基因往往就是值得后续验证的候选枢纽基因。5.4 灰色模块怎么处理grey模块是WGCNA里一个特殊的存在它是所有没能被划分到任何模块的基因集合。灰色模块通常不会和性状显著相关如果有也一样要关注一下可能意味着数据里还存在另一种独立的表达模式值得单独做一次亚聚类分析。多数情况下灰色模块基因不参与后续分析。6. 实操中的零散问题与排查建议6.1 样本量太小怎么办WGCNA对样本量的最低要求保守说至少要有15到20个样本如果少于这个数基因相关性估计会非常不稳定模块结果也很难重复。真遇到小样本数据我有几个折中经验把minModuleSize调低到10左右power选择稍微偏大一点让网络更稀疏同时使用signed TOM而不是unsigned往往会稳一些。但也要明确小样本条件下的WGCNA结果只能作为探索性分析不适合直接下强结论。6.2 计算太慢、内存爆掉基因数量接近两万时blockwiseModules默认会做分块计算因为一次性计算全部基因的TOM矩阵对内存压力很大。如果仍然卡死建议先确认R是64位版本然后适当降低保留基因数比如从15000降到10000。还可以开启多线程# 启用多线程靠CPU核心数决定 enableWGCNAThreads()注意这条命令要在加载WGCNA后、构建网络之前运行而且Windows系统下多线程支持不如Linux/macOS稳定。6.3 结果可重复性问题WGCNA里有一些步骤依赖随机性尤其是样本量不特别大的时候模块识别结果可能会有轻微波动。想保证结果可重复建议在脚本开头设置随机种子set.seed(2024)另外分析过程中所有中间结果都要及时保存save(datExpr, datTraits, net, mergedColors, mergedMEs, file wgcna_step1.rda)这样即使后面改参数也不需要从最开始的读取文件重新跑一遍。6.4 包版本差异WGCNA包这些年更新不算频繁但不同小版本的默认参数可能略有差异比如blockwiseModules里的某些参数在新版本提示deprecated。如果碰到函数调用报错先看包自带的NEWS文档再看官方教程。网上很多老教程用的代码在新版本下可能会跑不通这不是你写错了需要根据报错信息微调函数名或参数。最后再分享一个偷懒技巧整套流程里最耗时的往往是不同power下的网络构建对比。我会在跑正式分析前先拿一小部分基因比如随机抽3000个基因做一个快速测试看看模块数量是否合理、聚类树是否稳定。参数基本满意后再用全部基因跑正式版本。这个小技巧能省下大量反复调试的时间。另外blockwiseModules生成的TOM文件占空间非常大分析做完如果没有特殊需要记得及时清理不然一个项目下来几百GB一点也不夸张。把模块结果、基因颜色、特征基因这些核心结果保留好就足够了。这一篇把WGCNA从数据准备到模块-性状关联的完整前半程梳理完了代码基本可以直接照着改路径和数据跑通。下一篇我会继续写模块内部的可视化细节包括基因网络导出到Cytoscape、hub gene筛选、以及怎么把模块基因批量提交给富集分析工具。跑代码的过程中遇到具体的报错和诡异结果欢迎照着这篇的排查思路先自己试一圈多数问题都出在数据格式和样本匹配上。
WGCNA共表达网络分析全流程：从表达矩阵到基因模块挖掘
拿到转录组表达矩阵以后很多人第一反应是跑差异表达筛出一堆显著基因然后再去补一个GO/KEGG富集分析。但真到解释生物学意义的时候往往会卡住差异基因太多太散彼此之间到底怎么协同工作哪些基因模块跟临床性状、表型数据最相关哪些基因是模块里的“核心枢纽”这些问题如果不借助共表达网络很难回答得让人信服。WGCNAWeighted Gene Co-expression Network Analysis加权基因共表达网络分析解决的就是这个环节。这是一套在R语言里完成的系统性分析流程从表达矩阵清洗、样本质检到软阈值选择、网络构建、模块识别再到模块和性状的关联分析一步扣一步。作为“全步骤”系列的第一篇我会把从零开始跑通WGCNA上游核心流程的完整思路写清楚所有代码都附上逐行解读不只是让你能复制粘贴更要让你知道每一行在干什么、为什么这么干。如果你正拿着转录组表达数据不知道下一步怎么挖或者刚接触WGCNA看着官方教程一头雾水这篇应该能帮你省掉不少折腾的时间。1. 先搞清楚WGCNA到底在解决什么问题1.1 共表达网络的基本逻辑基因不是孤立工作的我们可以把一个样本里成千上万个基因的表达量想象成一群人的工作状态。某些基因的表达量在所有样本里总是同涨同跌说明它们大概率在同一个生物学通路里协同工作或者受同一个上游转录因子调控。WGCNA把这种“同步变化”的关系抽象成网络每个基因是一个节点两个基因之间有没有边、边的权重有多大取决于它们的表达量在所有样本里的相关性强不强。传统做法是直接计算两两基因间的Pearson相关系数然后设一个阈值比如相关系数大于0.8就认为有连接。这样做的最大问题是阈值是拍脑袋定的而且硬性切割会把很多弱但真实存在的共表达关系丢掉。WGCNA的不同之处在于它不设硬阈值而是用一个软阈值soft-thresholding power对相关系数做幂指数加权让网络尽量符合无标度拓扑特征——简单说就是让网络里少数基因拥有大量连接大多数基因只有少量连接这是很多真实生物网络具备的特征。1.2 WGCNA和普通相关分析的核心差别普通相关分析是在“基因对”层面做检验几万个基因会产生上亿个基因对多重检验校正会让结果变得非常保守而且你很难从海量相关关系里看出整体结构。WGCNA的处理方式是先构建全基因网络再用层次聚类的方法把网络划分成若干个模块module每个模块里的基因表达模式高度相似。这样分析单元就从“单个基因”升级成了“基因模块”后续做模块与性状关联、筛选枢纽基因、富集分析逻辑上都更清晰。1.3 一篇WGCNA能回答哪些问题在实际项目中WGCNA最常见的用途有三个找到与某个性状比如疾病/正常、药物响应、发育阶段、产量高低显著相关的基因模块。在每个模块内部筛选高连接度的枢纽基因hub gene作为后续实验验证的候选靶点。把模块基因拿去做功能富集解释模块代表的生物学通路。这套流程不是万能的它对样本量有一定要求对输入数据质量也很敏感但这些细节放在后面实操部分再讲。先把逻辑框架搭好后面每一步代码才不会变成“黑盒操作”。2. 数据准备两个文件的格式和你最容易踩的坑2.1 表达矩阵的格式行是样本列是基因WGCNA的输入数据格式非常明确表达矩阵必须是“行样本列基因”的二维数据框。我们经常拿到的原始表达矩阵是“行基因列样本”所以读入后通常要做一次转置。这也是新手第一处容易出错的地方很多报错都源自行列搞反了。library(WGCNA) # 关闭字符串自动转因子避免后面矩阵运算出幺蛾子 options(stringsAsFactors FALSE) # 读入表达矩阵行名是基因ID列名是样本名 datExpr_raw - read.csv(expression_data.csv, row.names 1, check.names FALSE) # 转置成 行样本列基因 datExpr0 - as.data.frame(t(datExpr_raw))这里有个细节check.names FALSE是为了防止R把样本名里的横杠、空格之类改写成点号否则后面和性状文件匹配样本名时经常会莫名其妙对不上。2.2 表达量数据到底该用什么很多人问WGCNA能不能直接输入DESeq2得到的标准化counts理论上可以但我更推荐用FPKM、TPM这类已经校正过基因长度和测序深度的表达量而且要注意数据的分布形态。标准的RNA-seq表达矩阵数值跨度很大直接拿来做相关性计算会被高表达基因主导一般都需要做log2(x 1)变换让数据更接近正态分布。# 如果数据还没有log变换做一个log2(x1)变换 datExpr0 - log2(datExpr0 1)如果你是拿芯片数据或者qPCR数据来跑也要先确认数值范围合理。表达量差距在几个数量级、又没有做变换的数据跑出来模块结构通常会非常碎而且很难解释。2.3 基因过滤挑表达量稳定且非零的基因WGCNA官网教程里的示例数据是已经预处理好的但我们自己处理的数据里经常有大量低表达基因。这类基因的“表达变化”很多是测序噪声计算相关矩阵时会把网络搅乱还会显著拖慢计算速度。我一般会分两步过滤第一步用goodSamplesGenes检查缺失值和标准差为零的基因把它们剔除。# 检查样本和基因的基本质量 gsg - goodSamplesGenes(datExpr0, verbose 3) # 如果存在不合格的基因或样本直接过滤掉 if (!gsg$allOK) { datExpr0 - datExpr0[gsg$goodSamples, gsg$goodGenes] }第二步过滤表达量太低的基因以及所有样本里表达量基本不变的基因并用中位绝对偏差MAD筛出表达量变化最明显的那些基因。# 计算每个基因在所有样本中的中位数和MAD gene_median - apply(datExpr0, 2, median) gene_mad - apply(datExpr0, 2, mad) # 保留中位数表达量较高、且变异较大的基因 keep - gene_median 1 gene_mad 0.5 datExpr - datExpr0[, keep]是不是一定要做这一步如果基因数量在2万左右机器配置也足够网络构建其实也能跑完。但过滤之后模块会更稳定后续模块-性状关联的显著性也更容易出现。这里需要把握一个度如果过滤太狠可能会丢掉重要基因如果不过滤噪声又会影响结果。比较稳妥的做法是先做常规低表达过滤再用MAD筛掉后20%的基因保留约8000到15000个基因进行网络构建。2.4 样本聚类肉眼剔除异常样本这一步很多人会跳过但我觉得它比后面的参数调优更关键。样本聚类树可以直观地暴露问题如果一个样本和其他样本离得非常远说明它的整体表达模式异常可能是实验批次差异、样品污染或者数据预处理出了问题。# 样本聚类 sampleTree - hclust(dist(datExpr), method average) # 画图观察 pdf(sample_cluster.pdf, width 12, height 6) plot(sampleTree, main Sample clustering to detect outliers, sub , xlab , cex.lab 1.5, cex.axis 1.5, cex.main 2) abline(h 100, col red) dev.off()如果图里出现明显的离群样本直接用下标把它剔除# 假设样本聚类图中第3个样本明显离群 datExpr - datExpr[-3, ]cutHeight阈值不是固定不变的要看聚类树的高度分布来定。我的习惯是先用目测选一个把大多数样本聚在一起、只把极少数离群样本切出去的高度不用刻意追求统一标准。2.5 性状数据文件怎么整理性状数据是WGCNA做模块关联分析的核心输入。格式要求是“行样本列性状”样本名要和表达矩阵的行名完全一致。性状可以是连续变量比如年龄、血压、药物浓度也可以是分组变量但需要转成数值型比如疾病组1、对照组0或者多分组设计使用0/1哑变量编码。# 读取性状数据 datTraits - read.csv(sample_traits.csv, row.names 1, check.names FALSE) # 查看两个数据集的样本名交集 common_samples - intersect(rownames(datExpr), rownames(datTraits)) datExpr - datExpr[common_samples, ] datTraits - datTraits[common_samples, , drop FALSE]这里要特别注意必须保证表达矩阵和性状文件的样本一一对应。我一开始跑的时候就因为表达矩阵和性状文件样本顺序不一致直接栽过跟头。用intersect统一样本后后续分析就不会出现张冠李戴的问题。3. 软阈值选择别只会默认选9要学会看两张图3.1 无标度拓扑准则到底是什么WGCNA里最核心的一个概念是“软阈值”power。简单理解它就是一个加权系数把基因间的相关系数取绝对值的power次方得到基因之间的邻接权重。power越大弱相关被压制得越厉害网络的稀疏程度也越高。但power并不是越大越好WGCNA选择power的依据是让网络尽可能地符合无标度拓扑我们希望网络中存在少数连接度极高的“枢纽节点”而大部分节点连接度较低。衡量网络是否符合这个特征的指标是SFT.R.sq也就是无标度拓扑拟合指数。通常来说我们希望这个值尽量高尤其是要超过0.85。3.2 pickSoftThreshold的代码和输出选择软阈值不需要自己瞎试WGCNA包提供了pickSoftThreshold函数# 选择一系列候选power值 powers - c(1:10, seq(from 12, to 30, by 2)) # 计算不同power下的无标度拟合指数和平均连接度 sft - pickSoftThreshold(datExpr, powerVector powers, verbose 5) # 把结果整理成一个数据框查看 fit - sft$fitIndices print(fit[, c(Power, SFT.R.sq, mean.k., median.k.)])输出结果会看到一张类似下面的表格PowerSFT.R.sqmean.k.median.k.10.0812145.31942.120.2121023.7856.440.553366.8245.360.742178.496.280.853103.647.3100.91266.224.2120.94845.113.6眼睛不要只盯着哪个power大要同时看SFT.R.sq和mean.k.。判断标准是选择最小的、让SFT.R.sq首次进入0.85以上区间的power同时平均连接度不能太低否则网络会过度稀疏模块识别失去意义。3.3 组合图怎么看通常还会画一张双面板的组合图左边是power和拟合指数SFT.R.sq的折线右边是power和平均连接度的趋势线pdf(soft_threshold.pdf, width 10, height 5) par(mfrow c(1, 2)) cex1 - 0.9 # 左图SFT.R.sq plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], xlab Soft Threshold (power), ylab Scale Free Topology Model Fit (R^2), type n, main Scale independence) text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], labels powers, col red, cex cex1) abline(h 0.85, col red, lty 2) # 右图平均连接度 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, cex cex1) dev.off()这里有一个WGCNA代码里比较反直觉的地方左图纵轴为什么是-sign(x[,3]) * x[,2]因为当power比较低时无标度拟合指数可能是负斜率乘上负号能让所有点都朝上显示方便观察。看的时候只要看绝对值就可以了。3.4 如果R平方一直上不了0.85怎么办这是实际分析中最常见的问题之一。遇到这种情况先不要急着怀疑数据有问题。有几种处理思路检查是否没有过滤低表达基因噪声太大影响了网络结构先回去重新做数据清洗。尝试更广的power范围比如seq(1, 40, by 2)。对表达矩阵再做一次更严格的标准差筛选。样本量太小比如少于15个样本时SFT.R.sq可能很难达到0.85此时可以考虑退而求其次选择曲线转折点附近的power同时结合网络生物学可解释性来选定。我个人在临床小样本数据上遇到过几次R平方只能到0.7的情况。我的做法是选取一个中间偏上的power比如根据曲线趋势选10或者12然后看后续模块是否稳定、模块与性状的关联是否有意义。如果模块结果乱七八糟再回头调整power。4. 核心代码块网络构建与模块识别4.1 blockwiseModules函数参数逐个拆解选定power之后就到了整个WGCNA分析的重头戏——构建网络并识别模块。标准代码并不复杂核心就是blockwiseModules# 用选定的power构建网络 net - blockwiseModules( datExpr, power 9, TOMType unsigned, minModuleSize 30, reassignThreshold 0, mergeCutHeight 0.25, numericLabels TRUE, pamRespectsDendro FALSE, saveTOMs TRUE, saveTOMFileBase TOM-blockwise, verbose 3 )这个函数到底做了什么它可以拆成三层理解先计算基因间的相关性矩阵再换算成邻接矩阵然后计算TOM相似度拓扑重叠最后基于TOM相异度做层次聚类和动态剪枝得到模块。下面逐个参数说明每个参数我都会结合自己踩过的坑来讲。power就是软阈值。选9还是选12取决于第3节pickSoftThreshold的结果不要无脑用默认值。TOMType推荐用unsigned。它的意思是把负相关也当作共表达关系因为基因调控网络里存在大量负调控一个转录因子抑制下游基因时它们的表达模式是负相关的。如果改用signed则只保留正相关基因之间的共表达关系模块会更保守。minModuleSize最小模块基因数。默认值是30但是基因数少的平台数据可以降到20甚至10。设置太小会导致模块数量爆炸很多只有一两个基因的“迷你模块”没有生物学意义。mergeCutHeight相似模块合并的阈值。默认0.25意思是模块间特征基因相关性高于0.75就需要考虑合并。如果后续发现模块多而碎可以适当调低mergeCutHeight到0.2或者在后期手动做一次模块合并。numericLabels如果为TRUE模块会以数字编号命名输出的是0、1、2…… 如果为FALSE会用颜色命名比如turquoise、blue、brown。我个人建议代码阶段先用数字导出绘图时再映射成颜色这样在后续匹配基因时更方便。pamRespectsDendro是否在动态剪枝后继续用PAMPartitioning Around Medoids对模块边界做进一步细化。如果设为TRUE生成的模块更紧凑但有时会过度切割如果设为FALSE会更容易得到比较大的、稳健的模块。默认是FALSE我一般也保持FALSE。saveTOMs和saveTOMFileBase是否把TOM矩阵保存到本地。TOM矩阵是WGCNA里最大的中间产物会占很大磁盘空间但保存下来可以避免后续重复计算尤其是在做不同power对比时能节省大量时间。4.2 从net对象里能拿到什么blockwiseModules跑完之后net对象里包含了做后续分析所需的几乎所有关键结果# 查看模块划分结果 table(net$colors) # 查看模块特征基因 net$MEs # 查看模块数量 length(unique(net$colors))net$colors是一个和基因列表长度相同的向量记录了每个基因被分配到了哪个模块0表示未进入任何模块的灰色基因。net$MEs是每个样本在每个模块上的特征基因表达值后面做模块-性状关联时主要靠它。4.3 把聚类树画出来看看整体结构模块识别之后最直观的可视化方式是画聚类树和模块颜色条pdf(dendrogram_modules.pdf, width 12, height 8) plotDendroAndColors( net$dendrograms[[1]], moduleColors net$colors, groupLabels Module colors, main Gene dendrogram and module colors, dendroLabels FALSE, addGuide TRUE, hang 0.03 ) dev.off()这张图是整个WGCNA结果里我最先看的一张。理想情况下树状图会分成几个明显的大分支每个分支对应一个颜色模块。如果看到颜色条像斑马线一样反复横跳或者模块数量特别多就要考虑调整minModuleSize、mergeCutHeight或者回去检查数据预处理。刚接触WGCNA的时候我拿到模块结果第一反应是赶紧去看模块和性状的关联结果画出来的聚类树乱得没法看。后来才意识到聚类树的稳定程度本身就能反映数据质量如果这一层就有问题后面所有统计都会跟着出问题。4.4 基因数量与模块大小的合理范围一个合理的网络通常会有10到30个模块其中最大的模块基因数量占10%到30%同时会有几个中等模块和少量小模块。如果所有基因都堆在同一个巨大模块里说明power选得太低模块没有分开如果模块数量超过50个要么过滤不充分要么minModuleSize设得太小。这些标准不是绝对死线但能帮助快速判断结果是否靠谱。5. 模块合并与模块-性状关联分析5.1 为什么要做模块合并blockwiseModules里的mergeCutHeight参数会自动完成相似模块的合并。但不同版本代码或者手动调整后可能还需要自己再检查一遍。相似模块指的是两个模块的特征基因MEModule Eigengene高度相关它们本质上是同一群基因的不同变体合并后更容易解释也减少后续多重检验的次数。如果需要手动合并代码是这样的# 计算模块特征基因 MEs - moduleEigengenes(datExpr, net$colors)$eigengenes # 对模块进行合并 merge - mergeCloseModules(datExpr, net$colors, cutHeight 0.25, verbose 3) # 合并后的模块颜色 mergedColors - merge$colors # 合并后的模块特征基因 mergedMEs - merge$newMEs合并后一定要重新把模块颜色画一遍确认不要直接跳过可视化。5.2 模块-性状关联的计算方式每个模块的ME是这个模块所有基因表达模式的第一主成分代表整个模块在样本间的主要变化趋势。把ME和每个性状做相关分析就能得到模块与性状的关联矩阵。# 计算模块特征基因与性状的相关性 nSamples - nrow(datExpr) moduleTraitCor - cor(mergedMEs, datTraits, use p) moduleTraitPvalue - corPvalueStudent(moduleTraitCor, nSamples)corPvalueStudent是WGCNA为这个场景专门封装好的函数它会基于相关系数和样本量给出p值比手算方便很多。模块与性状关联的结果通常画成热图横轴是性状纵轴是模块每个格子里的数字是相关系数和括号里的p值pdf(module_trait_heatmap.pdf, width 8, height 10) textMatrix - paste(signif(moduleTraitCor, 2), \n(, signif(moduleTraitPvalue, 1), ), sep ) dim(textMatrix) - dim(moduleTraitCor) labeledHeatmap( Matrix moduleTraitCor, xLabels colnames(datTraits), yLabels colnames(mergedMEs), ySymbols colnames(mergedMEs), colorLabels FALSE, colors blueWhiteRed(50), textMatrix textMatrix, setStdMargins FALSE, cex.text 0.6, main Module-trait relationships ) dev.off()重点看两个指标相关系数的绝对值大小以及p值的显著性水平。比如某个模块与“疾病状态”的相关系数是0.72p值小于0.001那这个模块就是你后续要重点挖掘的候选模块。5.3 从模块里挑基因基因显著性与模块成员度找到一个与性状显著相关的模块后还需要知道模块里哪些基因起主导作用。WGCNA提供了两个经典指标基因显著性GS, Gene Significance基因表达量与性状之间的相关绝对值表示该基因与性状的关联强度。模块成员度MM, Module Membership基因与该模块特征基因的相关性表示该基因在模块内的核心程度。# 选择一个性状比如第1列 trait_column - 1 # 计算每个基因与性状的相关性 geneTraitSignificance - as.data.frame(cor(datExpr, datTraits[, trait_column], use p)) GS - as.numeric(geneTraitSignificance$V1) # 计算每个基因与模块特征基因的相关性 MM - as.data.frame(cor(datExpr, mergedMEs, use p))然后可以拿GS和MM做散点图筛选兼具模块核心地位且和性状显著相关的基因。比如模块内MM 0.8且GS 0.2的基因往往就是值得后续验证的候选枢纽基因。5.4 灰色模块怎么处理grey模块是WGCNA里一个特殊的存在它是所有没能被划分到任何模块的基因集合。灰色模块通常不会和性状显著相关如果有也一样要关注一下可能意味着数据里还存在另一种独立的表达模式值得单独做一次亚聚类分析。多数情况下灰色模块基因不参与后续分析。6. 实操中的零散问题与排查建议6.1 样本量太小怎么办WGCNA对样本量的最低要求保守说至少要有15到20个样本如果少于这个数基因相关性估计会非常不稳定模块结果也很难重复。真遇到小样本数据我有几个折中经验把minModuleSize调低到10左右power选择稍微偏大一点让网络更稀疏同时使用signed TOM而不是unsigned往往会稳一些。但也要明确小样本条件下的WGCNA结果只能作为探索性分析不适合直接下强结论。6.2 计算太慢、内存爆掉基因数量接近两万时blockwiseModules默认会做分块计算因为一次性计算全部基因的TOM矩阵对内存压力很大。如果仍然卡死建议先确认R是64位版本然后适当降低保留基因数比如从15000降到10000。还可以开启多线程# 启用多线程靠CPU核心数决定 enableWGCNAThreads()注意这条命令要在加载WGCNA后、构建网络之前运行而且Windows系统下多线程支持不如Linux/macOS稳定。6.3 结果可重复性问题WGCNA里有一些步骤依赖随机性尤其是样本量不特别大的时候模块识别结果可能会有轻微波动。想保证结果可重复建议在脚本开头设置随机种子set.seed(2024)另外分析过程中所有中间结果都要及时保存save(datExpr, datTraits, net, mergedColors, mergedMEs, file wgcna_step1.rda)这样即使后面改参数也不需要从最开始的读取文件重新跑一遍。6.4 包版本差异WGCNA包这些年更新不算频繁但不同小版本的默认参数可能略有差异比如blockwiseModules里的某些参数在新版本提示deprecated。如果碰到函数调用报错先看包自带的NEWS文档再看官方教程。网上很多老教程用的代码在新版本下可能会跑不通这不是你写错了需要根据报错信息微调函数名或参数。最后再分享一个偷懒技巧整套流程里最耗时的往往是不同power下的网络构建对比。我会在跑正式分析前先拿一小部分基因比如随机抽3000个基因做一个快速测试看看模块数量是否合理、聚类树是否稳定。参数基本满意后再用全部基因跑正式版本。这个小技巧能省下大量反复调试的时间。另外blockwiseModules生成的TOM文件占空间非常大分析做完如果没有特殊需要记得及时清理不然一个项目下来几百GB一点也不夸张。把模块结果、基因颜色、特征基因这些核心结果保留好就足够了。这一篇把WGCNA从数据准备到模块-性状关联的完整前半程梳理完了代码基本可以直接照着改路径和数据跑通。下一篇我会继续写模块内部的可视化细节包括基因网络导出到Cytoscape、hub gene筛选、以及怎么把模块基因批量提交给富集分析工具。跑代码的过程中遇到具体的报错和诡异结果欢迎照着这篇的排查思路先自己试一圈多数问题都出在数据格式和样本匹配上。
