从Seurat聚类到功能模块：手把手教你用hdWGCNA挖掘单细胞数据中的基因“朋友圈”

从Seurat聚类到功能模块：手把手教你用hdWGCNA挖掘单细胞数据中的基因“朋友圈”

单细胞转录组技术让我们能够以前所未有的分辨率观察细胞间的异质性，但仅仅知道"有哪些细胞类型"往往不能满足研究需求。当我们在Seurat中完成细胞聚类和注释后，一个更深入的问题自然浮现：在特定细胞亚群内部，哪些基因正在协同工作，共同执行特定的生物学功能？这正是hdWGCNA（high-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使用哪些基因进行分析。有三种主要选择策略：

# 使用在至少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, ncol=2)

理想的软阈值通常满足以下标准：

• 无标度拓扑拟合指数(R²) > 0.8
• 平均连接度适中(通常10-100)
• 选择第一个达到平台期的阈值点

3.3 构建共表达网络

选定软阈值后，我们可以构建基因共表达网络并识别模块：

seurat_obj <- ConstructNetwork( seurat_obj, tom_name = 'INH', soft_power = 8 # 根据测试结果选择 ) # 可视化基因模块 PlotDendrogram(seurat_obj, main='INH 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)

这种分析能揭示调控网络如何随时间重组，为发育机制提供见解。