在生物信息学领域，处理高通量基因表达数据时，我们常常面临一个核心挑战：如何从成千上万个基因中，识别出具有生物学意义的、协同变化的基因模块，并揭示这些模块与特定性状或疾病状态之间的关联。传统的差异表达分析虽然能找出单个基因的变化，但忽略了基因之间复杂的相互作用网络。这正是加权基因共表达网络分析（Weighted Gene Co-expression Network Analysis, WGCNA）大显身手的地方。WGCNA通过构建一个无尺度网络，将表达模式相似的基因聚类成模块，并计算模块与外部性状之间的相关性，从而系统地解析基因表达的调控模式。
本文旨在为生物信息学初学者或需要快速应用WGCNA的研究者，提供一个从零开始、可复现的完整分析流程。我们将使用R语言环境，从数据预处理、网络构建、模块识别，到模块与性状关联分析、核心基因筛选，一步步拆解WGCNA的核心步骤。即使你之前没有接触过网络分析，按照本文的步骤和解释，也能独立完成一次标准的WGCNA分析，并获得可用于后续验证和深入研究的可靠结果。
1. 理解WGCNA：从概念到工作流程
在动手写代码之前，必须先理解WGCNA要解决什么问题，以及它背后的核心逻辑。这能帮助你在后续步骤中做出正确的参数选择，并在结果出现异常时知道从哪里排查。
1.1 什么是加权基因共表达网络？
简单来说，共表达网络就是将基因视为网络中的“节点”，如果两个基因的表达模式在所有样本中高度相似（例如，总是同时上调或下调），那么它们之间就存在一条“边”。WGCNA的“加权”体现在，它并不简单地将基因关系二分为“相关”或“不相关”，而是根据基因表达相关性（通常是皮尔逊相关系数）的绝对值，通过一个幂函数（Power）进行加权转换，使得强相关的连接权重更高，弱相关的连接权重趋近于零。这种转换旨在使最终的网络符合“无尺度”拓扑特性，即网络中大部分节点连接数较少，但存在少数高度连接的枢纽节点。
1.2 WGCNA的核心分析步骤
一个标准的WGCNA分析流程通常包含以下几个关键阶段：
- 数据输入与预处理：准备基因表达矩阵和样本性状数据，并进行初步的质量控制，如去除低表达基因、处理异常样本。
- 软阈值功率（Soft Thresholding Power）选择：这是构建网络最关键的一步，目的是确定一个幂指数（β），使得网络尽可能接近无尺度拓扑结构。
- 网络构建与模块识别：基于选定的软阈值，计算基因间的邻接关系，进而得到拓扑重叠矩阵（TOM）。然后利用层次聚类和动态树切割法，将基因划分为不同的共表达模块。
- 模块与性状关联分析：计算每个模块的特征向量基因（Module Eigengene, ME）与外部样本性状（如疾病分期、临床指标）之间的相关性，找出与目标性状显著相关的模块。
- 核心基因筛选与网络可视化：在感兴趣的模块内，根据基因与模块的相关性（模块成员度，MM）和基因与性状的相关性（基因显著性，GS），筛选出模块内的核心（Hub）基因，并进行网络可视化。
1.3 分析前的关键决策点
开始前，你需要明确：
- 数据类型：WGCNA主要针对芯片或RNA-seq得到的基因表达矩阵（行是基因，列是样本）。
- 样本量：WGCNA需要一定的样本量来稳定地估计基因间的相关性。通常建议样本数不少于15-20个。样本量过小可能导致网络不稳定，结果不可靠。
- 性状数据：你需要准备与样本一一对应的性状数据，可以是连续型（如血压值）或分类型（如健康/患病）。这是后续关联分析的基础。
2. 环境准备与数据加载
我们将在一个干净的R环境中完成所有分析。请确保你已安装R（建议版本4.0以上）和RStudio。
2.1 安装必要的R包
WGCNA分析主要依赖WGCNA和flashClust包。此外，我们还会用到一些数据处理和可视化的辅助包。在R控制台或脚本中执行以下命令：
# 设置CRAN镜像，加速下载（可选，根据你的网络环境选择） # options(repos = c(CRAN = "https://mirrors.tuna.tsinghua.edu.cn/CRAN/")) # 安装BiocManager，用于安装生物信息学相关包 if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") # 安装WGCNA及其依赖 BiocManager::install("WGCNA") # 安装其他有用的包 install.packages(c("tidyverse", "reshape2", "RColorBrewer", "corrplot", "pheatmap"))安装完成后，加载这些包：
library(WGCNA) library(tidyverse) library(reshape2) library(RColorBrewer) library(corrplot) # 允许并行计算以加速TOM计算（对于大数据集非常有效） enableWGCNAThreads()2.2 准备输入数据
WGCNA需要两个核心数据文件：
- 表达矩阵（Expression Data）：一个数据框或矩阵，行名是基因标识符（如Gene Symbol或Ensembl ID），列名是样本ID。值通常是经过标准化（如FPKM、TPM）的表达量。我们假设你有一个名为
gene_expr.csv的文件。 - 性状数据（Trait Data）：一个数据框，行名是样本ID（必须与表达矩阵的列名完全一致），列是各种性状指标。我们假设你有一个名为
sample_traits.csv的文件。
让我们加载并查看数据：
# 1. 加载表达数据 expr_data <- read.csv("gene_expr.csv", row.names = 1, check.names = FALSE) # 查看数据维度：基因数 x 样本数 dim(expr_data) # 查看前几行和前几列 head(expr_data[, 1:5]) # 2. 加载性状数据 trait_data <- read.csv("sample_traits.csv", row.names = 1, check.names = FALSE) # 确保样本顺序与表达矩阵一致 trait_data <- trait_data[colnames(expr_data), ] dim(trait_data) head(trait_data)注意：
check.names = FALSE参数可以防止R将列名（样本名）中的特殊字符（如“-”）修改，确保样本ID的一致性。
2.3 数据预处理与过滤
原始表达数据中可能存在大量低表达或在不同样本间无变化的基因，这些基因对构建有生物学意义的网络贡献很小，且会极大增加计算负担。我们需要进行过滤。
# 方法一：根据均值或方差过滤（常用） # 计算每个基因在所有样本中的平均表达量 gene_mean <- apply(expr_data, 1, mean) # 计算每个基因在所有样本中的表达量方差 gene_var <- apply(expr_data, 1, var) # 绘制分布图，帮助确定阈值 par(mfrow = c(1,2)) hist(gene_mean, breaks=100, main="Gene Mean Expression", xlab="Mean") hist(gene_var, breaks=100, main="Gene Variance", xlab="Variance") # 例如，保留平均表达量大于1且方差大于0.1的基因 filtered_expr <- expr_data[gene_mean > 1 & gene_var > 0.1, ] dim(filtered_expr) # 查看过滤后的基因数量 # 方法二：使用WGCNA内置的goodSamplesGenes函数进行快速检查 gsg <- goodSamplesGenes(filtered_expr, verbose = 3) gsg$allOK # 如果为TRUE，说明数据格式基本合格 # 如果不OK，可以查看并移除有问题的基因和样本 if (!gsg$allOK) { # 打印有问题的基因和样本 if (sum(!gsg$goodGenes) > 0) printFlush(paste("Removing genes:", paste(names(filtered_expr)[!gsg$goodGenes], collapse = ", "))) if (sum(!gsg$goodSamples) > 0) printFlush(paste("Removing samples:", paste(rownames(filtered_expr)[!gsg$goodSamples], collapse = ", "))) # 移除问题行和列 filtered_expr <- filtered_expr[gsg$goodSamples, gsg$goodGenes] }预处理后，我们得到了一个相对干净的表达矩阵filtered_expr，用于后续的网络构建。
3. 构建共表达网络与识别基因模块
这是WGCNA最核心的部分，我们将通过选择软阈值、计算邻接矩阵、TOM矩阵，最终将基因聚类成模块。
3.1 选择软阈值功率（β）
软阈值功率的选择目标是使构建的网络近似无尺度拓扑。我们通过检查不同β值下网络的拓扑结构拟合指数（scale-free topology fit index, R^2）和平均连接度（mean connectivity）来决定。
# 设置一组候选的软阈值功率 powers <- c(1:10, seq(12, 30, by=2)) # 调用pickSoftThreshold函数 sft <- pickSoftThreshold(t(filtered_expr), powerVector = powers, verbose = 5, networkType = "unsigned") # 可视化结果 par(mfrow = c(1,2)) cex1 = 0.9 # 图1：拟合指数与功率的关系 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 = paste("Scale independence")) text(sft$fitIndices[,1], -sign(sft$fitIndices[,3]) * sft$fitIndices[,2], labels = powers, cex = cex1, col = "red") abline(h = 0.85, col = "red") # 通常以R^2 > 0.85作为参考线 # 图2：平均连接度与功率的关系 plot(sft$fitIndices[,1], sft$fitIndices[,5], xlab = "Soft Threshold (power)", ylab = "Mean Connectivity", type = "n", main = paste("Mean connectivity")) text(sft$fitIndices[,1], sft$fitIndices[,5], labels = powers, cex = cex1, col = "red")选择β值的原则是：在满足无尺度拓扑拟合指数（R^2）足够高（通常>0.85或0.9）的前提下，选择最小的功率。因为功率越大，网络越稀疏（连接越少），可能会丢失一些有意义的弱连接。从上图示例中，假设当power=6时，R^2首次超过0.85，且平均连接度尚可，那么我们就可以选择softPower <- 6。
3.2 一步法构建网络与识别模块
WGCNA提供了blockwiseModules函数，可以高效地一次性完成邻接矩阵、TOM计算、聚类和模块识别。对于基因数不是特别多（如<20000）的数据集，我们可以直接使用。
# 设置软阈值功率 softPower <- 6 # 设置最小模块大小（基因数），通常建议在30-100之间 minModuleSize <- 30 # 设置合并相似模块的阈值（切割树状图后，特征向量相关性高于此值的模块将被合并） mergeCutHeight <- 0.25 # 执行一步法网络构建和模块识别 net <- blockwiseModules(t(filtered_expr), power = softPower, TOMType = "unsigned", # 网络类型：无符号 minModuleSize = minModuleSize, mergeCutHeight = mergeCutHeight, numericLabels = TRUE, # 模块用数字标签（TRUE）或颜色标签（FALSE） pamRespectsDendro = FALSE, # 聚类时是否尊重树状图结构 saveTOMs = TRUE, # 保存TOM矩阵，用于后续分析 saveTOMFileBase = "MyNetworkTOM", # TOM文件前缀 verbose = 3 # 输出详细信息 )运行完成后，net对象包含了模块识别结果。最重要的两个元素是：
net$colors：一个向量，长度等于输入基因数，每个基因被分配了一个模块标签（数字或颜色）。net$MEs：模块特征向量基因（Module Eigengenes, MEs）矩阵，行是样本，列是模块。
3.3 可视化模块识别结果
# 将数字标签转换为颜色标签，便于可视化 moduleColors <- labels2colors(net$colors) # 查看模块大小（每个颜色包含的基因数） table(moduleColors) # 绘制模块聚类树状图 par(mfrow = c(1,1)) plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]], "Module colors", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)这张图的上半部分是基因的层次聚类树状图，下半部分是每个基因所属模块的颜色。理想情况下，树状图下方的颜色块应该是连续的，表明聚类效果良好。
4. 关联模块与外部性状
识别出模块后，下一步是找出哪些模块与我们关心的性状（如疾病状态、治疗反应）显著相关。
4.1 计算模块与性状的相关性
首先，我们需要计算每个模块的特征向量基因（ME）与每个性状之间的相关性。
# 确保性状数据的样本顺序与表达数据完全一致 trait_data_aligned <- trait_data[colnames(filtered_expr), ] # 计算模块特征向量基因（MEs） MEs0 <- moduleEigengenes(t(filtered_expr), moduleColors)$eigengenes # 对MEs进行排序，使其与模块颜色顺序一致 MEs <- orderMEs(MEs0) # 计算MEs与性状的相关性及p值 moduleTraitCor <- cor(MEs, trait_data_aligned, use = "p") moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(trait_data_aligned)) # 可视化相关性热图 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(trait_data_aligned), 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值（括号内）。例如，0.86\n(1e-10)表示相关系数为0.86，p值为1e-10。颜色越红表示正相关越强，越蓝表示负相关越强。通过这张图，你可以快速定位到与目标性状最相关的模块（例如，与“疾病严重程度”最正相关的“蓝色”模块）。
4.2 深入分析目标模块
假设我们发现“蓝色”模块与“疾病状态”高度正相关。接下来，我们需要在这个模块内做进一步分析。
# 定义我们感兴趣的性状（例如数据框中名为“Disease_Stage”的列） trait_of_interest <- as.data.frame(trait_data_aligned$Disease_Stage) colnames(trait_of_interest) <- "DiseaseStage" # 定义我们感兴趣的模块（例如“blue”） module_of_interest <- "blue" # 获取该模块中所有基因的列索引 module_genes <- (moduleColors == module_of_interest) # 提取该模块的表达数据 module_expr <- filtered_expr[module_genes, ] # 计算模块内基因与模块特征向量基因的相关性（模块成员度，MM） ME_of_interest <- MEs[, paste0("ME", module_of_interest)] geneModuleMembership <- as.data.frame(cor(t(module_expr), ME_of_interest, use = "p")) colnames(geneModuleMembership) <- "MM" geneModuleMembership$p.MM <- corPvalueStudent(as.matrix(geneModuleMembership), nrow(trait_data_aligned)) # 计算模块内基因与目标性状的相关性（基因显著性，GS） geneTraitSignificance <- as.data.frame(cor(t(module_expr), trait_of_interest, use = "p")) colnames(geneTraitSignificance) <- "GS" geneTraitSignificance$p.GS <- corPvalueStudent(as.matrix(geneTraitSignificance), nrow(trait_data_aligned))5. 识别核心基因与结果导出
核心基因（Hub Genes）通常指在模块内连接度最高，且与目标性状相关性也较高的基因，它们往往是模块中的关键调控因子。
5.1 筛选核心基因
一个常用的策略是同时考虑模块成员度（MM）和基因显著性（GS）。
# 将MM和GS合并到一个数据框中 module_stats <- data.frame(Gene = rownames(geneModuleMembership), MM = geneModuleMembership$MM, p.MM = geneModuleMembership$p.MM, GS = geneTraitSignificance$GS, p.GS = geneTraitSignificance$p.GS) # 筛选标准：例如，MM > 0.8 且 GS > 0.5 hub_genes <- module_stats[module_stats$MM > 0.8 & module_stats$GS > 0.5, ] # 按GS降序排列 hub_genes <- hub_genes[order(-hub_genes$GS), ] head(hub_genes)5.2 可视化模块内关系
我们可以绘制MM与GS的散点图，直观展示模块内基因与模块及性状的关系。
par(mfrow = c(1,1)) verboseScatterplot(geneModuleMembership$MM, geneTraitSignificance$GS, xlab = paste("Module Membership in", module_of_interest, "module"), ylab = paste("Gene significance for", colnames(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_of_interest) abline(h = 0.5, v = 0.8, col = "red", lty = 2) # 添加筛选阈值线图中右上角的点（高MM且高GS）就是我们筛选出的潜在核心基因。
5.3 导出结果用于下游分析
最后，将关键结果保存到文件，便于后续进行功能富集分析（如GO、KEGG）或其他验证。
# 1. 导出所有基因的模块分配信息 gene_module_annotation <- data.frame(GeneID = rownames(filtered_expr), ModuleColor = moduleColors) write.csv(gene_module_annotation, file = "WGCNA_Gene_Module_Assignment.csv", row.names = FALSE) # 2. 导出模块-性状相关性矩阵 module_trait_result <- data.frame(Module = gsub("ME", "", names(MEs)), moduleTraitCor, moduleTraitPvalue) write.csv(module_trait_result, file = "WGCNA_Module_Trait_Correlation.csv", row.names = FALSE) # 3. 导出特定模块的核心基因列表 write.csv(hub_genes, file = paste0("WGCNA_Hub_Genes_", module_of_interest, "_Module.csv"), row.names = FALSE) # 4. （可选）导出整个网络的连接度（可选，文件可能很大） # adj <- adjacency(t(filtered_expr), power = softPower, type = "unsigned") # write.csv(adj, file = "WGCNA_Adjacency_Matrix.csv") # 谨慎操作，矩阵可能巨大6. 常见问题排查与参数调优指南
WGCNA分析流程相对固定，但参数选择和数据处理中的细节常导致结果不理想。以下是几个典型问题及排查思路。
6.1 软阈值选择困难
| 问题现象 | 可能原因 | 检查与解决思路 |
|---|---|---|
| 无论选择哪个power，R^2始终低于0.8 | 1. 数据噪声过大。 2. 样本量太少，相关性估计不稳定。 3. 基因过滤过于宽松，包含大量不表达基因。 | 1. 检查数据标准化和质量控制流程。 2. 考虑增加样本量（如果可能）。 3. 尝试更严格的基因过滤（提高均值或方差阈值）。 4. 如果确实无法达到高标准，可适当降低R^2阈值（如0.7），但需在文章中说明。 |
| 平均连接度随power增加下降过快 | 选择的power可能过大，导致网络过于稀疏，丢失生物学信号。 | 在R^2达标的前提下，选择较小的power。可以观察pickSoftThreshold结果图中，平均连接度开始急剧下降的拐点，选择拐点之前的power。 |
6.2 模块数量过多或过少
| 问题现象 | 可能原因 | 检查与解决思路 |
|---|---|---|
| 模块数量非常多（如>50），且很多模块基因数很少 | minModuleSize参数设置过小。mergeCutHeight参数设置过小，模块合并不充分。 | 1. 适当增大minModuleSize（如从30调到50或100）。2. 适当增大 mergeCutHeight（如从0.25调到0.3），让相似模块更容易合并。可视化合并后的树状图，观察颜色块是否更紧凑。 |
| 模块数量非常少（如<5），模块很大 | minModuleSize参数设置过大。mergeCutHeight参数设置过大，导致不同模块被过度合并。 | 1. 适当减小minModuleSize。2. 适当减小 mergeCutHeight。 |
6.3 模块-性状无显著关联
| 问题现象 | 可能原因 | 检查与解决思路 |
|---|---|---|
| 所有模块与目标性状的相关系数都很低（绝对值<0.3），且p值不显著 | 1. 目标性状与基因表达模式确实无强关联。 2. 性状数据存在错误或类型不匹配（如将分类变量当作连续变量）。 3. 样本中存在未知批次效应，掩盖了真实关联。 | 1. 重新审视科学问题，性状选择是否合理。 2. 检查性状数据分布，分类变量应转换为因子或进行哑变量编码。 3. 检查并校正可能的批次效应（可使用 sva、limma等包）。4. 尝试不同的网络类型（ networkType参数，如signedvsunsigned）。 |
6.4 计算速度慢或内存不足
对于大型数据集（基因数>2万），一步法blockwiseModules可能消耗大量内存和时间。
- 解决方案：使用
blockwiseModules的分块（block-wise）计算功能。通过设置blocks参数，将基因集分成多个块分别计算TOM，再合并。net_bw <- blockwiseModules(t(filtered_expr), power = softPower, TOMType = "unsigned", minModuleSize = minModuleSize, mergeCutHeight = mergeCutHeight, numericLabels = TRUE, nThreads = 4, # 设置使用的CPU线程数 maxBlockSize = 5000, # 每个块的最大基因数，根据内存调整 saveTOMs = TRUE, saveTOMFileBase = "MyNetworkTOM_blockwise", verbose = 3) - 内存管理：分析完成后，使用
rm()命令及时清除中间产生的大型对象（如adjacency,TOM），并使用gc()触发垃圾回收。
7. 生产环境分析与最佳实践
在科研生产中，WGCNA分析不应是一次性的脚本运行，而应是可追溯、可复现的分析流程的一部分。
7.1 分析流程固化与版本控制
- 脚本化：将上述所有步骤整合到一个或多个R脚本中，并添加详细的注释。
- 参数外部化：将软阈值功率（
softPower）、最小模块大小（minModuleSize）等关键参数放在脚本开头的变量中，方便调整和记录。 - 版本控制：使用Git等工具管理你的分析脚本和关键结果文件。每次重要的参数调整都应有一次提交记录。
- 记录会话信息：使用
sessionInfo()函数记录分析环境的R版本和所有包的版本，这是结果可复现性的关键。sink("WGCNA_Analysis_SessionInfo.txt") sessionInfo() sink()
7.2 结果解读与验证
- 生物学合理性：筛选出的核心基因列表，必须通过文献查阅或功能富集分析（如DAVID、clusterProfiler）进行验证，看其是否富集在与研究性状相关的通路上。
- 独立数据集验证：如果条件允许，应在另一个独立的队列数据中验证核心基因的表达模式及其与性状的关联。
- 不要过度解读：WGCNA揭示的是相关性，而非因果性。共表达模块中的核心基因是重要的候选分子，但需要后续实验（如敲除、过表达）来验证其功能。
7.3 性能与鲁棒性考量
- 并行计算：始终使用
enableWGCNAThreads()开启多线程支持，大幅提升TOM计算速度。 - 随机种子：WGCNA中的层次聚类等步骤可能受随机数影响。使用
set.seed(12345)固定随机种子，确保每次运行结果一致。 - 敏感性分析：尝试微调关键参数（如
softPower± 1，mergeCutHeight± 0.05），观察核心模块和核心基因列表是否稳定。如果结果对参数极度敏感，需要谨慎解释。
7.4 下一步工作方向
完成本次基础WGCNA分析后，你可以根据研究方向深入以下工作：
- 网络可视化：使用Cytoscape等工具，导入模块内基因的拓扑重叠权重，绘制更美观、可交互的子网络图，直观展示核心基因与其他基因的互作关系。
- 模块功能分析：对每个模块的基因进行GO、KEGG功能富集分析，解读模块的生物学功能。
- 构建调控网络：整合转录因子（TF）与靶基因（TG）数据库，分析模块中是否富集特定转录因子的靶基因，推测上游调控机制。
- 与其他组学数据整合：例如，将共表达模块与甲基化模块、蛋白互作网络进行关联分析，实现多组学层面的数据整合。
WGCNA是一个强大的探索性工具，它能从全局视角提炼出基因表达的协同规律。成功的分析不仅依赖于代码的正确运行，更取决于对生物学问题的深刻理解、严谨的数据预处理、合理的参数选择以及对结果的审慎生物学解释。将本文的代码框架作为起点，结合你的具体数据不断调试和思考，才能真正让WGCNA成为你解决科研问题的利器。
