企业官网建设流程全解析

用WGCNA在R中挖掘RNA-seq数据的基因社交网络：从入门到实战

当我们拿到一份RNA-seq差异表达分析结果时，常常会好奇这些基因之间是否存在某种"社交关系"——哪些基因喜欢一起"活动"？哪些基因是社交圈中的"核心人物"？这正是WGCNA（加权基因共表达网络分析）能够回答的问题。想象一下，如果把基因比作人群，WGCNA就是帮我们找出哪些人经常一起出现（共表达），并识别出每个朋友圈中的关键人物（hub基因）。下面我们就用R语言一步步实现这个分析过程。

1. 准备工作与环境搭建

在开始分析之前，我们需要确保R环境和必要的包已经准备就绪。WGCNA分析对计算资源有一定要求，特别是处理大规模基因表达数据时。建议使用至少16GB内存的计算机进行分析。

首先安装必要的R包：

if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install(c("WGCNA", "DESeq2", "ggplot2"))

加载所需库：

library(WGCNA) library(DESeq2) library(ggplot2)

提示：WGCNA包在Bioconductor上维护，建议使用BiocManager进行安装以确保版本兼容性。

对于RNA-seq数据，我们通常从原始计数数据开始。假设你已经有了一个DESeq2的数据对象dds，我们可以先进行方差稳定转换：

vsd <- vst(dds, blind=FALSE) expr_data <- assay(vsd)

2. 数据预处理与质量检查

良好的数据预处理是WGCNA成功的关键。我们需要确保输入数据的质量，并进行适当的过滤和转换。

数据过滤建议：

• 去除在超过90%样本中表达量极低的基因
• 保留表达变异较大的基因（如标准差排名前50%的基因）
• 检查样本聚类是否存在明显离群值

# 计算基因表达的标准差 gene_sd <- apply(expr_data, 2, sd) # 保留标准差大于中位数的基因 keep_genes <- gene_sd > median(gene_sd) expr_filtered <- expr_data[, keep_genes] # 样本聚类检查 sampleTree <- hclust(dist(expr_filtered), method = "average") plot(sampleTree, main = "Sample clustering", sub="", xlab="")

注意：如果样本聚类图中出现明显离群的样本，需要检查实验记录或技术因素，考虑是否移除该样本。

3. 构建基因共表达网络

WGCNA的核心是构建加权基因共表达网络。这一步的关键是选择合适的软阈值（soft threshold）β值，它决定了基因间相关性强度的转换方式。

软阈值选择步骤：

1. 测试一系列β值（通常1-20）
2. 评估网络是否符合无标度拓扑特征
3. 选择使网络达到无标度拓扑的最小β值

# 测试软阈值 powers <- c(1:20) sft <- pickSoftThreshold(expr_filtered, powerVector = powers, verbose = 5) # 绘制结果 par(mfrow = c(1,2)) plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)", ylab="Scale Free Topology Model Fit", type="n", main = paste("Scale independence")) text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers, col="red") abline(h=0.85, col="red") 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")

选择使无标度拓扑拟合指数（左图）首次达到0.85以上的最小β值。右图展示了不同β值下网络的平均连接度，β值增大通常会导致连接度下降。

4. 模块识别与可视化

确定了合适的软阈值后，我们可以进行模块识别。WGCNA使用动态树切割算法将基因聚类到不同的模块中。

# 一步法构建网络和识别模块 net <- blockwiseModules(expr_filtered, power = sft$powerEstimate, TOMType = "signed", minModuleSize = 30, reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE, pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = "geneTOM", verbose = 3) # 查看模块数量及大小 table(net$colors)

模块识别后，我们可以用热图展示模块特征基因（eigengene）之间的关系：

# 计算模块特征基因 MEs <- net$MEs # 计算模块特征基因间的相关性 ME_cor <- cor(MEs) # 绘制热图 heatmap(ME_cor, symm = TRUE, margins = c(10,10), col = colorRampPalette(c("blue", "white", "red"))(100))

5. 关联模块与表型数据

识别出的基因模块是否有生物学意义？我们需要将模块与样本的表型数据关联起来分析。假设我们有一个包含临床特征的数据框traitData。

# 计算模块特征基因与表型的相关性 moduleTraitCor <- cor(MEs, traitData, use = "p") moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(expr_filtered)) # 可视化 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(traitData), yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE, colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.5, zlim = c(-1,1), main = "Module-trait relationships")

热图中每个单元格显示了模块特征基因与表型的相关系数及显著性（括号内p值）。红色表示正相关，蓝色表示负相关。

6. 识别关键基因（Hub Genes）

在每个模块中，连接度最高的基因被称为hub基因，它们通常是模块中最重要的调控因子。我们可以计算基因的模块成员度（MM）和基因显著性（GS）来识别hub基因。

# 计算所有基因的模块成员度 geneModuleMembership <- as.data.frame(cor(expr_filtered, MEs, use = "p")) MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nrow(expr_filtered))) # 计算基因与表型的相关性（基因显著性） geneTraitSignificance <- as.data.frame(cor(expr_filtered, traitData, use = "p")) GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nrow(expr_filtered))) # 选择与目标表型最相关的模块 module <- "blue" # 替换为你感兴趣的模块名称 column <- match(module, names(MEs)) moduleGenes <- net$colors == module # 绘制模块成员度与基因显著性的关系 par(mfrow = c(1,1)) verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]), abs(geneTraitSignificance[moduleGenes, 1]), xlab = paste("Module Membership in", module, "module"), ylab = "Gene significance for target trait", main = paste("Module membership vs. gene significance\n"), cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

图中右上角的基因既是模块的核心成员（高MM值），又与目标表型高度相关（高GS值），是潜在的hub基因候选。

7. 结果解读与生物学意义挖掘

获得hub基因后，下一步是探索它们的生物学功能。我们可以使用以下方法：

1. 功能富集分析：对每个模块的基因进行GO或KEGG通路富集
2. 蛋白互作网络分析：将hub基因导入STRING等数据库构建蛋白互作网络
3. 文献挖掘：查询hub基因在相关疾病或表型中的已知功能

# 提取特定模块的基因列表 blue_module_genes <- colnames(expr_filtered)[net$colors == "blue"] write.table(blue_module_genes, file = "blue_module_genes.txt", quote = FALSE, row.names = FALSE, col.names = FALSE) # 提取hub基因（模块内连接度最高的基因） intramodular_connectivity <- intramodularConnectivity(adjacency, net$colors) hub_genes <- rownames(intramodular_connectivity)[order(intramodular_connectivity$kWithin, decreasing = TRUE)[1:10]]

8. 高级技巧与常见问题解决

内存优化技巧：

• 对于大型数据集，使用blockwiseModules的分块计算功能
• 设置maxBlockSize参数控制每块处理的基因数量
• 使用saveTOMFileBase参数保存中间结果避免重复计算

网络类型选择：

• signed：只保留正相关，区分强正相关和弱相关
• unsigned：考虑绝对相关性，不区分正负
• signed hybrid：正相关加权处理，负相关简单处理

常见问题处理：

1. 样本量不足：WGCNA推荐至少15-20个样本，样本量过少会导致网络不稳定
2. 批次效应：使用sva::ComBat等工具校正批次效应
3. 模块过多或过少：调整mergeCutHeight和minModuleSize参数
4. 无标度拓扑拟合不佳：尝试更高的软阈值或检查数据质量

# 分块计算示例 net <- blockwiseModules(expr_filtered, power = sft$powerEstimate, maxBlockSize = 5000, # 每块处理5000个基因 TOMType = "signed", minModuleSize = 30, saveTOMs = TRUE, saveTOMFileBase = "geneTOM", verbose = 3)

9. 完整流程整合与自动化脚本

为了便于重复分析，我们可以将整个流程整合到一个R脚本中。以下是一个简化的流程框架：

# WGCNA完整分析流程框架 run_wgcna_analysis <- function(expr_data, trait_data, output_dir) { # 1. 数据预处理 expr_filtered <- filter_low_expr_genes(expr_data) # 2. 软阈值选择 sft <- pick_soft_threshold(expr_filtered) # 3. 网络构建与模块识别 net <- build_coexp_network(expr_filtered, sft$powerEstimate) # 4. 模块-表型关联 module_trait_cor <- correlate_modules_traits(net$MEs, trait_data) # 5. Hub基因识别 hub_genes <- identify_hub_genes(expr_filtered, net, trait_data) # 6. 结果保存 save_results(net, module_trait_cor, hub_genes, output_dir) return(list(net=net, module_trait_cor=module_trait_cor, hub_genes=hub_genes)) }

实际项目中，你可能需要根据具体需求调整参数和分析步骤。WGCNA的强大之处在于其灵活性，可以根据不同的生物学问题和数据类型进行定制化分析。