企业官网建设流程全解析

1. 先搞清楚WGCNA到底能帮你做什么，别急着跑代码

如果你手头有一堆基因表达数据，想找出哪些基因是“一伙的”、一起干活，或者想从海量基因里找到和某个性状（比如疾病、抗逆性）最相关的核心模块，那WGCNA（Weighted 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], 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, col="red") # 通常建议选择R^2首次达到0.8或0.9以上的最小power值 # 图2：平均连接度（Mean Connectivity）随软阈值的变化 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, col="red")

如何选择？优先看左图（Scale independence）。选择第一个使R^2达到0.8或0.9以上的power值。同时参考右图，平均连接度不能下降得太快（不能接近0）。如果power值需要设置得非常大（比如>20）才能达到0.8，可能意味着你的数据不太适合构建无尺度网络，需要回头检查数据质量。通常，power值在6到12之间比较常见。

3.4 一步法构建网络与识别模块

确定了软阈值（假设我们选择power=10）后，就可以构建网络并将基因划分到模块中。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是个不错的起点。
• numericLabels：TRUE时模块用数字（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(mfrow=c(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, mode="undirected", weighted=TRUE, diag=FALSE) g <- simplify(g) # 去除可能的重复边 # 可以按权重过滤边，例如只保留权重前500的边 edge_weights <- E(g)$weight threshold <- sort(edge_weights, decreasing=TRUE)[500] g_sub <- delete_edges(g, E(g)[weight < threshold]) # 将边列表和节点列表写入文件 write_graph(g_sub, file="blue_module_network.graphml", format="graphml") # 也可以导出为简单的边列表 edge_list <- get.edgelist(g_sub) edge_list_with_weight <- cbind(edge_list, E(g_sub)$weight) write.table(edge_list_with_weight, file="blue_module_edgelist.txt", sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)

在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矩阵非常消耗内存。
• 解决：
  1. 使用blockwiseModules函数，它本身就是为大数据设计的。
  2. 在函数内设置maxBlockSize参数，指定每个块的最大基因数（例如5000）。函数会自动分块计算。
  3. 在第一步就进行更严格的基因过滤，只保留方差最大的几千个基因进行初步分析。
  4. 使用更高内存的计算机或服务器。

5.4 得到的模块太多或太少

• 调整minModuleSize：增加此值（如从30调到50）会减少小模块，使模块总数变少。
• 调整mergeCutHeight：增加此值（如从0.25调到0.3）会使更多相似模块被合并，模块总数变少；减小则相反。
• 调整deepSplit参数：在blockwiseModules中，deepSplit控制树状图切割的深度，取值0-4。值越大，切割越细，模块越多。可以尝试设置为2或3。

5.5 进阶思考：WGCNA结果的可靠性

WGCNA是一个强大的探索性工具，但它给出的结果是“相关关系”，而非“因果关系”。模块与性状相关，不代表该模块的基因直接调控该性状。后续必须通过：

1. 实验验证：对筛选出的枢纽基因进行敲除、过表达等实验。
2. 与其他数据整合：如与ChIP-seq（转录因子结合）、ATAC-seq（染色质开放性）数据整合，寻找上游调控证据。
3. 因果推断方法：如使用孟德尔随机化等方法进行因果探索。

最后，也是最关键的建议：不要只满足于跑通流程和画出漂亮的图。把分析脚本、中间文件、参数选择理由都记录清楚。对于关键结果（如枢纽基因列表），用独立的数据集（如果有）或通过文献检索进行交叉验证。WGCNA是一个起点，它帮你从数据海洋中捞出几条“大鱼”，但判断这些鱼到底是什么、怎么吃，还需要更深入的生物学知识和后续工作。