从Seurat聚类到功能模块手把手教你用hdWGCNA挖掘单细胞数据中的基因“朋友圈”单细胞转录组技术让我们能够以前所未有的分辨率观察细胞间的异质性但仅仅知道有哪些细胞类型往往不能满足研究需求。当我们在Seurat中完成细胞聚类和注释后一个更深入的问题自然浮现在特定细胞亚群内部哪些基因正在协同工作共同执行特定的生物学功能这正是hdWGCNAhigh-dimensional Weighted Gene Co-expression Network Analysis大显身手的舞台。想象一下你刚完成了一个大脑皮层单细胞数据集的分析通过FindClusters和cell_type注释已经识别出了兴奋性神经元(EX)、抑制性神经元(INH)等主要细胞类型。现在你想深入了解INH神经元内部是否存在功能亚群以及这些亚群如何通过基因共表达网络调控特定神经功能。传统差异表达分析只能告诉你哪些基因在不同组间表达量有差异而hdWGCNA则能揭示基因间复杂的协作关系帮你找到那些志同道合的基因朋友圈。1. 准备工作与环境配置在开始hdWGCNA分析前确保你已经完成了标准的单细胞分析流程。以下是一个典型的准备工作清单数据预处理使用Seurat完成了NormalizeData、FindVariableFeatures、ScaleData等基本步骤降维与聚类执行了RunPCA/RunHarmony降维和FindClusters细胞聚类细胞类型注释基于已知标记基因完成了cell_type注释R环境准备安装必要的R包并加载# 安装必要R包如果尚未安装 if (!require(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(c(WGCNA, hdWGCNA)) # 加载所需库 library(Seurat) library(WGCNA) library(hdWGCNA) library(tidyverse) library(cowplot) library(patchwork) # 设置随机种子保证结果可重复 set.seed(12345) # 启用多线程加速计算 enableWGCNAThreads(nThreads 8)提示hdWGCNA对计算资源要求较高建议在服务器或高性能工作站上运行大型数据集。对于人脑皮层这类包含数万个细胞的数据集8-16GB内存是基本要求。2. 从Seurat对象到元细胞构建单细胞数据的一个主要挑战是数据稀疏性——单个细胞中许多基因的检测计数为零。hdWGCNA通过构建元细胞(Metacells)来缓解这一问题将相似的细胞聚合成更具代表性的表达谱。2.1 初始化hdWGCNA设置首先我们需要告诉hdWGCNA使用哪些基因进行分析。有三种主要选择策略基因选择方法适用场景参数设置fraction大多数情况选择在至少X%细胞中表达的基因variable关注高变基因使用Seurat已识别的高变基因custom特定研究需求提供自定义基因列表# 使用在至少5%细胞中表达的基因 seurat_obj - SetupForWGCNA( seurat_obj, gene_select fraction, fraction 0.05, wgcna_name tutorial )2.2 构建元细胞元细胞的构建需要考虑两个关键因素细胞类型(cell_type)和样本来源(sample)。这种双重分组确保了元细胞既反映生物学差异又考虑技术批次效应。seurat_obj - MetacellsByGroups( seurat_obj seurat_obj, group.by c(cell_type, sample), reduction harmony, k 25, # 每个元细胞包含约25个原始细胞 max_shared 10, ident.group cell_type ) # 标准化元细胞表达矩阵 seurat_obj - NormalizeMetacells(seurat_obj)注意参数k控制每个元细胞包含的细胞数量需要根据数据集大小调整。较大的k值会增加表达信号的稳定性但可能掩盖稀有细胞亚群的特征。3. 构建基因共表达网络有了元细胞表达矩阵后我们就可以开始探索基因间的共表达关系了。这个过程类似于社交网络分析——寻找那些表达模式高度相关的基因好友。3.1 设置目标细胞群体假设我们特别关注抑制性神经元(INH)内部的基因调控网络seurat_obj - SetDatExpr( seurat_obj, group_name INH, group.by cell_type, assay RNA, layer data )3.2 选择软阈值软阈值的选择是WGCNA分析中最关键的步骤之一它决定了基因间相关性强度的转换方式。太低的阈值会保留过多噪声而太高的阈值可能导致信息丢失。# 测试不同软阈值的影响 seurat_obj - TestSoftPowers( seurat_obj, networkType signed # 考虑相关性的方向 ) # 可视化结果 plot_list - PlotSoftPowers(seurat_obj) wrap_plots(plot_list, ncol2)理想的软阈值通常满足以下标准无标度拓扑拟合指数(R²) 0.8平均连接度适中(通常10-100)选择第一个达到平台期的阈值点3.3 构建共表达网络选定软阈值后我们可以构建基因共表达网络并识别模块seurat_obj - ConstructNetwork( seurat_obj, tom_name INH, soft_power 8 # 根据测试结果选择 ) # 可视化基因模块 PlotDendrogram(seurat_obj, mainINH hdWGCNA Dendrogram)在生成的树状图中不同颜色代表不同的基因模块。灰色模块包含未被归类的基因通常在下游分析中排除。4. 模块分析与生物学解释识别出基因模块后我们需要回答两个关键问题(1)这些模块代表什么生物学功能(2)它们与细胞表型有何关联4.1 模块特征基因与功能富集每个模块可以用其特征基因(eigengene)来代表这是模块内基因表达的第一主成分。我们可以计算模块与外部特征的关联# 计算模块特征基因 seurat_obj - ModuleEigengenes(seurat_obj) # 进行GO/KEGG功能富集分析 enrich_results - RunEnrichment( seurat_obj, db GO # 也可选择KEGG或其他数据库 ) # 可视化富集结果 PlotEnrichment(enrich_results, top_terms 5)典型的抑制性神经元模块可能富集在GABA合成与释放突触传递调控神经元发育与分化4.2 模块-性状关联分析如果我们有额外的表型数据(如疾病状态、治疗反应等)可以探索模块与这些性状的关联# 假设metadata中有疾病状态列 module_trait_cor - ModuleTraitCorrelation( seurat_obj, traits disease_status, cor_method pearson ) # 热图展示关联结果 PlotModuleTraitCorrelation(module_trait_cor)强相关的模块可能指示了与疾病相关的关键调控网络为机制研究提供线索。5. 高级应用与技巧掌握了基本流程后下面这些技巧可以帮助你从hdWGCNA中获得更多洞见5.1 跨细胞类型的比较网络分析# 同时分析INH和EX神经元 seurat_obj - SetDatExpr( seurat_obj, group_name c(INH, EX), group.by cell_type ) # 构建比较网络 comparative_net - CompareNetworks( seurat_obj, group.1 INH, group.2 EX ) # 可视化共享和特有模块 PlotComparativeNetwork(comparative_net)这种分析可以揭示不同神经元类型间保守和特异的调控模式。5.2 网络可视化与关键基因识别使用Cytoscape等工具对重要模块进行网络可视化# 提取特定模块的基因连接信息 module_genes - GetModuleGenes(seurat_obj, module blue) # 计算基因连接度 connectivity - ModuleConnectivity(seurat_obj, module blue) # 保存为Cytoscape可读格式 ExportNetworkToCytoscape( seurat_obj, module blue, top_genes 50, file INH_blue_module.txt )在网络中高度连接的hub基因往往是调控关键值得特别关注。5.3 时间序列数据的动态网络分析对于发育或时间序列数据可以追踪模块的动态变化# 按发育阶段分组构建元细胞 seurat_obj - MetacellsByGroups( seurat_obj, group.by c(cell_type, developmental_stage) ) # 分阶段构建网络 stage_networks - DynamicNetworkAnalysis( seurat_obj, time_var developmental_stage ) # 可视化模块演变 PlotDynamicModules(stage_networks)这种分析能揭示调控网络如何随时间重组为发育机制提供见解。