WGCNA 共表达网络分析完整流程从 count 矩阵到 hub gene 清单做了差异分析找到一堆 DEG然后呢想挖掘共表达模块、找关键调控基因WGCNA 是绕不过去的经典方法。这篇文章把完整流程从零跑一遍代码可直接复制。一、什么时候需要 WGCNA差异分析给你一张 DEG 基因列表但基因之间是孤立看待的。WGCNAWeighted Gene Co-expression Network Analysis的核心思路不同它不看单个基因而是找一起变化的基因群——共表达模块module。典型应用场景你有分组数据比如肿瘤 vs 正常想找和表型最相关的基因模块你想从模块里筛hub gene核心调控基因做下游实验验证你的 RNA-seq 样本量 15有多样本重复注意样本量太小 10做 WGCNA 结果不可靠相关性估计不稳定。建议至少 15 个样本30 更好。二、WGCNA 核心概念三句话讲清软阈值Soft Threshold, β用幂函数把基因间相关系数放大差异。β 越大弱相关越被压制网络越聚焦到强共表达关系上。目标是让网络逼近无尺度拓扑scale-free topology, R^2 0.85。模块Module用层次聚类把基因按表达模式分组。同一模块的基因高度共表达用颜色编号如 turquoise、blue、brown。模块特征值Module Eigengene, ME每个模块第一主成分代表该模块所有基因的整体表达趋势。拿 ME 和你的表型分组、临床指标等做相关性就能找到与表型最相关的模块。三、数据准备3.1 安装并加载包# WGCNA 官方包CRAN if (!requireNamespace(WGCNA, quietly TRUE)) { install.packages(WGCNA) } library(WGCNA) # 多线程加速可选 allowWGCNAThreads() # 辅助包 library(dplyr) library(ggplot2) library(pheatmap) # 关闭 R 的 multithreading 警告 options(stringsAsFactors FALSE)3.2 读取 count 矩阵# 假设你的 count 矩阵格式行为基因名列为样本 counts - read.csv(counts_matrix.csv, row.names 1) # 检查维度 dim(counts) # [1] 20000 30 # 20000 基因30 个样本 head(counts)[, 1:5]3.3 过滤低表达基因这一步很关键。低表达基因的方差极小会干扰共表达网络的构建。# 过滤标准至少在一半样本中 CPM 1 # 先算 CPMcounts per million cpm_mat - cpm(counts) # 保留满足条件的基因 keep - rowSums(cpm_mat 1) ncol(counts) / 2 expr - counts[keep, ] cat(过滤前:, nrow(counts), 基因\n) cat(过滤后:, nrow(expr), 基因\n) # 通常过滤后剩 8000-12000 个基因 # 转为 log2(CPM 1) 作为输入 expr_log - log2(cpm(expr) 1)为什么用 CPM 而不是直接 log2(count 1)CPM 做了文库大小校正不同样本间的表达量才可比。也有人用 DESeq2 的 vst() 变换效果类似。3.4 样本聚类剔除离群样本# 基于所有基因做样本间聚类 sampleTree - hclust(dist(t(expr_log)), method average) # 可视化样本聚类树 plot(sampleTree, main Sample Clustering Dendrogram, xlab , sub ) # 如果有明显的离群样本如某分支特别长手动剔除 # 比如剔除编号为 3 的样本 # expr_log - expr_log[, -3]实际操作中如果样本树状图某分支明显独立且分支长度远超其他考虑剔除。但如果只有 2-3 个样本偏离一点不建议剔除——样本量本来就宝贵。四、选择软阈值Soft Threshold β这是 WGCNA 最关键的一步。β 选大了模块太少太粗选小了模块太多太碎。4.1 计算不同 β 值的无尺度拓扑拟合度# 计算一系列 β 对应的 scale-free topology fit index powers - c(1:20) sft - pickSoftThreshold(expr_log, powerVector powers, verbose 5) # 查看结果 sft$fitIndices[, 1:3]4.2 可视化选择# Scale-free topology fit index vs β 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) text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], labels powers, cex 0.9, col red) abline(h 0.90, col blue, lty 2)选择标准R^2 首次≥ 0.85对应的 β 值最低要求最好≥ 0.90推荐同时看连接度均值不能太低右侧第二个图一般选 R^2 达标后偏小的值常见选择β 6 或 7 是最常用的。如果 R^2 在 6 就到了 0.90 以上就不要再选更大的值了。# 假设选 β 6 soft_power - 6五、构建网络与检测模块5.1 计算邻近矩阵adjacency和 TOM# 计算加权邻近矩阵用选好的 softPower adjacency - adjacency(expr_log, power soft_power) # 转为 TOMTopological Overlap Matrix矩阵 TOM - TOMsimilarity(adjacency) dissTOM - 1 - TOMTOM 是什么两个基因就算各自都和第三个基因相关但彼此不一定相关。TOM 把这种间接共享邻居的关系考虑进去比单纯的相关系数更稳健。5.2 层次聚类 动态剪切分模块# 基因层次聚类 geneTree - hclust(as.dist(dissTOM), method average) # 动态树切割核心参数deepSplit, minModuleSize, mergeCutHeight dynamicMods - cutreeDynamic(dendro geneTree, distM dissTOM, deepSplit 2, minModuleSize 30, pamRespectsDendro FALSE) # 查看模块数量和大小 table(dynamicMods)5.3 合并相似模块# 计算每个模块的特征值ME MEList - moduleEigengenes(expr_log, colors dynamicMods) MEs - MEList$eigengenes # 合并特征值相关性 0.75 的模块阈值可调 mergeCutHeight - 0.25 # 1 - 0.75 0.25 merge - mergeCloseModules(expr_log, dynamicMods, cutHeight mergeCutHeight) # 合并后的模块颜色 mergedColors - merge$colors mergedMEs - merge$newMEs cat(合并前模块数:, length(unique(dynamicMods)), \n) cat(合并后模块数:, length(unique(mergedColors)), \n)5.4 可视化模块聚类树# 绘制聚类树用颜色标注模块 plotDendroAndColors(geneTree, mergedColors[groupLabels], Module colors, dendroLabels FALSE, hang 0.03, addGuide TRUE, guideHang 0.05)六、模块与表型关联找到和你的分组/表型最相关的模块。6.1 准备表型数据# 假设你有这样的样本信息表 # sample condition stage score # S01 tumor III 8.5 # S02 normal NA 2.1 # ... sample_info - read.csv(sample_info.csv, row.names 1) # 确保样本顺序一致 sample_info - sample_info[match(colnames(expr_log), rownames(sample_info)), ] # 把分类变量转为数值WGCNA 需要数值型 sample_info$tumor_numeric - ifelse(sample_info$condition tumor, 1, 0)6.2 计算模块-表型相关性# 构建表型矩阵 traits - as.data.frame(sample_info$tumor_numeric) colnames(traits) - Tumor # 计算 ME 与表型的相关性 moduleTraitCor - cor(MEs, traits, use p) moduleTraitPvalue - corPvalueStudent(moduleTraitCor, nSamples ncol(expr_log)) # 可视化热力图 labeledHeatmap(Matrix moduleTraitCor, xLabels Tumor, yLabels names(MEs), ySymbols names(MEs), colorLabels FALSE, colors blueWhiteRed(50), textMatrix signif(moduleTraitCor, 2), setStdMargins FALSE, cex.text 0.7, zlim c(-1, 1))解读看每个模块右侧的相关系数和 p 值。比如 brown 模块与 Tumor 相关系数 r 0.85, p 1e-6说明 brown 模块基因在肿瘤中整体高表达是最值得关注的目标模块。七、提取目标模块基因 hub gene 筛选7.1 提取目标模块基因列表# 假设 brown 模块与你的表型最相关 module_genes - names(expr_log) target_module - brown module_gene_list - module_genes[mergedColors target_module] cat(target_module, 模块包含, length(module_gene_list), 个基因\n) # 导出基因列表 write.csv(data.frame(gene module_gene_list), file paste0(target_module, _genes.csv), row.names FALSE)7.2 计算 GS 和 MMGSGene Significance单个基因与表型的相关性MMModule Membership单个基因与所属模块特征值ME的相关性# 计算 GS GS - as.numeric(cor(expr_log, traits, use p)) # 计算 MM基因与模块特征值的相关性 module_colors - mergedColors ME - mergedMEs[, paste0(ME, target_module)] MM - as.numeric(cor(expr_log, ME, use p))7.3 GS vs MM 散点图# 只画目标模块的基因 inModule - module_colors target_module plot(MM[inModule], GS[inModule], xlab paste0(Module Membership (, target_module, )), ylab Gene Significance (Tumor), main paste0(target_module, module: GS vs MM), col ifelse(inModule, red, grey50), pch 16, cex 0.6) # 添加拟合线 abline(lm(GS[inModule] ~ MM[inModule]), col blue, lwd 2)Hub gene 判定右上角的基因——同时具有高 GS与表型强相关和高 MM与模块核心表达模式一致。这些就是最可能的关键调控基因。7.4 筛选并导出 hub gene# Hub gene 标准MM 0.8 且 GS 0.5阈值可调 hub_genes - data.frame( gene module_genes, MM MM, GS GS, module module_colors ) %% filter(module target_module) %% filter(MM 0.8, GS 0.5) %% arrange(desc(MM)) print(hub_genes) # 导出 write.csv(hub_genes, file paste0(target_module, _hub_genes.csv), row.names FALSE)八、GO/KEGG 功能富集推荐下游找到目标模块的基因列表后下一步通常做功能富集分析看这个模块的基因参与了哪些通路。# 使用 clusterProfiler 做 GO 富集 if (!requireNamespace(clusterProfiler, quietly TRUE)) BiocManager::install(clusterProfiler) if (!requireNamespace(org.Hs.eg.db, quietly TRUE)) BiocManager::install(org.Hs.eg.db) library(clusterProfiler) library(org.Hs.eg.db) # 基因名转 Entrez ID entrez_ids - bitr(module_gene_list, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db) # GO 富集 go_res - enrichGO(gene entrez_ids$ENTREZID, OrgDb org.Hs.eg.db, ont BP, pAdjustMethod BH, qvalueCutoff 0.05) # KEGG 富集 kegg_res - enrichKEGG(gene entrez_ids$ENTREZID, organism hsa, pAdjustMethod BH, qvalueCutoff 0.05) # 可视化 top 10 dotplot(go_res, showCategory 10) ggtitle(GO BP Enrichment) dotplot(kegg_res, showCategory 10) ggtitle(KEGG Pathway)九、常见报错处理报错原因解决方案vector size cannot be NA通常是基因名包含 NA检查rownames(expr)是否有缺失值用na.omit()或重新命名allowWGCNAThreads() 无效Windows 系统多线程限制Windows 下需要用parallel::makeCluster()手动设置或忽略该警告模块数量为 0minModuleSize 设置过大降低 minModuleSize如从 30 降到 20或降低 deepSplit模块全被合并mergeCutHeight 设置过小把 mergeCutHeight 从 0.25 提高到 0.3 或 0.4内存不足基因数量过多先过滤到 10000 基因以内或使用blockwiseModules()代替手动构建 TOM内存不够时的替代方案如果你的基因数量超过 15000手动构建 TOM 会非常吃内存。直接用 WGCNA 提供的blockwiseModules()一行代码替代第五步的所有操作net - blockwiseModules(expr_log, power soft_power, TOMType unsigned, minModuleSize 30, mergeCutHeight 0.25, verbose 3) # 结果直接包含模块颜色和 ME module_colors - net$colors MEs - net$MEs十、完整代码汇总下面是从数据读取到 hub gene 导出的完整流程可以直接复制运行记得修改文件路径# # WGCNA 完整流程 # # 1. 加载包 library(WGCNA) library(dplyr) library(ggplot2) allowWGCNAThreads() options(stringsAsFactors FALSE) # 2. 读取数据 counts - read.csv(counts_matrix.csv, row.names 1) # 3. 过滤低表达基因 log2 变换 cpm_mat - cpm(counts) keep - rowSums(cpm_mat 1) ncol(counts) / 2 expr_log - log2(cpm(counts[keep, ]) 1) # 4. 样本聚类检查 sampleTree - hclust(dist(t(expr_log)), method average) plot(sampleTree, main Sample Clustering, xlab , sub ) # 5. 选软阈值 powers - c(1:20) sft - pickSoftThreshold(expr_log, powerVector powers, verbose 5) soft_power - 6 # 根据图选 # 6. 构建网络内存充足时 adjacency - adjacency(expr_log, power soft_power) TOM - TOMsimilarity(adjacency) dissTOM - 1 - TOM geneTree - hclust(as.dist(dissTOM), method average) dynamicMods - cutreeDynamic(dendro geneTree, distM dissTOM, deepSplit 2, minModuleSize 30, pamRespectsDendro FALSE) merge - mergeCloseModules(expr_log, dynamicMods, cutHeight 0.25) mergedColors - merge$colors mergedMEs - merge$newMEs # 7. 模块-表型相关性 traits - as.data.frame(ifelse(sample_info$condition tumor, 1, 0)) colnames(traits) - Tumor moduleTraitCor - cor(mergedMEs, traits, use p) moduleTraitPvalue - corPvalueStudent(moduleTraitCor, ncol(expr_log)) # 8. 提取 hub gene以最相关模块为例 target_module - brown module_genes - names(expr_log) GS - as.numeric(cor(expr_log, traits, use p)) ME - mergedMEs[, paste0(ME, target_module)] MM - as.numeric(cor(expr_log, ME, use p)) hub_genes - data.frame(gene module_genes, MM MM, GS GS, module mergedColors) %% filter(module target_module, MM 0.8, GS 0.5) %% arrange(desc(MM)) write.csv(hub_genes, file paste0(target_module, _hub_genes.csv), row.names FALSE) cat(Hub genes saved.\n)写在最后WGCNA 的学习曲线比差异分析陡不少但它的产出——共表达模块 hub gene 候选——对后续实验设计qPCR 验证、敲低/过表达实验有直接的指导意义。几个实战建议样本量是硬门槛15 个以下慎做20 再考虑软阈值不要盲目选大的R^2 到 0.85 就够了过大会损失模块分辨率deepSplit 参数影响模块粒度值越大模块越多越细一般 2 是个平衡点hub gene 要结合生物学判断GS/MM 只是统计指标最终要回到通路和文献中去验证如果你在实际操作中遇到了具体报错欢迎在评论区留言我帮你排查。