单细胞数据跑WGCNA这件事，我一开始是相当抗拒的。WGCNA当年是按bulk转录组设计的，输入矩阵是n个样本×m个基因，n一般只有几十、最多上百个；单细胞数据动辄几万个细胞，要是真把每个细胞当一个样本，内存先撑不住，统计假设更撑不住。更麻烦的是单细胞表达矩阵里60%到90%的位点都是零值，dropout效应把基因间的相关系数往0那边稀释，跑出来的共表达网络不是糊成一团就是被噪声带偏。后来有人把细胞按类型平均成伪bulk再跑，虽然算是能出结果了，但组内异质性被抹得干干净净，那些真正由亚群驱动的调控模块完全看不见。
hdWGCNA就是冲着这个矛盾来的。它不放弃单细胞分辨率，也不假装每个细胞都是独立样本，而是先把高度相似的细胞聚合成元细胞（meta-cell），再用元细胞当样本去做加权基因共表达网络分析。单细胞转录组和空间转录组都适用。这篇文章我会从原理讲到实操，把meta-cell是怎么构建的、软阈值为什么在单细胞里调得更高、hub基因怎么筛才靠谱、空间转录组又有哪些不一样的坑，一次性说清楚。
1. 为什么传统WGCNA撒到单细胞数据上就崩了
1.1 样本自由度是虚高的，显著性是骗人的
WGCNA的核心计算对象是基因与基因的相似度矩阵，理论上基因数量固定在两万上下时，不管样本是30个还是30000个，内存压力差别其实没有想象中大。真正出问题的是统计推断这一层。WGCNA算基因间相关性时，每个相关系数都对应一个p值，而p值取决于“样本量”。在bulk数据里样本量是几十个，基因相关性的检验本身已经很激进；在单细胞数据里样本量直接变成几万个细胞，自由度膨胀到上千，结果就是随便两个基因的p值都会被压缩到10的负几十次方，显著性筛选这个环节直接失效。
这里有个更隐蔽的问题：细胞压根不是独立抽样单位。测序拿到的那几万个细胞，来自同一个组织块、同一批建库，很多细胞之间本来就有极高的转录相似性。你可以把它们想象成对同一个人反复测量了十次血压，然后再拿这十次数据去做“十个人”的统计，自由度完全是被自己骗出来的。WGCNA的模块显著性检验、hub基因筛选全都要依赖这些统计量，基础不成立，后面全是空中楼阁。
1.2 dropout零膨胀把相关性稀释到失效
单细胞测序的捕获效率有限，一个基因明明在细胞里表达了，但反转录和扩增环节没有捕获到，测出来就是0。这种技术性缺失造成的零值占单细胞矩阵的绝大多数。比如一个基因平均表达量很低，它在一万个细胞里只有20%的检测率，那么它与另一个基因的共表达信号，很大程度会被两个基因同时“检测不到”的样本主导，而不是真正“一起表达”的样本主导。
这种情况下用皮尔逊相关去算基因间关系，结果会被大量0-0配对拖向一个非常虚高的正相关。也就是说，很多共表达关系其实是“都测不到”造成的假象，而不是生物学上的共调控。我见过不少初跑单细胞WGCNA的人，模块热图漂亮得不行，一查marker发现全落在“低表达-低表达”的假相关里，这种模块拿去做功能富集，富出来一堆泛组织表达的基础代谢通路，毫无辨识度。
1.3 伪bulk是妥协，不是解法
既然细胞级不可靠，最常见的替代方案就是伪bulk：把同一个细胞类型、同一个样本来源的细胞表达量全部加和取平均，得到一张样本×基因的矩阵，再用经典WGCNA流程分析。这样做的好处是计算稳定、统计合理，坏处是丢失了单细胞数据最大的价值——异质性。
同一个细胞类型内部往往藏着若干功能状态或分化中间态，伪bulk把这些状态揉成一个平均值，相当于把调控网络的“细颗粒度”全部磨平。不同状态共享的基因可能被平均成一个模块，状态特异的基因因为信号太弱直接被过滤掉。最后你得到的模块，更接近这个细胞类型的“平均个性”而不是“调控逻辑”。
2. hdWGCNA的破局思路：Meta-cell到底在聚合什么
2.1 在PCA空间里做KNN聚合
hdWGCNA的核心策略很简单：与其把几万个细胞直接塞进网络，不如先让“长得像”的细胞抱团成元细胞，再用元细胞作为分析的基本单位。具体做法是在某个细胞类型内部，先用PCA把高维表达谱压缩到低维空间，然后对每个细胞找它的k个最近邻，将这k个细胞的表达谱平均成一条元细胞表达谱。
拿k=25来举例：一个包含1万个细胞的类型，大约会产生400个左右的元细胞。如果你有八个细胞类型，全局元细胞数量大约三千个。这个量级对WGCNA算法来说非常舒服，基因相关矩阵的统计推断恢复了意义，同时计算量也能被普通工作站接受。
这个思路可以类比成民意调查。普查几万人的意见既费时又容易引入各种操作误差，更聪明的做法是按地理或人群特征先分层，再在每一层里抽有代表性的样本。元细胞就是这个采样单元，区别在于它是由数据自己通过邻近关系决定的，而不是由研究者拍脑袋划分的。
2.2 元细胞和伪bulk的根本区别
很多人第一次听到元细胞，会觉得这不就是缩小版的伪bulk吗？还真不是。伪bulk是对整个细胞类型做一次平均，得到的是该类型的一张“标准脸”；元细胞是在类型内部对局部邻域做平均，得到的是若干张“局部脸谱”。
k=25意味着每个元细胞只覆盖该细胞类型里极小的一个局部区域。如果这个类型内部存在从静息态到活化态的连续谱系，伪bulk会把中间过程平均掉，元细胞则能保留从静息端到活化端的一系列中间状态。共表达网络恰恰需要这些中间状态来“拉开”基因表达的变化幅度，否则相关性计算没有足够动态范围，模块结构就出不来。
下面这个对比表格应该能帮你更直观地理解两者的差异：
| 维度 | 伪bulk | hdWGCNA元细胞 |
|---|---|---|
| 聚合范围 | 整个细胞类型 | 类型内部的局部邻域 |
| 样本数量 | 类型数×样本数，通常只有几十 | 几百到几千，取决于k值 |
| 异质性保留 | 几乎不保留 | 保留局部状态和亚群梯度 |
| 对dropout的抑制 | 强（大量细胞平均） | 中等偏上（25个左右细胞平均） |
| 下游网络分析 | 可直接用经典WGCNA | 仍用WGCNA核心算法，但参数需调整 |
2.3 为什么聚合之后共表达信号反而更清晰
你可能觉得平均会让信号变模糊，实际上在单细胞场景里恰恰相反。dropout是随机的技术损失，一个基因在25个相似细胞中有15个检测到、10个丢失，平均之后表达值趋近于真实水平，技术噪声被压下去了。生物学共表达信号是系统性的，不会因为平均而消失；技术零值是随机性的，平均之后被大幅削弱，信噪比自然就上来了。
还有一层容易被忽略的好处：基因表达量经元细胞聚合后，动态范围扩大了，分布形态也更接近bulk RNA-seq。WGCNA默认假设表达数据近似正态或至少具有稳定的方差结构，原始单细胞数据的零膨胀分布几乎不符合这个前提，聚合后的元细胞矩阵让这个前提重新成立。这就是为什么hdWGCNA识别出的模块往往比直接拿伪bulk做出来的模块更有生物学解释力。
3. 实操流程：从Seurat对象到WGCNA输入矩阵
3.1 安装和前置条件
hdWGCNA目前没有登录CRAN和Bioconductor，需要从官方GitHub仓库安装开发版本，安装前先确保Seurat、WGCNA、tidyverse、matrixStats、pheatmap、igraph这些依赖库已经是可用状态。如果你是在RStudio里使用，也可以直接下载官方发布页的压缩包，通过“Install from: Package Archive File”手动安装，效果一样。
运行hdWGCNA的前提是你的Seurat对象至少满足两个条件：一是标准化已经完成，assay里存的是normalized data而不是raw counts；二是细胞类型注释列已经存在，后面meta-cell聚合全靠这个分组信息。我习惯在分析前把低质量细胞、双细胞和线粒体基因比例过高的细胞都清干净，这个步骤做不好，后面模块质量一定会受影响。
3.2 SetupForWGCNA：用表达比例先把基因筛一遍
第一行代码通常是SetupForWGCNA。这里最重要的参数是gene_select和fraction。绝大多数教程推荐的gene_select="fraction"意思是按“有多少比例的细胞检测到了这个基因”来过滤。比如fraction=0.05，就要求至少5%的元细胞里有该基因的检出，才能进入后续分析。
这个过滤逻辑是必要的。单细胞矩阵里有大量只在极少数细胞中“闪现”的基因，这些基因的表达跨度几乎全被零值主导，放进网络里只会产生大量虚假相关，还会拖慢计算速度。用fraction来筛比直接用平均表达量筛更合理，因为它保证留下来的基因在群体中具有足够的检测覆盖率，相关系数才有计算基础。
library(hdWGCNA) # 假设 pbmc 是已经做过标准化和注释的 Seurat 对象 pbmc <- SetupForWGCNA( pbmc, gene_select = "fraction", fraction = 0.05, wgcna_name = "pbmc_hdwgcna" )wgcna_name是一个标签字符串，本套分析的所有中间结果都会挂在pbmc@misc$pbmc_hdwgcna下面。这意味着你完全可以在同一个Seurat对象上跑好几套不同参数的hdWGCNA，互不覆盖，后续对比时再逐个取出来。
3.3 MetacellsByGroups：决定元细胞长什么样
真正决定分析成败的是MetacellsByGroups这一步。group.by指定的是Seurat对象里的细胞类型注释列，例如“cell_type”。hdWGCNA会在每个类型内部独立做KNN聚合，不同细胞类型不会混在一起产生跨群元细胞，这样既保持了细胞类型的隔离，又让类型内部的局部结构得到保留。
k值控制每个元细胞包含多少细胞，官方文档里默认是25，我自己的经验是15到40之间都算合理区间。k太小，比如5，dropout压制不足，元细胞和单细胞差别不大；k太大，比如100，元细胞数量骤降，局部异质性又开始被抹平。如果你不确定，可以跑两个k值,分别做一次完整分析，然后比较模块的稳定性和富集结果，选解释力更好的那套。
max_shared=10这个参数容易被忽略：它限制两个元细胞之间最多可以重叠10%的细胞，防止过度重叠导致元细胞之间不独立。你可以把max_shared理解成“相邻民意调查样本之间允许的重复受访者比例”，太高会让后续相关性的自由度假性膨胀。
pbmc <- MetacellsByGroups( pbmc, group.by = "cell_type", k = 25, max_shared = 10, ident.group = "cell_type" )跑完后可以检查一下生成的元细胞对象：到底产生了多少个元细胞、每个类型各占多少、表达矩阵稀疏度降到什么水平。这些检查看着琐碎，但能帮你提前发现注释列错误或者k值不合理的迹象。
3.4 SetDatExpr和TestSoftPowers：把矩阵喂给网络
SetDatExpr负责真正构建用于WGCNA的元细胞表达矩阵。这里最关键的是group_name参数，它决定这一轮分析聚焦哪一类细胞。hdWGCNA允许你只挑某个感兴趣的细胞类型单独建网络，这种做法通常比全员一起跑更合理，因为不同细胞类型的调控逻辑差异很大，强行放在同一张网络里只会得到一堆“平均调控模式”的模块。
pbmc <- SetDatExpr( pbmc, group_name = "CD4_T", group.by = "cell_type", assay = "RNA", slot = "data" ) pbmc <- TestSoftPowers(pbmc, networkType = "signed")TestSoftPowers会在不同的软阈值候选值下逐个测试，计算每个阈值对应的无标度拟合指数，为下一步的网络构建提供参考。它输出的结果里包含每个软阈值对应的拟合R²、平均连接度等指标，选阈值这件事我们下一章专门展开，这里先知道这个函数的目的就够了。
4. 软阈值、拓扑重叠与模块识别：解读中间结果
4.1 软阈值β到底在干什么
WGCNA不是直接拿基因表达相关系数当边的权重，而是先对相关系数做幂次变换：两个基因间的邻接度等于它们相关系数绝对值的β次方。这么做是想让网络逼近无标度网络结构——少数hub基因连接大量节点，多数基因只连接少数节点，这种结构在真实生物网络中非常普遍。
软阈值β越小，权重分配越均匀，网络越像密集的“毛线团”；β越大，只有相关性足够强的基因对才能保留有意义的边权重，网络就越稀疏、模块边界越清楚。单细胞和元细胞数据里的相关系数整体比bulk数据低，因为残留的技术噪声仍然会拉低相关性，所以通常需要比bulk更高的β才能把无标度特征拟合出来。bulk数据常见β在4到8之间，hdWGCNA分析里β到6到12都不算夸张。
4.2 怎么读无标度拟合曲线
GetPowerTable可以取出TestSoftPowers的完整结果。每一行对应一个候选软阈值，有个关键列是拟合无标度拓扑的R²，另外还有平均连接度。选择原则我总结成一句话：在R²达到0.8以上的候选β里，选最小那个。
这个原则的原因很实际：β越高，对相关性的压制越狠，许多真实的弱共表达关系会被直接砍掉，导致模块数量增多、每个模块都很小、解释起来非常零碎。R²刚过0.8时，网络既保留了足够的真实结构，又没有过度稀疏化。有些数据集整体噪声大，R²怎么都到不了0.8，这种情况我一般会接受0.7以上的值，与此同时回上游检查数据质量。
fit <- GetPowerTable(pbmc) head(fit)还有个小技巧：观察平均连接度那一列。好的软阈值应该让平均连接度保持在一个合理的中间水平，每个基因平均连接几十个节点左右。如果平均连接度已经不到10，网络基本上被切成了大量孤立的“牙齿”，那说明β选得太激进了。
4.3 模块分割和树状图的读法
选定软阈值只是开始，ConstructNetwork这一步会做三件事：计算拓扑重叠矩阵，基于TOM做层次聚类，再用动态树切割算法把聚类树切成模块。
这里有必要解释一下TOM。普通的相关系数只考虑两个基因直接相关有多强，TOM还额外考虑了它们“共同邻居”的数量：两个基因共享的邻居越多，即使直接相关性不算强，也倾向于被划进同一个模块。这个思路非常贴近真实调控网络，因为同一个通路里的基因往往不一定两两直接相关，但会共同受到上游调控因子的影响，共享大量间接关联。
ConstructNetwork跑完后，用PlotDendrogram看聚类树，树的分支颜色块就是初步的模块。我拿到树状图第一眼会看两件事：一是模块数量是否合理，通常几到十几个是正常的，如果切出三四十个模块，基本可以判定软阈值偏高或基因筛选参数偏严；二是模块大小是否悬殊，如果有个模块占了全部基因的六成以上，往往是网络过于松散的信号，考虑降低β重跑。
另外，WGCNA允许同一模块内部再检查模块特征基因之间的相关性，如果两个“模块特征基因”之间的相关性接近1，说明这两个模块本质上应该合并。ModuleEigengenes跑完之后可以顺手看一眼这个矩阵，别让高度冗余的模块白白分散你的解释精力。
5. 从模块到hub基因：特征基因与连通度的完整计算链
5.1 模块特征基因：让几百个基因变成一个综合指标
当你有了稳定的模块划分，下一步就是计算模块特征基因ME。每个模块的ME可以理解为该模块内所有基因表达的“主成分合成的代表向量”，它浓缩了这个模块在每一个元细胞里的整体活跃程度。有了ME，你就可以把“某模块活跃”变成一个可以跟任何表型、细胞状态、空间位置做关联的数值指标。
举个例子：如果你发现某个模块的ME在活化状态的细胞里显著升高，在静息细胞里降低，那么整个模块几十个基因就可以被整体解读为“活化相关表达程序”，而不是你要把几十个基因的表达模式挨个描述一遍。这种降维能力是共表达网络分析最实用的地方。
5.2 kME是筛选hub基因的正确姿势
模块划好之后，另一个绕不开的问题是：模块里哪些基因是核心调控者？这里不能用“表达量高不高”来判断，而应该用kME，也就是模块成员基于ME的连通度。通俗地说，kME衡量的是某个基因与它所在模块整体表达模式的相关程度，相关度越高，说明这个基因越能代表模块的核心行为。
常规筛选hub基因的阈值可以定在kME绝对值大于0.6到0.8之间，同时模块内排名靠前。我个人的习惯是先用0.6做一道粗筛，再取每个模块kME排名前10的基因看是否符合预期。粗筛之后要做的不是直接采信，而是人工检查hub基因在单细胞里的表达分布，确认它不是因为某个稀有亚群的高表达而虚高。
pbmc <- ModuleEigengenes(pbmc, group.by = "cell_type") pbmc <- ModuleConnectivity(pbmc) PlotKMEs(pbmc)PlotKMEs会展示所有基因在不同模块里的kME分布，如果你发现某个模块里几乎所有基因的kME都在0.3以下，这个模块很可能是个“垃圾模块”，由一些彼此关联极弱的基因被硬凑在一起，应当谨慎解释。
5.3 功能富集很容易被背景基因集带偏
模块拿到手，大多数人第一件事就是做GO或KEGG富集。这里有一个非常普遍的错误：使用全基因组做背景。单细胞hdWGCNA分析里，输入网络的基因本来就经过fraction过滤，只保留了表达覆盖率足够高的基因，如果用全基因组当背景，富集结果的p值会被系统性低估，很多假阳性的通路会被顶上来。
正确做法是把你实际纳入网络的基因列表作为背景集。很多富集工具允许自定义背景，花一分钟设置好，结果的可信度会高很多。另外，富集结果出现多个高度重叠的通路时，优先看那些有明确细胞类型或状态指向性的通路，而不是泛泛的“RNA代谢过程”。模块的生物学身份最终要靠marker基因和已知调控因子来锚定，富集只是辅助判断。
6. 空间转录组：meta-cell思路可以搬到空间邻域上
6.1 空间数据到底特殊在哪
空间转录组数据的本质和单细胞完全不同：每个空间点位捕获的往往是多个细胞的混合转录组，基因检出更加稀疏，而且相邻点位之间的表达天然高度相似。这种空间自相关会把“共表达”和“空间邻近”混在一起，如果直接把点位当独立样本跑WGCNA，模块很有可能反映的是组织切片上的解剖分区，而不是真正的基因调控程序。
所以处理空间数据时，meta-cell的构建逻辑要调整。你可以用两种思路：一种是把空间上的相邻点位聚合起来，得到“空间元细胞”，这种元细胞能有效提升信号稳定性，同时保留空间位置信息；另一种是先按无监督聚类把空间点位分成若干结构域，再以结构域为分组变量构建元细胞，这样模块识别出来的就是“结构域特异表达程序”而非随机空间噪声。
6.2 空间对象的hdWGCNA流程
整体流程和单细胞版本类似，区别主要在于输入对象和分组策略。SetupForWGCNA可以直接接受Seurat的Spatial对象，然后你需要在group.by里指定空间聚类结果或者解剖注释列。如果没有任何注释信息，先用Seurat默认的空间聚类流程跑一遍，把空间聚类结果作为分组变量，再做元细胞聚合。
空间数据里我一般会把k适当调大，比如30到50。原因是每个空间点位的检测深度比单细胞更低，单次平均覆盖的细胞数也不确定，稍微大一点的聚合窗口能更有效地压制技术噪声。但注意，k太大也会让模块丧失空间分辨能力，模块的区域特异性会下降。
6.3 模块的空间可视化与区室富集判断
模块算完之后，最有价值的一个操作是把每个模块的ME映射到组织切片上。如果某个模块的ME集中在特定解剖区域，例如某个类型的组织结构周围，那么这个模块极有可能对应一种区域性的细胞协作程序。你完全可以用Seurat里现成的空间可视化函数，把模块ME当成一个“自定义特征”展示在组织坐标上，切片上一眼就能看出模块的组织学意义。
除了主观观察，我建议顺手做一个简单的空间富集检验：把空间点位按结构域注释分组，比较不同组别之间的ME分布差异。这样既能验证模块的区域特异性，又能避免“看起来富集”实际上只是几个极端点贡献的视觉错觉。
7. 我在实战里踩过的坑和验证套路
7.1 别拿SCT后的数据直接喂进去
SCTransform是个好工具，它的目标是同时完成归一化和方差稳定化，还能消除测序深度的影响。但问题在于SCT回归掉的那部分技术变异里，有一部分其实携带共表达信息。我在自己的数据上对比过，用常规LogNormalize流程处理的矩阵跑hdWGCNA，模块的结构更稳定，特征基因也更符合已知的marker注释；换成SCT后的矩阵，模块明显更松散，不少已知调控关系被拆散了。
官方对这个问题的建议也比较保守，默认只接受常规标准化之后的表达矩阵。如果你确实需要用SCT数据，至少也要用SCT的“corrected counts”而不是“scale data”作为输入，并且对结果做更严格的稳定性验证。
7.2 模块稳定性检验不能省
共表达模块是数据驱动的结果，天然存在一定程度的随机性。单细胞数据采样噪声又大，很可能换个随机种子或者换个批次数据，模块边界就变了。我的习惯是做一次重采样验证：从数据里随机抽取八成细胞，再完整跑一遍元细胞聚合和网络构建，看同一个模块里的核心hub基因是否还能稳定聚在一起。
如果一轮重采样下来，模块核心成员保存率不足六成，那么这个模块就很可疑，需要谨慎对待。你可以把它标记为“候选低稳定模块”，只做描述不做过度的机制解读。与其被别人审稿人问住，不如自己先把筛选做严。
7.3 空间数据要警惕“组织结构假模块”
空间转录组还有个隐形坑：相邻点位高度相似，两个基因如果恰好富集在同一个解剖区室，即使它们根本没有调控关系，也会因为空间共定位而表现出很强的共表达。这种模块反映的是“组织结构”而非“调控网络”。
我的判断方法是做空间自相关检验。把模块特征基因在空间上的分布算一个空间自相关指标，如果某个模块的指标高到离谱，同时它在各个解剖区域都有分布但彼此隔离，那么它更可能是结构性模块。此时可以下修k值，缩小聚合窗口，再看模块是否被拆散；如果模块依然稳定存在，才有底气把它当成真实的调控程序来解读。
还有一点经验：越是复杂的组织样本，越不要一刀切地把所有点位丢进同一个网络。按组织区域或者明确的细胞组成先分层，再分别构建网络，得到的结果通常更有针对性，也更容易定位到具体的细胞互作信号。
我现在做hdWGCNA分析，最常做的一件事其实是回头检查那些“太漂亮的模块”。如果一个模块的marker基因特别齐整、富集通路特别显著、hub基因特别有道理，我都会忍不住怀疑是不是某个批次效应或者空间伪影在背后撑着。基因共表达网络的魅力在于它能帮你从数据里提炼调控假设，但它的可信度永远建立在数据质量和验证策略之上。把聚合参数调清楚、把模块稳定性查清楚、把背景基因集选清楚，这三件事做扎实了，hdWGCNA给出的模块才能真正成为你下一步实验设计的起点。
