如果你正在处理高通量基因表达数据比如RNA-seq或芯片数据面对成千上万个基因你是否感到无从下手传统的差异表达分析只能告诉你哪些基因“变了”却无法揭示这些基因之间如何“协同工作”。这正是许多生物信息学新手在完成差异分析后遇到的第一个瓶颈知道了“谁在变”却不知道“他们为什么一起变”。这时WGCNAWeighted Gene Co-Expression Network Analysis加权基因共表达网络分析就成为了破局的关键。它不是一个简单的统计检验而是一套强大的系统生物学工具能将海量的基因表达数据转化为一张清晰的“社交网络图”帮你找到那些行为高度一致的基因模块Module并挖掘模块与关键性状如疾病分期、药物处理、表型数据之间的深层关联。网络上关于WGCNA的教程很多但往往陷入两个极端要么是充斥着R语言代码和数学公式的“天书”让初学者望而却步要么是过于简化的流程隐藏了关键的参数调整和结果解读细节导致你即使跑通了代码也对结果一知半解。本文的目的就是带你穿透迷雾实现“一个视频学会”的深度与效率。我们不会停留在代码的简单罗列而是会深入解释为什么WGCNA能解决共表达问题每一步关键参数如软阈值power的选择背后有何玄机如何从生成的网络中提取出有生物学意义的结论更重要的是我们会手把手带你走完从数据预处理、网络构建、模块识别到性状关联的全流程并提供可复现的完整代码和避坑指南。无论你是生物信息学的入门者还是希望系统掌握WGCNA的进阶者这篇文章都将是你案头必备的实战手册。1. WGCNA要解决的核心问题从“差异”到“关联”在深入代码之前我们必须先理解WGCNA究竟要做什么。这决定了我们后续所有分析步骤的目标。假设你有一项研究比较健康组和癌症组的基因表达谱。差异表达分析如DESeq2, edgeR会给你一份长长的“差异基因列表”DEGs。这份列表很重要但它存在几个局限信息丢失成百上千个DEGs难以形成整体认知。功能孤立列表中的基因看似独立但它们可能隶属于同一个信号通路或生物学过程。驱动因素不明哪些基因是核心调控者哪些模块与癌症恶性程度最相关WGCNA的智慧在于它转换了视角不关注单个基因的差异而是关注基因与基因之间表达模式的相似性共表达。表达模式高度相似的基因更可能被共同调控参与相同的生物学功能。因此WGCNA分析的核心产出是基因共表达网络一个以基因作为节点以基因间表达相似性加权值作为连接边权重的网络。基因模块通过网络聚类将基因划分为若干个内部高度连通、外部连接稀疏的“社区”每个社区就是一个模块如“蓝色模块”、“棕色模块”。模块-性状关联计算每个模块的“代表值”模块特征基因ME并与样本的性状数据临床信息进行关联分析找出与目标性状最相关的模块。核心基因在每个关键模块内部通过计算连通性等指标找出位于网络中心位置的“枢纽基因”Hub Genes它们往往是调控的关键。所以当你启动WGCNA分析时你真正要解决的问题是在我的数据中哪些基因倾向于“抱团”行动这些“团体”中哪一个与我最关心的疾病表型或实验处理关联最强这个“团体”里的核心成员枢纽基因是谁2. 核心概念与原理为什么是“加权”网络理解几个核心概念是避免后续操作沦为“黑箱”的关键。2.1 共表达相似性与邻接矩阵WGCNA的基础是计算任意两个基因在所有样本中表达量的相关性。通常使用斯皮尔曼Spearman或皮尔逊Pearson相关系数。假设有基因i和基因j它们的相关系数为 ( s_{ij} )范围在-1到1之间。然而直接使用相关系数构建网络即硬阈值法存在弊端你需要设定一个阈值如 |r| 0.8高于阈值的连接设为1低于的设为0。这种方法武断地切断了弱连接而生物学中许多有意义的调控关系可能是中低强度的。2.2 软阈值与加权网络WGCNA的“加权”Weighted精髓就在于引入了软阈值Soft Thresholding。它通过一个幂函数将相关系数转化为邻接矩阵Adjacency Matrix的权重 [ a_{ij} |s_{ij}|^\beta ] 其中( \beta ) 就是软阈值功率soft thresholding power。这个变换有两大好处保留弱连接即使相关系数不高经过幂运算后仍会有一个小权重而不是被直接归零。强化无尺度拓扑通过选择合适的 ( \beta )可以使最终网络的连接度分布近似服从无尺度分布即大部分节点连接少少数枢纽节点连接极多。这在生物学网络中很常见。如何选择 ( \beta ) 这是WGCNA第一个关键步骤。通常我们会绘制不同 ( \beta ) 值下网络的“无尺度拓扑拟合指数”scale-free topology fit index图选择使该指数达到较高水平如 0.8的最小 ( \beta ) 值。同时还要兼顾平均连接度不能太低以保证网络信息量。2.3 拓扑重叠矩阵与模块识别直接使用邻接矩阵进行聚类仍会受噪声影响。WGCNA进一步计算了拓扑重叠矩阵Topological Overlap Matrix, TOM。TOM不仅考虑两个基因是否直接相关还考虑它们是否共享相似的邻居。这能更稳健地衡量基因在网络中的功能相似性。基于TOM距离1-TOM采用层次聚类Hierarchical Clustering和动态树切割Dynamic Tree Cut方法将基因划分成不同的模块。动态树切割算法能智能地根据树枝形状确定切割高度比固定高度切割更优。2.4 模块特征基因与性状关联每个模块可以用一个“代表”来概括其表达模式即模块特征基因Module Eigengene, ME。ME本质上是该模块基因表达矩阵的第一主成分。然后计算ME与样本性状如疾病评分、生存时间、处理组别的相关性得到模块-性状关联热图直观展示哪些模块与目标性状最相关。2.5 基因显著性、模块成员与枢纽基因基因显著性Gene Significance, GS单个基因与目标性状的相关性绝对值。GS越高说明该基因与性状关联越强。模块成员Module Membership, MM也称为kME指单个基因表达谱与其所在模块ME的相关性。MM越高说明该基因在其模块内的“代表性”越强。枢纽基因Hub Gene通常指在一个模块内同时具有高GS与性状相关和高MM在模块内核心的基因。它们是后续实验验证的首要候选。3. 环境准备与数据要求3.1 R与RStudio环境WGCNA是一个R语言包因此你需要安装最新版的 R 。安装 RStudio 推荐非必须。安装必要的R包。在R中运行以下命令# 设置CRAN镜像加速下载以清华镜像为例 options(repos c(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN/)) # 安装WGCNA核心包及其依赖这可能需要几分钟 if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(c(WGCNA, impute, preprocessCore)) install.packages(c(matrixStats, Hmisc, foreach, doParallel, fastcluster, dynamicTreeCut, survival)) # 安装常用的辅助包用于数据处理和可视化 install.packages(c(tidyverse, ggplot2, pheatmap, corrplot, reshape2))3.2 数据要求与格式WGCNA对输入数据有明确要求准备不当是失败的主要原因。表达矩阵Expression Data格式一个数据框data.frame或矩阵matrix行是基因列是样本。内容通常是标准化后的表达量如FPKM、TPMRNA-seq或标准化后的信号强度芯片数据。强烈建议使用方差稳定变换如log2(TPM1)或标准化后的数据。清洗去除低表达或缺失值过多的基因。例如保留在至少80%的样本中表达量大于1的基因。样本量WGCNA需要一定样本量来稳定估计相关性。通常建议至少15-20个样本越多越好。性状数据Trait Data格式一个数据框行是样本与表达矩阵列名顺序一致列是性状。内容可以是连续型如肿瘤大小、血压值、二分类0/1代表健康/患病或有序分类数据。需要预先处理好。一个常见的数据结构示例假设你有20个样本测量了10000个基因的表达量。exp_data: 一个10000行 x 20列的矩阵。trait_data: 一个20行 x 3列的数据框三列分别是SampleID、DiseaseStageI, II, III、SurvivalTime。4. 完整WGCNA分析流程拆解含代码下面我们用一个模拟的、但结构真实的数据集走通整个WGCNA流程。你可以将代码中的路径和数据替换成你自己的。4.1 步骤一加载包与数据预处理# 文件路径WGCNA_Analysis.R # 1. 加载必要的包 library(WGCNA) library(tidyverse) library(pheatmap) options(stringsAsFactors FALSE) # 避免字符自动转因子 enableWGCNAThreads() # 启用多线程加速如果支持 # 2. 模拟加载表达数据实际中请替换为你的数据读取代码 # 假设我们有一个名为“gene_expression_matrix.csv”的文件 # exp_data - read.csv(gene_expression_matrix.csv, row.names 1) # 这里我们创建一个模拟数据以便演示 set.seed(123) nGenes - 2000 nSamples - 30 exp_data - matrix(rnorm(nGenes * nSamples, mean10, sd2), nrownGenes, ncolnSamples) rownames(exp_data) - paste0(Gene, 1:nGenes) colnames(exp_data) - paste0(Sample, 1:nSamples) # 模拟一些基因具有共表达模式 exp_data[1:100, ] - exp_data[1:100, ] rnorm(100*nSamples, mean0, sd0.5) # 模块1 exp_data[101:200, ] - exp_data[101:200, ] rnorm(100*nSamples, mean1, sd0.5) # 模块2 # 3. 数据清洗过滤低表达基因 gsg - goodSamplesGenes(exp_data, verbose 3) gsg$allOK # 如果为TRUE则无需过滤如果为FALSE需要移除不符合条件的基因和样本 if (!gsg$allOK){ # 移除不符合条件的基因 exp_data - exp_data[gsg$goodGenes, gsg$goodSamples] print(paste(Removed, sum(!gsg$goodGenes), genes and, sum(!gsg$goodSamples), samples.)) } # 4. 加载性状数据实际中请替换 # trait_data - read.csv(trait_data.csv, row.names 1) # 创建模拟性状数据 trait_data - data.frame( SampleID colnames(exp_data), DiseaseStage sample(c(StageI, StageII, StageIII), nSamples, replace TRUE), TumorSize rnorm(nSamples, mean5, sd1.5), Response sample(c(CR, PR, SD), nSamples, replace TRUE) ) rownames(trait_data) - trait_data$SampleID trait_data$SampleID - NULL # 将分类性状转换为数值WGCNA关联分析需要数值 trait_data_numeric - model.matrix(~0., datatrait_data) colnames(trait_data_numeric) - gsub(^DiseaseStage|^Response, , colnames(trait_data_numeric)) # 检查样本顺序是否一致 print(all(rownames(trait_data_numeric) colnames(exp_data)))4.2 步骤二样本聚类与异常值检测在构建网络前检查样本是否有明显异常这会影响相关性计算。# 基于表达数据对样本进行聚类 sampleTree - hclust(dist(t(exp_data)), method average) # 绘制样本聚类树 par(cex 0.6) par(mar c(0,4,2,0)) plot(sampleTree, main Sample clustering to detect outliers, sub, xlab, cex.lab 1.5, cex.axis 1.5, cex.main 2) # 如果发现明显远离其他样本的离群点可以考虑手动移除 # 例如假设我们想剪掉高度200的枝 # cutHeight - 200 # clust - cutreeStatic(sampleTree, cutHeight cutHeight, minSize 10) # keepSamples - (clust1) # exp_data - exp_data[, keepSamples] # trait_data_numeric - trait_data_numeric[keepSamples, ]4.3 步骤三选择软阈值功率β这是最关键的一步直接影响网络性质。# 定义一组候选的软阈值功率 powers - c(1:10, seq(12, 20, by2)) # 调用函数进行网络拓扑分析 sft - pickSoftThreshold(exp_data, powerVector powers, verbose 5, networkType unsigned) # 网络类型通常选择 unsigned (只考虑相关性绝对值), signed (区分正负相关), signed hybrid # 绘制结果图 par(mfrow c(1,2)) cex1 - 0.9 # 图1无尺度拓扑拟合指数 vs. 软阈值 plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlabSoft Threshold (power), ylabScale Free Topology Model Fit, signed R^2, typen, main paste(Scale independence)) text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], labelspowers, cexcex1, colred) abline(h0.85, colred) # 通常以0.85作为参考线 # 图2平均连接度 vs. 软阈值 plot(sft$fitIndices[,1], sft$fitIndices[,5], xlabSoft Threshold (power), ylabMean Connectivity, typen, main paste(Mean connectivity)) text(sft$fitIndices[,1], sft$fitIndices[,5], labelspowers, cexcex1, colred) # 查看sft$fitIndices表格选择使R^2 0.85且平均连接度不过低通常10的最小power print(sft$fitIndices) # 假设我们根据图表选择 power 6 softPower - 64.4 步骤四一步法构建网络与识别模块WGCNA提供了blockwiseModules函数可以高效地一次性完成网络构建和模块识别特别适合基因数很多5000的情况。对于中小规模数据也可以用blockwiseModules。# 设置网络构建参数 net - blockwiseModules(exp_data, power softPower, # 上一步选择的软阈值 TOMType unsigned, # 与pickSoftThreshold时一致 minModuleSize 30, # 最小模块基因数可根据数据调整 deepSplit 2, # 控制切割灵敏度0-4越大模块越多越小 pamRespectsDendro FALSE, # 通常设为FALSE mergeCutHeight 0.25, # 合并相似模块的阈值越小合并越少 numericLabels TRUE, # 模块用数字0,1,2...标记0通常代表未分组的基因 saveTOMs TRUE, # 保存TOM矩阵供后续分析 saveTOMFileBase MyNetworkTOM, # TOM文件前缀 verbose 3) # 查看模块数量及大小 table(net$colors) # net$colors 是一个向量长度等于基因数每个基因被分配了一个模块颜色数字代码4.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. 查看模块特征基因MEs MEs - net$MEs # MEs是一个数据框行是样本列是各模块的ME如ME1, ME2... head(MEs)4.6 步骤六关联模块与外部性状这是将网络与生物学意义连接起来的关键步骤。# 1. 计算模块特征基因与性状的相关性及p值 moduleTraitCor - cor(MEs, trait_data_numeric, use p) moduleTraitPvalue - corPvalueStudent(moduleTraitCor, nSamples) # 2. 可视化模块-性状关联热图 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 colnames(trait_data_numeric), yLabels names(MEs), ySymbols names(MEs), colorLabels FALSE, colors blueWhiteRed(50), textMatrix textMatrix, setStdMargins FALSE, cex.text 0.7, zlim c(-1,1), main paste(Module-trait relationships)) # 热图中每个单元格显示相关系数和p值。颜色越深红表示正相关越强越深蓝表示负相关越强。 # 例如你可能发现“MEblue”模块与“TumorSize”强正相关r0.8, p0.01。4.7 步骤七深入关键模块寻找枢纽基因假设我们发现“蓝色模块”对应数字标签比如是ME1与“肿瘤大小”最相关。# 1. 定义我们感兴趣的性状例如“TumorSize” trait_of_interest - as.data.frame(trait_data_numeric[,TumorSize]) colnames(trait_of_interest) - TumorSize # 2. 计算基因显著性GS GS - as.numeric(cor(exp_data, trait_of_interest, use p)) # 给GS命名 names(GS) - rownames(exp_data) # 3. 计算模块成员MM即基因与模块特征基因的相关性 # 首先确定蓝色模块对应的ME列名。假设ME1是蓝色模块。 module - blue # 模块颜色 module_genes - (moduleColors module) # 逻辑向量标记属于该模块的基因 modME - MEs[, paste0(ME, module)] # 获取该模块的ME # 计算这些基因的MM MM - as.numeric(cor(exp_data[, module_genes], modME, use p)) names(MM) - rownames(exp_data)[module_genes] # 4. 绘制GS vs MM 散点图查看相关性 par(mfrowc(1,1)) verboseScatterplot(MM, GS[module_genes], xlab paste(Module Membership in, module, module), ylab paste(Gene significance for, colnames(trait_of_interest)), main paste(MM vs. GS\n), cex.main 1.2, cex.lab 1.2, cex.axis 1.2, col module) # 通常GS和MM呈正相关说明该模块中与性状相关的基因也往往是模块的核心成员。 # 5. 提取枢纽基因可以定义同时满足GS和MM高阈值的基因 # 例如选择GS 0.5 且 MM 0.8 的基因 hub_genes - names(GS[module_genes])[abs(GS[module_genes]) 0.5 abs(MM) 0.8] print(paste(Number of hub genes in, module, module:, length(hub_genes))) print(head(hub_genes))4.8 步骤八结果导出与可视化将关键结果导出用于后续分析和绘图。# 1. 将基因与模块颜色、GS等信息整合到一个数据框 geneInfo - data.frame( GeneID rownames(exp_data), ModuleColor moduleColors, GeneSignificanceTumorSize GS ) # 可以进一步合并每个基因在其所属模块内的MM # 这需要循环计算每个模块的MM这里提供一个简化示例 # 实际中可以使用WGCNA的signedKME函数批量计算 kME - signedKME(exp_data, MEs) colnames(kME) - gsub(kME, MM., colnames(kME)) geneInfo - cbind(geneInfo, kME) # 2. 保存结果 write.csv(geneInfo, file WGCNA_Gene_Module_Info.csv, row.names FALSE) write.csv(moduleTraitCor, file Module_Trait_Correlation.csv) write.csv(moduleTraitPvalue, file Module_Trait_Pvalue.csv) # 3. 可视化模块内基因的互作网络以TOM为基础展示部分核心基因 # 选择蓝色模块的TOM子矩阵 module - blue module_genes - (moduleColors module) modTOM - TOM[module_genes, module_genes] # 注意需要从保存的TOM文件中加载或重新计算 dimnames(modTOM) - list(rownames(exp_data)[module_genes], rownames(exp_data)[module_genes]) # 导出为Cytoscape等网络可视化软件可读的格式 cyt - exportNetworkToCytoscape(modTOM, edgeFile paste(CytoscapeEdge-, module, .txt, sep), nodeFile paste(CytoscapeNode-, module, .txt, sep), weighted TRUE, threshold 0.02, # 只导出权重高于阈值的边控制网络规模 nodeNames rownames(exp_data)[module_genes], nodeAttr moduleColors[module_genes])5. 运行结果解读与验证运行完上述代码你会得到一系列文件和图表。如何判断分析是否成功并解读关键结果软阈值选择图检查选择的power值是否使左图的“Scale independence”指标R^2达到一个较高的平台如0.85同时右图的平均连接度没有骤降到极低。如果R^2始终很低可能数据本身不适合构建无尺度网络可考虑降低标准或使用“signed”网络类型。模块聚类树状图这是最直观的结果。观察树状图下方的颜色条成功的分析应该显示出几个颜色分明、边界清晰的色块模块。灰色或数字0通常代表未被分配到任何模块的基因。模块数量不宜过多如30或过少如5可通过调整minModuleSize和mergeCutHeight参数控制。模块-性状关联热图这是核心生物学发现的来源。重点关注高相关系数绝对值大例如MEblue与TumorSize的相关系数为0.82p0.001这意味着蓝色模块的整体表达模式与肿瘤大小高度正相关。模块内基因可能参与促进肿瘤生长的通路。高显著性p值小p值通常标注在括号内。p0.05表示关联显著但经过多重检验校正如FDR后更可靠。模式识别有时一个模块可能与多个性状相关这提示该模块可能处于核心调控地位。基因信息表geneInfo这张表是你的“基因花名册”。通过筛选ModuleColor和GeneSignificance你可以快速定位到关键模块中与性状最相关的基因列表用于后续的GO/KEGG富集分析、生存分析或实验验证。枢纽基因从hub_genes列表中获得的基因是后续功能研究和生物标志物开发的首选目标。建议在String数据库https://string-db.org/中检查这些基因的已知蛋白互作关系验证其网络中心性。6. 常见问题与排查思路WGCNA分析流程长参数多新手极易出错。下表总结了最常见的问题及解决方法问题现象可能原因排查方式解决方案pickSoftThreshold报错或R^2始终很低 (0.8)1. 数据未标准化或存在极端值。2. 样本量太少。3. 基因表达量变化太小噪音大。4. 数据本身就不符合无尺度网络特征。1. 检查表达矩阵summary(exp_data)看分布。2. 绘制样本聚类图看是否有异常样本。3. 尝试对数据做log2(x1)变换。1. 严格进行数据预处理标准化、去批次、过滤低表达基因。2. 增加样本量如果可能。3. 尝试使用networkType signed或signed hybrid。4. 如果R^2在0.7-0.8之间也可酌情继续或参考平均连接度选择power。模块数量过多30或过少3minModuleSize设置过小或过大deepSplit和mergeCutHeight参数设置不当。查看table(net$colors)输出观察模块大小分布。1. 调整minModuleSize常用30-100。2. 调整deepSplit0-4越大切割越细。3. 调整mergeCutHeight0.1-0.3越小越不易合并。建议多次尝试选择生物学上合理的模块数量。模块-性状关联全部不显著p值很大1. 性状数据与表达数据不匹配。2. 性状本身在样本间变异很小。3. 模块划分未能捕捉到与性状相关的共表达模式。1. 检查样本顺序all(rownames(trait_data) colnames(exp_data))。2. 检查性状数据的分布summary(trait_data)。3. 检查关键模块的ME是否在性状组间有差异箱线图。1. 确保性状与表达矩阵样本严格对应。2. 重新审视实验设计所选性状是否合理。3. 尝试不同的网络构建参数如power或使用其他聚类方法。运行blockwiseModules时内存不足或极慢基因数量太多20000一次性计算TOM矩阵内存消耗巨大。监控任务管理器中的内存使用。1. 使用blockwiseModules的分块计算功能默认已启用。2. 在函数中设置maxBlockSize如10000来分块。3. 预先过滤掉低方差或低表达的基因减少基因数。4. 使用高性能服务器或云计算资源。枢纽基因列表为空或很少GS和MM的筛选阈值设置过高。绘制verboseScatterplot查看GS和MM的分布情况。1. 降低GS和MM的筛选阈值如从 0.8 降到 0.6。2. 可以分别按GS和MM排序取前N个基因作为候选。导入Cytoscape的边文件过大软件卡死threshold设置过低导致导出的边数量过多成千上万。查看导出的*Edge.txt文件行数。1. 提高exportNetworkToCytoscape中的threshold参数如从0.02提高到0.1或0.15。2. 仅导出枢纽基因之间的连接子网络。7. 最佳实践与高级技巧掌握了基础流程后以下几点能让你的WGCNA分析更稳健、更具洞察力数据预处理是王道WGCNA对输入数据质量非常敏感。务必做好标准化使用limma、DESeq2的vst或rlog变换处理RNA-seq数据。去批次效应如果数据来自不同批次使用ComBatsva包等方法校正。过滤基因过滤在大部分样本中低表达或零表达的基因。goodSamplesGenes是基础检查但更严格的过滤如按方差能提升信噪比。参数选择不是一蹴而就power、minModuleSize、deepSplit、mergeCutHeight共同决定了模块的形态。没有“黄金参数”。建议进行参数敏感性分析固定其他参数微调其中一个观察模块数量、大小和性状关联稳定性的变化。选择结果稳健、生物学解释性强的参数组合。利用blockwiseModules的保存与加载功能网络构建耗时很长。使用saveTOMFileBase参数保存TOM矩阵。下次想用不同参数切割模块时可以使用recutBlockwiseTrees函数直接加载已保存的TOM快速重新聚类无需重复计算网络。模块的功能注释至关重要得到关键模块后下一步就是将基因列表进行功能富集分析GO、KEGG。推荐使用clusterProfilerR包。这能直接回答“这个与性状相关的模块主要参与什么生物学过程”。与差异表达分析结合WGCNA和差异表达分析DEG是互补的。可以取交集找出既是DEG又位于关键模块中的基因。这些基因很可能扮演着更重要的角色。时间序列或复杂性状分析WGCNA可以处理时间序列数据将时间点作为性状也可以分析复杂性状如生存数据需要用到WGCNA的corPvalueStudent等函数进行适配。这需要更深入的学习。结果的可视化与报告除了内置函数使用ggplot2定制更精美的图表如模块特征基因表达趋势图、枢纽基因表达热图等能让你的文章或报告更加出彩。通过本文近万字的梳理你应该已经对WGCNA从原理、实战到排错有了系统的认识。记住WGCNA是一个强大的“探索性”工具它的价值在于从复杂数据中生成可验证的假设。真正的生物学结论还需要后续的实验进行验证。现在就打开RStudio载入你的数据开始构建你的第一个基因共表达网络吧。建议将本文代码保存为脚本并结合你的数据边运行边理解这才是“一个视频学会”背后的真谛——在动手实践中融会贯通。