1. WGCNA分析前的准备工作WGCNAWeighted Gene Co-expression Network Analysis是一种强大的生物信息学分析方法特别适合处理高通量基因表达数据。在开始分析之前我们需要做好充分的准备工作。首先确保你已经安装了最新版本的R语言环境。我推荐使用R 4.0或更高版本因为新版本对内存管理和计算效率都有显著提升。安装完R后我们需要安装WGCNA包及其依赖包。这里有个小技巧如果你使用的是R 3.5及以上版本可以使用BiocManager来安装相关包这样能自动解决依赖关系问题。install.packages(BiocManager) BiocManager::install(c(AnnotationDbi, impute, GO.db, preprocessCore)) install.packages(c(matrixStats, Hmisc, foreach, doParallel, fastcluster, dynamicTreeCut, survival, WGCNA, stringr, reshape2))安装完成后记得加载WGCNA包library(WGCNA)在实际项目中我经常遇到包版本冲突的问题。比如WGCNA包中的cor函数会覆盖stats包中的同名函数。为了避免混淆我习惯在调用特定函数时使用双冒号语法比如WGCNA::cor。2. 数据导入与清洗2.1 加载基因表达数据数据质量是WGCNA分析成功的关键。我们先从加载基因表达数据开始。假设我们有一个CSV格式的表达矩阵文件通常第一列是基因ID后续各列是不同样本的表达值。options(stringsAsFactors FALSE) # 避免自动将字符串转为因子 femData - read.csv(LiverFemale3600.csv)加载数据后我们需要检查数据结构。使用dim()查看数据维度head()预览前几行数据。在实际操作中我经常发现原始数据包含很多注释列这些列在后续分析中并不需要我们可以先将其去除。datExpr0 - as.data.frame(t(femData[, -c(1:8)])) # 删除前8列注释信息并转置 names(datExpr0) - femData$geneID # 设置列名为基因ID rownames(datExpr0) - colnames(femData)[-c(1:8)] # 设置行名为样本名2.2 数据质量检查与清洗数据清洗是WGCNA分析中最容易被忽视但至关重要的步骤。我们需要检查数据中的缺失值和异常值。gsg - goodSamplesGenes(datExpr0, verbose 3) if (!gsg$allOK) { # 打印并移除有问题的基因和样本 if (sum(!gsg$goodGenes)0) printFlush(paste(Removing genes:, paste(names(datExpr0)[!gsg$goodGenes], collapse , ))) if (sum(!gsg$goodSamples)0) printFlush(paste(Removing samples:, paste(rownames(datExpr0)[!gsg$goodSamples], collapse , ))) datExpr0 - datExpr0[gsg$goodSamples, gsg$goodGenes] }接下来我们需要检查样本是否存在异常。通过层次聚类可以帮助我们识别离群样本sampleTree - hclust(dist(datExpr0), method average) plot(sampleTree, main Sample clustering to detect outliers)如果发现明显偏离的样本比如在聚类树中单独成枝的样本可以考虑将其移除。在我的经验中这一步往往能发现实验过程中的批次效应或技术问题。2.3 加载临床特征数据表型数据是WGCNA分析的另一重要组成部分。我们需要确保表型数据与表达数据的样本顺序一致。traitData - read.csv(ClinicalTraits.csv) # 清理不需要的列 allTraits - traitData[, -c(31, 16)] allTraits - allTraits[, c(2, 11:36)] # 匹配样本 femaleSamples - rownames(datExpr0) traitRows - match(femaleSamples, allTraits$Mice) datTraits - allTraits[traitRows, -1] rownames(datTraits) - allTraits[traitRows, 1]为了直观展示样本聚类与表型特征的关系我们可以绘制热图sampleTree2 - hclust(dist(datExpr0), method average) traitColors - numbers2colors(datTraits, signed FALSE) plotDendroAndColors(sampleTree2, traitColors, groupLabels names(datTraits))3. 构建共表达网络与模块识别3.1 软阈值的选择软阈值的选择是WGCNA分析中最关键的步骤之一。它决定了基因间相关性强度的转换方式。我们需要选择一个合适的幂值power来构建无标度网络。powers - c(c(1:10), seq(from 12, to 20, by 2)) sft - pickSoftThreshold(datExpr0, powerVector powers, verbose 5)选择标准是使网络达到近似无标度拓扑结构scale-free topology。通常我们选择使拟合指数R²达到0.9以上的最小power值。plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab Soft Threshold (power), ylab Scale Free Topology Model Fit) abline(h 0.90, col red)在实际分析中我建议不仅要看R²值还要考虑平均连接度mean connectivity。过高的power值会导致网络过于稀疏丢失有价值的信息。3.2 一步法构建网络与识别模块确定好power值后我们可以一步完成网络构建和模块识别net - blockwiseModules(datExpr0, power sft$powerEstimate, TOMType unsigned, minModuleSize 30, reassignThreshold 0, mergeCutHeight 0.25, numericLabels TRUE, pamRespectsDendro FALSE, saveTOMs TRUE, saveTOMFileBase femaleMouseTOM, verbose 3)这个函数会输出检测到的模块数量和各模块包含的基因数。我们可以用以下代码查看结果table(net$colors)为了直观展示模块划分结果我们可以绘制基因树状图及模块颜色mergedColors - labels2colors(net$colors) plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]], Module colors, dendroLabels FALSE, hang 0.03, addGuide TRUE, guideHang 0.05)3.3 模块特征基因计算模块特征基因Module Eigengene, ME是每个模块中基因表达模式的第一主成分代表了该模块的整体表达模式。MEs - net$MEs geneTree - net$dendrograms[[1]] moduleLabels - net$colors moduleColors - labels2colors(net$colors)保存这些中间结果是个好习惯可以避免重复计算save(MEs, moduleLabels, moduleColors, geneTree, file FemaleLiver-networkConstruction-auto.RData)4. 模块与表型的关联分析4.1 量化模块-表型关联这是WGCNA分析的核心部分我们将探索哪些基因模块与关注的表型特征相关。moduleTraitCor - cor(MEs, datTraits, use p) moduleTraitPvalue - corPvalueStudent(moduleTraitCor, nrow(datExpr0))为了直观展示结果我们可以绘制热图textMatrix - paste(signif(moduleTraitCor, 2), \n(, signif(moduleTraitPvalue, 1), ), sep ) dim(textMatrix) - dim(moduleTraitCor) 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 Module-trait relationships)4.2 基因重要性与模块成员分析对于与表型显著相关的模块我们可以进一步分析模块内的基因。weight - as.data.frame(datTraits$weight_g) names(weight) - weight geneModuleMembership - as.data.frame(cor(datExpr0, MEs, use p)) MMPvalue - as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nrow(datExpr0))) modNames - substring(names(MEs), 3) names(geneModuleMembership) - paste(MM, modNames, sep ) names(MMPvalue) - paste(p.MM, modNames, sep ) geneTraitSignificance - as.data.frame(cor(datExpr0, weight, use p)) GSPvalue - as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nrow(datExpr0))) names(geneTraitSignificance) - paste(GS., names(weight), sep ) names(GSPvalue) - paste(p.GS., names(weight), sep )4.3 模块内关键基因的识别我们可以通过基因重要性Gene Significance, GS和模块成员Module Membership, MM来识别模块中的关键基因。module - brown # 选择感兴趣的模块 column - match(module, modNames) moduleGenes - moduleColors module verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]), abs(geneTraitSignificance[moduleGenes, 1]), xlab paste(Module Membership in, module, module), ylab Gene significance for body weight, main paste(Module membership vs. gene significance\n), cex.main 1.2, cex.lab 1.2, cex.axis 1.2, col module)这个散点图展示了模块成员与基因重要性之间的关系。通常我们会关注右上角的基因它们既与模块高度相关又与表型高度相关。5. 结果可视化与解读5.1 网络热图可视化虽然全基因网络的TOM热图计算量很大但能提供全局视角dissTOM - 1 - TOMsimilarityFromExpr(datExpr0, power 6) plotTOM - dissTOM^7 diag(plotTOM) - NA TOMplot(plotTOM, geneTree, moduleColors, main Network heatmap plot, all genes)对于大数据集我们可以随机选择部分基因来绘制热图nSelect - 400 set.seed(10) select - sample(ncol(datExpr0), size nSelect) selectTOM - dissTOM[select, select] selectTree - hclust(as.dist(selectTOM), method average) selectColors - moduleColors[select] plotDiss - selectTOM^7 diag(plotDiss) - NA TOMplot(plotDiss, selectTree, selectColors, main Network heatmap plot, selected genes)5.2 模块关系网络可视化我们可以通过特征基因网络来展示模块之间的关系MET - orderMEs(cbind(MEs, weight)) plotEigengeneNetworks(MET, , marDendro c(0,4,1,2), marHeatmap c(3,4,1,2))这个图分为两部分上部是特征基因的聚类树状图下部是特征基因间的相关性热图。它可以帮助我们识别元模块meta-modules即一组高度相关的模块。5.3 结果整理与输出最后我们需要将分析结果整理成表格方便后续研究annot - read.csv(file GeneAnnotation.csv) probes - names(datExpr0) probes2annot - match(probes, annot$geneID) geneInfo0 - data.frame(geneID probes, geneSymbol annot$gene_symbol[probes2annot], moduleColor moduleColors, geneTraitSignificance, GSPvalue) modOrder - order(-abs(cor(MEs, weight, use p))) for (mod in 1:ncol(geneModuleMembership)) { oldNames - names(geneInfo0) geneInfo0 - data.frame(geneInfo0, geneModuleMembership[, modOrder[mod]], MMPvalue[, modOrder[mod]]) names(geneInfo0) - c(oldNames, paste(MM., modNames[modOrder[mod]], sep ), paste(p.MM., modNames[modOrder[mod]], sep )) } geneOrder - order(geneInfo0$moduleColor, -abs(geneInfo0$GS.weight)) geneInfo - geneInfo0[geneOrder, ] write.csv(geneInfo, file geneInfo.csv)这个结果文件包含了每个基因的模块归属信息、与目标表型的相关性以及统计显著性是后续功能分析和实验验证的重要依据。