尧图网页设计 
1. 先搞清楚WGCNA到底能帮你做什么别急着跑代码如果你手头有一堆基因表达数据想找出哪些基因是“一伙的”、一起干活或者想从海量基因里找到和某个性状比如疾病、抗逆性最相关的核心模块那WGCNAWeighted Gene Co-Expression Network Analysis就是你绕不开的工具。它不是一个简单的差异分析而是通过计算基因之间的表达相关性把成千上万个基因归类成不同的“社团”模块然后把这些模块和你关心的性状关联起来帮你从系统层面理解生物学问题。很多新手一上来就照着教程跑代码但经常卡在第一步我的数据到底适不适合做WGCNA跑出来的模块一个基因都没有或者关联分析结果全是零怎么办这篇文章不会只给你代码我会先带你理清几个关键判断点确保你的分析能跑通、结果有意义。WGCNA最核心的价值是从无监督的角度发现基因的共表达模式并建立模块与性状的量化关联。它特别适合样本量相对较多建议至少15-20个以上、想探索未知调控网络的研究场景。2. 跑通WGCNA前必须检查的三件事数据、样本和性状在安装任何包、运行任何函数之前先把这三件事确认好。很多分析失败的根本原因就在这里而不是代码写错了。2.1 数据矩阵格式、缺失值和表达量你的输入数据通常是一个矩阵行是基因列是样本。WGCNA对数据质量要求比较高。格式确保是数值型矩阵data.frame或matrix不能有字符、因子。基因名行名要唯一样本名列名要清晰。缺失值理论上WGCNA要求不能有缺失值。如果你的原始数据有缺失需要提前用适当方法填补或删除缺失严重的基因/样本。常见的做法是删除在太多样本中表达量为0或NA的基因。表达量过滤低表达或几乎不变异的基因对构建网络没有贡献反而会增加计算负担和噪音。通常可以先过滤掉在所有样本中表达量都很低例如在所有样本中CPM或FPKM小于1的基因或者方差很小的基因。注意不要一上来就用全部几万个基因去做。可以先根据方差排序选择排名靠前例如前5000或8000个变化最显著的基因进行初步分析这样能大幅缩短计算时间快速验证流程。等流程跑通后再考虑用全部基因或更大的子集。2.2 样本数量和质量是关键WGCNA构建的网络稳定性高度依赖于样本量。样本数量虽然官方没有绝对下限但普遍经验是至少需要15-20个样本。样本太少基因间的相关性估计不可靠很难得到稳定的模块。样本越多网络越稳健。样本分组与批次效应如果你的样本来自不同批次、不同处理需要特别注意。强烈的批次效应会主导共表达模式让你找到的是“批次相关模块”而非“生物学功能模块”。在分析前强烈建议先检查是否存在批次效应并使用ComBat或sva等包进行校正但需谨慎避免过度校正。输入WGCNA的数据最好是已经去除主要批次效应的数据。2.3 性状数据关联分析的“靶子”这是WGCNA分析出彩的地方也是容易出问题的地方。你需要准备一个与表达矩阵样本顺序完全一致的性状数据框data frame。性状类型可以是连续型如血压值、产量、二分类如疾病/健康或多分类。对于分类性状需要转换为0/1数值或设计矩阵。性状与样本的对应这是最常见的错误来源。务必、务必、务必检查你的表达矩阵列名样本名与性状数据框的行名是否一一对应且顺序一致。一个简单的检查命令是all(colnames(gene_data) rownames(trait_data))结果必须是TRUE。性状的数量可以同时分析多个性状。WGCNA会计算每个模块与每个性状的相关性Module-Trait Relationship帮你发现哪个模块与哪个性状最相关。3. 从零开始WGCNA标准流程拆解与实操要点假设你的数据已经通过了上一章的检查我们现在开始一步步跑流程。我会用R语言环境为例重点解释每个步骤的目的和关键参数。3.1 环境准备与数据加载首先安装并加载必要的R包。# 安装WGCNA可能需要从Bioconductor安装 if (!require(WGCNA, quietly TRUE)) { install.packages(BiocManager) BiocManager::install(WGCNA) } library(WGCNA) # 启用多线程加速计算可选但推荐 enableWGCNAThreads(nThreads 4)然后加载你的表达数据datExpr和性状数据datTraits。假设你的数据是CSV格式。# 读取表达矩阵确保第一列是基因ID并设为行名 datExpr0 - read.csv(your_expression_matrix.csv, row.names 1) # 转换为矩阵并确保所有值为数值 datExpr0 - as.matrix(datExpr0) # 读取性状数据确保第一列是样本ID并设为行名 datTraits - read.csv(your_trait_data.csv, row.names 1) # 再次检查样本顺序一致性 stopifnot(all(colnames(datExpr0) rownames(datTraits)))3.2 数据预处理与离群样本检测这一步的目的是清理数据并剔除严重偏离群体的样本这些样本会破坏网络构建。# 1. 检查缺失值过多的基因和样本 gsg - goodSamplesGenes(datExpr0, verbose 3) gsg$allOK # 如果为FALSE需要处理 if (!gsg$allOK) { # 剔除不符合要求的基因和样本 datExpr0 - datExpr0[gsg$goodGenes, gsg$goodSamples] datTraits - datTraits[gsg$goodSamples, ] } # 2. 样本聚类检测离群样本 sampleTree - hclust(dist(datExpr0), method average) # 绘制聚类树目视检查是否有单独分支很长的样本离群 par(cex 0.6) plot(sampleTree, main Sample clustering to detect outliers, sub, xlab) # 3. 如果发现离群样本可以手动设定一个高度阈值进行切除 # 例如设定高度为200这个值需要根据你的图调整 cutHeight - 200 clust - cutreeStatic(sampleTree, cutHeight cutHeight, minSize 10) # clust 0 的样本就是被判定为离群的样本 keepSamples - (clust 1) datExpr - datExpr0[, keepSamples] datTraits - datTraits[keepSamples, ]3.3 确定软阈值功率Soft Thresholding Power这是WGCNA最核心的一步决定了基因间相关性的加权程度。目的是让基因连接度分布接近无尺度网络scale-free network即大部分基因连接少少数基因是高度连接的枢纽hub gene。# 选择一组软阈值进行测试 powers - c(1:20) # 调用网络拓扑分析函数 sft - pickSoftThreshold(datExpr, powerVector powers, verbose 5, networkType signed) # 绘制结果图辅助选择 par(mfrow c(1,2)) # 图1无尺度拓扑拟合指数Scale Free Topology Model Fit随软阈值的变化 plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlabSoft Threshold (power), ylabScale Free Topology Model Fit,signed R^2, typen, main paste(Scale independence)) text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labelspowers, colred) # 通常建议选择R^2首次达到0.8或0.9以上的最小power值 # 图2平均连接度Mean Connectivity随软阈值的变化 plot(sft$fitIndices[,1], sft$fitIndices[,5], xlabSoft Threshold (power), ylabMean Connectivity, typen, main paste(Mean connectivity)) text(sft$fitIndices[,1], sft$fitIndices[,5], labelspowers, colred)如何选择优先看左图Scale independence。选择第一个使R^2达到0.8或0.9以上的power值。同时参考右图平均连接度不能下降得太快不能接近0。如果power值需要设置得非常大比如20才能达到0.8可能意味着你的数据不太适合构建无尺度网络需要回头检查数据质量。通常power值在6到12之间比较常见。3.4 一步法构建网络与识别模块确定了软阈值假设我们选择power10后就可以构建网络并将基因划分到模块中。blockwiseModules函数是主流方法尤其适合基因数很多5000的情况它能分块计算节省内存。# 设置随机种子保证结果可重复 set.seed(12345) # 一步法构建网络和模块 net - blockwiseModules(datExpr, power 10, # 替换为你选定的power值 TOMType signed, # 使用有符号的TOM矩阵 minModuleSize 30, # 模块最小基因数可根据数据调整通常30-100 mergeCutHeight 0.25, # 模块合并阈值值越小模块越不易合并 numericLabels TRUE, # 模块用数字标签 pamRespectsDendro FALSE, saveTOMs TRUE, # 保存TOM矩阵用于后续分析 saveTOMFileBase MyNetworkTOM, verbose 3) # 查看模块数量及大小 table(net$colors)关键参数解释minModuleSize模块最少包含的基因数。设太小会产生很多琐碎的小模块设太大会合并掉有生物学意义的小模块。可以从30开始尝试。mergeCutHeight模块树状图切割高度用于合并相似度高的模块。值越小模块越不容易被合并得到的模块数可能越多。默认0.25是个不错的起点。numericLabelsTRUE时模块用数字0,1,2...表示其中0代表未被分配到任何模块的基因灰色模块。3.5 可视化模块与关联性状得到模块后我们需要看看它们长什么样以及和性状的关系。# 1. 将数字标签转换为颜色标签便于可视化 moduleColors - labels2colors(net$colors) # 2. 绘制模块聚类树状图 plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], Module colors, dendroLabels FALSE, hang 0.03, addGuide TRUE, guideHang 0.05) # 3. 计算模块特征基因Module Eigengene, ME MEs - net$MEs # 为MEs命名将数字前缀改为“ME” moduleLabels - net$colors MEs0 - moduleEigengenes(datExpr, moduleColors)$eigengenes MEs - orderMEs(MEs0) # 4. 计算模块与性状的相关性 moduleTraitCor - cor(MEs, datTraits, use p) moduleTraitPvalue - corPvalueStudent(moduleTraitCor, nSamples ncol(datExpr)) # 5. 绘制模块-性状关系热图 textMatrix - paste(signif(moduleTraitCor, 2), \n(, signif(moduleTraitPvalue, 1), ), sep ) dim(textMatrix) - dim(moduleTraitCor) par(mar c(6, 8.5, 3, 3)) labeledHeatmap(Matrix moduleTraitCor, xLabels names(datTraits), yLabels names(MEs), ySymbols names(MEs), colorLabels FALSE, colors blueWhiteRed(50), textMatrix textMatrix, setStdMargins FALSE, cex.text 0.5, zlim c(-1,1), main paste(Module-trait relationships))热图中每个格子显示了相关性系数和p值括号内。颜色越红表示正相关越强越蓝表示负相关越强。重点关注那些与目标性状相关性高且p值显著的模块比如与“疾病严重程度”显著正相关的“蓝色模块”。4. 结果解读与下游分析找到核心基因与功能分析跑完了图也画了接下来才是真正产生价值的步骤解读。4.1 定位关键模块与核心枢纽基因Hub Genes假设我们发现“蓝色模块”MEblue与我们的目标性状最相关。# 定义我们感兴趣的性状例如数据框中名为“Disease_Score”的列 trait_of_interest - Disease_Score # 找到与该性状最相关的模块 modNames - substring(names(MEs), 3) # 去掉ME前缀 geneModuleMembership - as.data.frame(cor(datExpr, MEs, use p)) colnames(geneModuleMembership) - paste(MM, modNames, sep) geneTraitSignificance - as.data.frame(cor(datExpr, datTraits[[trait_of_interest]], use p)) colnames(geneTraitSignificance) - paste(GS., trait_of_interest, sep) # 提取蓝色模块的基因 module - blue moduleGenes - moduleColors module # 绘制基因重要性散点图模块成员度MM vs 基因性状显著性GS par(mfrowc(1,1)) verboseScatterplot(abs(geneModuleMembership[moduleGenes, paste0(MM, module)]), abs(geneTraitSignificance[moduleGenes, 1]), xlab paste(Module Membership in, module, module), ylab paste(Gene significance for, trait_of_interest), main paste(Module membership vs. gene significance\n), cex.main 1.2, cex.lab 1.2, cex.axis 1.2, col module)在这个散点图中右上角的基因既是模块的核心成员与模块特征基因高度相关又与目标性状高度相关它们就是潜在的枢纽基因Hub Genes是后续实验验证的优先候选。4.2 模块的功能富集分析知道了关键模块是哪个下一步是理解这个模块的基因集合在生物学上意味着什么。我们需要做功能富集分析GO、KEGG等。# 获取蓝色模块的所有基因ID假设行名是基因ID blue_gene_ids - rownames(datExpr)[moduleColors blue] # 将基因ID写入文件用于后续在线工具如DAVID、Metascape或R包如clusterProfiler分析 write.table(blue_gene_ids, file blue_module_genes.txt, quote FALSE, row.names FALSE, col.names FALSE)不要只做最相关模块。建议对排名前几的模块都进行富集分析有时一个性状可能由多个生物学通路协同影响。4.3 导出网络用于Cytoscape可视化如果你想深入探索模块内部的基因互作关系可以将拓扑重叠矩阵TOM导出用Cytoscape等软件进行网络可视化。# 重新计算整个网络的TOM如果之前保存了可以加载 load(net$TOMFiles[1]) # 加载TOM矩阵 TOM - as.matrix(TOM) # 选择蓝色模块的基因 blue_module_indices - which(moduleColors blue) blue_gene_names - rownames(datExpr)[blue_module_indices] # 提取蓝色模块对应的子TOM矩阵 blue_TOM - TOM[blue_module_indices, blue_module_indices] rownames(blue_TOM) - blue_gene_names colnames(blue_TOM) - blue_gene_names # 导出为Cytoscape可读的格式例如边列表 # 这里需要将TOM矩阵转换为边列表并设定一个连接度阈值例如TOM 0.1 library(igraph) blue_adj_matrix - blue_TOM diag(blue_adj_matrix) - 0 # 将对角线自连接设为0 # 创建图对象只保留权重较高的边 g - graph.adjacency(blue_adj_matrix, modeundirected, weightedTRUE, diagFALSE) g - simplify(g) # 去除可能的重复边 # 可以按权重过滤边例如只保留权重前500的边 edge_weights - E(g)$weight threshold - sort(edge_weights, decreasingTRUE)[500] g_sub - delete_edges(g, E(g)[weight threshold]) # 将边列表和节点列表写入文件 write_graph(g_sub, fileblue_module_network.graphml, formatgraphml) # 也可以导出为简单的边列表 edge_list - get.edgelist(g_sub) edge_list_with_weight - cbind(edge_list, E(g_sub)$weight) write.table(edge_list_with_weight, fileblue_module_edgelist.txt, sep\t, quoteFALSE, row.namesFALSE, col.namesFALSE)在Cytoscape中你可以直观地看到模块内部的连接情况找出处于网络中心位置连接数多的枢纽基因。5. 常见报错、排查与进阶思考跑WGCNA时你大概率会遇到下面这些问题。5.1 报错 “Error in cor(x, y, use ‘p’) : ‘y’ must be numeric”原因你的性状数据datTraits中包含了非数值列如字符型的样本分组。解决检查str(datTraits)。将分类变量转换为数值。例如将“Control”和“Case”转换为0和1或者使用model.matrix()创建设计矩阵。5.2 模块-性状热图全是灰色或不显著原因1样本量太少统计效力不足。解决增加样本量是根本。如果无法增加需谨慎解读结果或考虑使用其他更适合小样本的分析方法。原因2性状与基因表达确实没有强关联。解决这是可能的生物学事实。检查你的性状测量是否准确或者考虑其他类型的性状如临床指标、代谢物数据。原因3数据预处理不当批次效应过强。解决重新检查并校正批次效应。5.3 运行blockwiseModules时内存不足或时间过长原因基因数太多如20000一次性计算TOM矩阵非常消耗内存。解决使用blockwiseModules函数它本身就是为大数据设计的。在函数内设置maxBlockSize参数指定每个块的最大基因数例如5000。函数会自动分块计算。在第一步就进行更严格的基因过滤只保留方差最大的几千个基因进行初步分析。使用更高内存的计算机或服务器。5.4 得到的模块太多或太少调整minModuleSize增加此值如从30调到50会减少小模块使模块总数变少。调整mergeCutHeight增加此值如从0.25调到0.3会使更多相似模块被合并模块总数变少减小则相反。调整deepSplit参数在blockwiseModules中deepSplit控制树状图切割的深度取值0-4。值越大切割越细模块越多。可以尝试设置为2或3。5.5 进阶思考WGCNA结果的可靠性WGCNA是一个强大的探索性工具但它给出的结果是“相关关系”而非“因果关系”。模块与性状相关不代表该模块的基因直接调控该性状。后续必须通过实验验证对筛选出的枢纽基因进行敲除、过表达等实验。与其他数据整合如与ChIP-seq转录因子结合、ATAC-seq染色质开放性数据整合寻找上游调控证据。因果推断方法如使用孟德尔随机化等方法进行因果探索。最后也是最关键的建议不要只满足于跑通流程和画出漂亮的图。把分析脚本、中间文件、参数选择理由都记录清楚。对于关键结果如枢纽基因列表用独立的数据集如果有或通过文献检索进行交叉验证。WGCNA是一个起点它帮你从数据海洋中捞出几条“大鱼”但判断这些鱼到底是什么、怎么吃还需要更深入的生物学知识和后续工作。 WGCNA实战指南：从数据检查到模块分析，避免新手常见错误
1. 先搞清楚WGCNA到底能帮你做什么别急着跑代码如果你手头有一堆基因表达数据想找出哪些基因是“一伙的”、一起干活或者想从海量基因里找到和某个性状比如疾病、抗逆性最相关的核心模块那WGCNAWeighted Gene Co-Expression Network Analysis就是你绕不开的工具。它不是一个简单的差异分析而是通过计算基因之间的表达相关性把成千上万个基因归类成不同的“社团”模块然后把这些模块和你关心的性状关联起来帮你从系统层面理解生物学问题。很多新手一上来就照着教程跑代码但经常卡在第一步我的数据到底适不适合做WGCNA跑出来的模块一个基因都没有或者关联分析结果全是零怎么办这篇文章不会只给你代码我会先带你理清几个关键判断点确保你的分析能跑通、结果有意义。WGCNA最核心的价值是从无监督的角度发现基因的共表达模式并建立模块与性状的量化关联。它特别适合样本量相对较多建议至少15-20个以上、想探索未知调控网络的研究场景。2. 跑通WGCNA前必须检查的三件事数据、样本和性状在安装任何包、运行任何函数之前先把这三件事确认好。很多分析失败的根本原因就在这里而不是代码写错了。2.1 数据矩阵格式、缺失值和表达量你的输入数据通常是一个矩阵行是基因列是样本。WGCNA对数据质量要求比较高。格式确保是数值型矩阵data.frame或matrix不能有字符、因子。基因名行名要唯一样本名列名要清晰。缺失值理论上WGCNA要求不能有缺失值。如果你的原始数据有缺失需要提前用适当方法填补或删除缺失严重的基因/样本。常见的做法是删除在太多样本中表达量为0或NA的基因。表达量过滤低表达或几乎不变异的基因对构建网络没有贡献反而会增加计算负担和噪音。通常可以先过滤掉在所有样本中表达量都很低例如在所有样本中CPM或FPKM小于1的基因或者方差很小的基因。注意不要一上来就用全部几万个基因去做。可以先根据方差排序选择排名靠前例如前5000或8000个变化最显著的基因进行初步分析这样能大幅缩短计算时间快速验证流程。等流程跑通后再考虑用全部基因或更大的子集。2.2 样本数量和质量是关键WGCNA构建的网络稳定性高度依赖于样本量。样本数量虽然官方没有绝对下限但普遍经验是至少需要15-20个样本。样本太少基因间的相关性估计不可靠很难得到稳定的模块。样本越多网络越稳健。样本分组与批次效应如果你的样本来自不同批次、不同处理需要特别注意。强烈的批次效应会主导共表达模式让你找到的是“批次相关模块”而非“生物学功能模块”。在分析前强烈建议先检查是否存在批次效应并使用ComBat或sva等包进行校正但需谨慎避免过度校正。输入WGCNA的数据最好是已经去除主要批次效应的数据。2.3 性状数据关联分析的“靶子”这是WGCNA分析出彩的地方也是容易出问题的地方。你需要准备一个与表达矩阵样本顺序完全一致的性状数据框data frame。性状类型可以是连续型如血压值、产量、二分类如疾病/健康或多分类。对于分类性状需要转换为0/1数值或设计矩阵。性状与样本的对应这是最常见的错误来源。务必、务必、务必检查你的表达矩阵列名样本名与性状数据框的行名是否一一对应且顺序一致。一个简单的检查命令是all(colnames(gene_data) rownames(trait_data))结果必须是TRUE。性状的数量可以同时分析多个性状。WGCNA会计算每个模块与每个性状的相关性Module-Trait Relationship帮你发现哪个模块与哪个性状最相关。3. 从零开始WGCNA标准流程拆解与实操要点假设你的数据已经通过了上一章的检查我们现在开始一步步跑流程。我会用R语言环境为例重点解释每个步骤的目的和关键参数。3.1 环境准备与数据加载首先安装并加载必要的R包。# 安装WGCNA可能需要从Bioconductor安装 if (!require(WGCNA, quietly TRUE)) { install.packages(BiocManager) BiocManager::install(WGCNA) } library(WGCNA) # 启用多线程加速计算可选但推荐 enableWGCNAThreads(nThreads 4)然后加载你的表达数据datExpr和性状数据datTraits。假设你的数据是CSV格式。# 读取表达矩阵确保第一列是基因ID并设为行名 datExpr0 - read.csv(your_expression_matrix.csv, row.names 1) # 转换为矩阵并确保所有值为数值 datExpr0 - as.matrix(datExpr0) # 读取性状数据确保第一列是样本ID并设为行名 datTraits - read.csv(your_trait_data.csv, row.names 1) # 再次检查样本顺序一致性 stopifnot(all(colnames(datExpr0) rownames(datTraits)))3.2 数据预处理与离群样本检测这一步的目的是清理数据并剔除严重偏离群体的样本这些样本会破坏网络构建。# 1. 检查缺失值过多的基因和样本 gsg - goodSamplesGenes(datExpr0, verbose 3) gsg$allOK # 如果为FALSE需要处理 if (!gsg$allOK) { # 剔除不符合要求的基因和样本 datExpr0 - datExpr0[gsg$goodGenes, gsg$goodSamples] datTraits - datTraits[gsg$goodSamples, ] } # 2. 样本聚类检测离群样本 sampleTree - hclust(dist(datExpr0), method average) # 绘制聚类树目视检查是否有单独分支很长的样本离群 par(cex 0.6) plot(sampleTree, main Sample clustering to detect outliers, sub, xlab) # 3. 如果发现离群样本可以手动设定一个高度阈值进行切除 # 例如设定高度为200这个值需要根据你的图调整 cutHeight - 200 clust - cutreeStatic(sampleTree, cutHeight cutHeight, minSize 10) # clust 0 的样本就是被判定为离群的样本 keepSamples - (clust 1) datExpr - datExpr0[, keepSamples] datTraits - datTraits[keepSamples, ]3.3 确定软阈值功率Soft Thresholding Power这是WGCNA最核心的一步决定了基因间相关性的加权程度。目的是让基因连接度分布接近无尺度网络scale-free network即大部分基因连接少少数基因是高度连接的枢纽hub gene。# 选择一组软阈值进行测试 powers - c(1:20) # 调用网络拓扑分析函数 sft - pickSoftThreshold(datExpr, powerVector powers, verbose 5, networkType signed) # 绘制结果图辅助选择 par(mfrow c(1,2)) # 图1无尺度拓扑拟合指数Scale Free Topology Model Fit随软阈值的变化 plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlabSoft Threshold (power), ylabScale Free Topology Model Fit,signed R^2, typen, main paste(Scale independence)) text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labelspowers, colred) # 通常建议选择R^2首次达到0.8或0.9以上的最小power值 # 图2平均连接度Mean Connectivity随软阈值的变化 plot(sft$fitIndices[,1], sft$fitIndices[,5], xlabSoft Threshold (power), ylabMean Connectivity, typen, main paste(Mean connectivity)) text(sft$fitIndices[,1], sft$fitIndices[,5], labelspowers, colred)如何选择优先看左图Scale independence。选择第一个使R^2达到0.8或0.9以上的power值。同时参考右图平均连接度不能下降得太快不能接近0。如果power值需要设置得非常大比如20才能达到0.8可能意味着你的数据不太适合构建无尺度网络需要回头检查数据质量。通常power值在6到12之间比较常见。3.4 一步法构建网络与识别模块确定了软阈值假设我们选择power10后就可以构建网络并将基因划分到模块中。blockwiseModules函数是主流方法尤其适合基因数很多5000的情况它能分块计算节省内存。# 设置随机种子保证结果可重复 set.seed(12345) # 一步法构建网络和模块 net - blockwiseModules(datExpr, power 10, # 替换为你选定的power值 TOMType signed, # 使用有符号的TOM矩阵 minModuleSize 30, # 模块最小基因数可根据数据调整通常30-100 mergeCutHeight 0.25, # 模块合并阈值值越小模块越不易合并 numericLabels TRUE, # 模块用数字标签 pamRespectsDendro FALSE, saveTOMs TRUE, # 保存TOM矩阵用于后续分析 saveTOMFileBase MyNetworkTOM, verbose 3) # 查看模块数量及大小 table(net$colors)关键参数解释minModuleSize模块最少包含的基因数。设太小会产生很多琐碎的小模块设太大会合并掉有生物学意义的小模块。可以从30开始尝试。mergeCutHeight模块树状图切割高度用于合并相似度高的模块。值越小模块越不容易被合并得到的模块数可能越多。默认0.25是个不错的起点。numericLabelsTRUE时模块用数字0,1,2...表示其中0代表未被分配到任何模块的基因灰色模块。3.5 可视化模块与关联性状得到模块后我们需要看看它们长什么样以及和性状的关系。# 1. 将数字标签转换为颜色标签便于可视化 moduleColors - labels2colors(net$colors) # 2. 绘制模块聚类树状图 plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], Module colors, dendroLabels FALSE, hang 0.03, addGuide TRUE, guideHang 0.05) # 3. 计算模块特征基因Module Eigengene, ME MEs - net$MEs # 为MEs命名将数字前缀改为“ME” moduleLabels - net$colors MEs0 - moduleEigengenes(datExpr, moduleColors)$eigengenes MEs - orderMEs(MEs0) # 4. 计算模块与性状的相关性 moduleTraitCor - cor(MEs, datTraits, use p) moduleTraitPvalue - corPvalueStudent(moduleTraitCor, nSamples ncol(datExpr)) # 5. 绘制模块-性状关系热图 textMatrix - paste(signif(moduleTraitCor, 2), \n(, signif(moduleTraitPvalue, 1), ), sep ) dim(textMatrix) - dim(moduleTraitCor) par(mar c(6, 8.5, 3, 3)) labeledHeatmap(Matrix moduleTraitCor, xLabels names(datTraits), yLabels names(MEs), ySymbols names(MEs), colorLabels FALSE, colors blueWhiteRed(50), textMatrix textMatrix, setStdMargins FALSE, cex.text 0.5, zlim c(-1,1), main paste(Module-trait relationships))热图中每个格子显示了相关性系数和p值括号内。颜色越红表示正相关越强越蓝表示负相关越强。重点关注那些与目标性状相关性高且p值显著的模块比如与“疾病严重程度”显著正相关的“蓝色模块”。4. 结果解读与下游分析找到核心基因与功能分析跑完了图也画了接下来才是真正产生价值的步骤解读。4.1 定位关键模块与核心枢纽基因Hub Genes假设我们发现“蓝色模块”MEblue与我们的目标性状最相关。# 定义我们感兴趣的性状例如数据框中名为“Disease_Score”的列 trait_of_interest - Disease_Score # 找到与该性状最相关的模块 modNames - substring(names(MEs), 3) # 去掉ME前缀 geneModuleMembership - as.data.frame(cor(datExpr, MEs, use p)) colnames(geneModuleMembership) - paste(MM, modNames, sep) geneTraitSignificance - as.data.frame(cor(datExpr, datTraits[[trait_of_interest]], use p)) colnames(geneTraitSignificance) - paste(GS., trait_of_interest, sep) # 提取蓝色模块的基因 module - blue moduleGenes - moduleColors module # 绘制基因重要性散点图模块成员度MM vs 基因性状显著性GS par(mfrowc(1,1)) verboseScatterplot(abs(geneModuleMembership[moduleGenes, paste0(MM, module)]), abs(geneTraitSignificance[moduleGenes, 1]), xlab paste(Module Membership in, module, module), ylab paste(Gene significance for, trait_of_interest), main paste(Module membership vs. gene significance\n), cex.main 1.2, cex.lab 1.2, cex.axis 1.2, col module)在这个散点图中右上角的基因既是模块的核心成员与模块特征基因高度相关又与目标性状高度相关它们就是潜在的枢纽基因Hub Genes是后续实验验证的优先候选。4.2 模块的功能富集分析知道了关键模块是哪个下一步是理解这个模块的基因集合在生物学上意味着什么。我们需要做功能富集分析GO、KEGG等。# 获取蓝色模块的所有基因ID假设行名是基因ID blue_gene_ids - rownames(datExpr)[moduleColors blue] # 将基因ID写入文件用于后续在线工具如DAVID、Metascape或R包如clusterProfiler分析 write.table(blue_gene_ids, file blue_module_genes.txt, quote FALSE, row.names FALSE, col.names FALSE)不要只做最相关模块。建议对排名前几的模块都进行富集分析有时一个性状可能由多个生物学通路协同影响。4.3 导出网络用于Cytoscape可视化如果你想深入探索模块内部的基因互作关系可以将拓扑重叠矩阵TOM导出用Cytoscape等软件进行网络可视化。# 重新计算整个网络的TOM如果之前保存了可以加载 load(net$TOMFiles[1]) # 加载TOM矩阵 TOM - as.matrix(TOM) # 选择蓝色模块的基因 blue_module_indices - which(moduleColors blue) blue_gene_names - rownames(datExpr)[blue_module_indices] # 提取蓝色模块对应的子TOM矩阵 blue_TOM - TOM[blue_module_indices, blue_module_indices] rownames(blue_TOM) - blue_gene_names colnames(blue_TOM) - blue_gene_names # 导出为Cytoscape可读的格式例如边列表 # 这里需要将TOM矩阵转换为边列表并设定一个连接度阈值例如TOM 0.1 library(igraph) blue_adj_matrix - blue_TOM diag(blue_adj_matrix) - 0 # 将对角线自连接设为0 # 创建图对象只保留权重较高的边 g - graph.adjacency(blue_adj_matrix, modeundirected, weightedTRUE, diagFALSE) g - simplify(g) # 去除可能的重复边 # 可以按权重过滤边例如只保留权重前500的边 edge_weights - E(g)$weight threshold - sort(edge_weights, decreasingTRUE)[500] g_sub - delete_edges(g, E(g)[weight threshold]) # 将边列表和节点列表写入文件 write_graph(g_sub, fileblue_module_network.graphml, formatgraphml) # 也可以导出为简单的边列表 edge_list - get.edgelist(g_sub) edge_list_with_weight - cbind(edge_list, E(g_sub)$weight) write.table(edge_list_with_weight, fileblue_module_edgelist.txt, sep\t, quoteFALSE, row.namesFALSE, col.namesFALSE)在Cytoscape中你可以直观地看到模块内部的连接情况找出处于网络中心位置连接数多的枢纽基因。5. 常见报错、排查与进阶思考跑WGCNA时你大概率会遇到下面这些问题。5.1 报错 “Error in cor(x, y, use ‘p’) : ‘y’ must be numeric”原因你的性状数据datTraits中包含了非数值列如字符型的样本分组。解决检查str(datTraits)。将分类变量转换为数值。例如将“Control”和“Case”转换为0和1或者使用model.matrix()创建设计矩阵。5.2 模块-性状热图全是灰色或不显著原因1样本量太少统计效力不足。解决增加样本量是根本。如果无法增加需谨慎解读结果或考虑使用其他更适合小样本的分析方法。原因2性状与基因表达确实没有强关联。解决这是可能的生物学事实。检查你的性状测量是否准确或者考虑其他类型的性状如临床指标、代谢物数据。原因3数据预处理不当批次效应过强。解决重新检查并校正批次效应。5.3 运行blockwiseModules时内存不足或时间过长原因基因数太多如20000一次性计算TOM矩阵非常消耗内存。解决使用blockwiseModules函数它本身就是为大数据设计的。在函数内设置maxBlockSize参数指定每个块的最大基因数例如5000。函数会自动分块计算。在第一步就进行更严格的基因过滤只保留方差最大的几千个基因进行初步分析。使用更高内存的计算机或服务器。5.4 得到的模块太多或太少调整minModuleSize增加此值如从30调到50会减少小模块使模块总数变少。调整mergeCutHeight增加此值如从0.25调到0.3会使更多相似模块被合并模块总数变少减小则相反。调整deepSplit参数在blockwiseModules中deepSplit控制树状图切割的深度取值0-4。值越大切割越细模块越多。可以尝试设置为2或3。5.5 进阶思考WGCNA结果的可靠性WGCNA是一个强大的探索性工具但它给出的结果是“相关关系”而非“因果关系”。模块与性状相关不代表该模块的基因直接调控该性状。后续必须通过实验验证对筛选出的枢纽基因进行敲除、过表达等实验。与其他数据整合如与ChIP-seq转录因子结合、ATAC-seq染色质开放性数据整合寻找上游调控证据。因果推断方法如使用孟德尔随机化等方法进行因果探索。最后也是最关键的建议不要只满足于跑通流程和画出漂亮的图。把分析脚本、中间文件、参数选择理由都记录清楚。对于关键结果如枢纽基因列表用独立的数据集如果有或通过文献检索进行交叉验证。WGCNA是一个起点它帮你从数据海洋中捞出几条“大鱼”但判断这些鱼到底是什么、怎么吃还需要更深入的生物学知识和后续工作。