单细胞代谢分析这两年从锦上添花变成了很多免疫治疗课题的必答题。原因很直接同样一群T细胞表面标记物看起来一模一样但代谢状态可能天差地别——有的靠糖酵解猛冲有的靠脂肪酸氧化打持久战而这两类细胞对免疫检查点抑制剂的响应完全不同。问题在于常规的单细胞转录组分析很难直接读出代谢通路的活性因为代谢酶的mRNA丰度和实际通量之间并不是线性关系。COMPASS这个工具就是冲着这个痛点来的它不满足于某个代谢基因表达高低这种粗糙判断而是通过约束优化模型从单细胞表达谱推断每个细胞在几十条代谢通路上的通量潜力。如果你正在做肿瘤免疫、自身免疫病或者感染免疫的单细胞课题手头有一批处理过和未处理组的样本想找出代谢层面的差异靶点那这套思路值得花时间跑一遍。下面我会把COMPASS的底层逻辑、数据准备、运行细节、差异分析技巧以及我自己踩过的坑完整地拆开讲。1. 为什么常规代谢基因打分在单细胞里经常失灵1.1 代谢酶mRNA与代谢通量之间的脱节现象很多人第一反应是代谢通路活性不就是把这条通路里所有酶的基因表达加起来取平均吗AUCell、ssGSEA、Seurat的AddModuleScore都是这么干的。这个做法在bulk RNA-seq里勉强能用因为大量细胞平均之后转录调控和代谢需求的耦合关系会被放大。但到了单细胞层面问题就暴露了。单个细胞的mRNA总量只有几千到几万条代谢酶基因往往表达量很低很多酶的转录本在单个细胞里根本检测不到dropout率极高。你算出来的通路评分很大程度上反映的是这条通路里的基因有没有被检测到而不是这个细胞的代谢通量高不高。更麻烦的是代谢通量受翻译后修饰和别构调节控制酶的量不变活性可以翻好几倍。所以直接拿mRNA打分去推断代谢状态在单细胞分辨率下噪声极大。1.2 COMPASS的核心思路从表达量转向通量潜力COMPASSConstraint-based Modeling of Single-cell metabolic Activity的做法完全不同。它把每个细胞看作一个独立的代谢网络用基因组尺度代谢模型GEM通常是Recon系列模型作为骨架把单细胞表达数据映射到模型中的反应上然后通过约束优化求解每条代谢通路的最大可行通量。打个比方常规打分像是看一个工厂里有多少台机器酶的表达而COMPASS是看这些机器连成的生产线在原料供应和产能约束下最多能产出多少产品代谢通量。即使某台机器的数量检测不到只要上下游通路是通的模型仍然可能推断出这条通路有通量潜力。这就绕开了单细胞dropout的部分干扰。具体来说COMPASS对每个细胞求解一个类似FBA通量平衡分析的优化问题以代谢通路为单位最大化该通路的通量同时满足化学计量约束、热力学约束不可逆反应方向和表达约束酶表达低的反应通量上限被压低。最终每个细胞得到的是一个向量每个维度对应一条代谢通路的活性评分。1.3 COMPASS输出的是什么和普通打分有什么本质区别COMPASS的输出是一个细胞×代谢通路的矩阵值域通常在0到1之间经过归一化后表示该通路在该细胞中的相对通量潜力。和普通打分的本质区别有三点第一通路定义基于代谢模型而非基因集。COMPASS用的通路是代谢网络中的功能模块有明确的化学计量关系不是简单的基因列表。这意味着它考虑了代谢物在通路之间的流转而不是孤立地看基因表达。第二约束优化带来了非线性信息。一个反应的表达低不一定导致通路通量低因为可能有限速步骤之外的补偿。COMPASS通过优化框架捕捉了这种非线性关系。第三输出具有跨细胞可比性。因为每个细胞都用同一个代谢模型求解得到的通路评分在细胞之间是可比的这为后续的差异分析打下了基础。注意COMPASS不是万能的。它依赖代谢模型的质量而Recon等模型对免疫细胞的特异性代谢途径覆盖并不完美。如果你的课题聚焦某种特殊代谢比如色氨酸代谢在肿瘤微环境中的作用需要检查模型里相关反应是否完整。2. 跑COMPASS之前必须搞定的数据准备2.1 输入数据的格式要求与常见踩坑COMPASS的R包compassR对输入数据有明确要求。核心输入是一个基因表达矩阵行是基因Entrez ID或Symbol列是细胞值是标准化后的表达量。这里有几个坑我必须提前说坑一基因ID不匹配。COMPASS内部用的代谢模型是基于Entrez ID的如果你的表达矩阵用的是Symbol需要先转换。转换过程中会有大量基因丢失尤其是那些Symbol有多个Entrez ID对应的基因。我的建议是直接用Entrez ID从头跑上游定量或者在转换后用AnnotationDbi检查丢失比例如果超过15%要考虑换注释版本。坑二标准化方式选错。COMPASS对输入的标准化方式比较敏感。不建议用TPM或FPKM因为这两种标准化会引入基因长度偏差而代谢模型里的反应约束是基于相对表达量的。推荐用CPMcounts per million或者SCTransform后的校正counts。我自己测试下来CPM在COMPASS里的表现最稳定。坑三细胞数量太多导致计算爆炸。COMPASS对每个细胞求解一个优化问题细胞数超过5000之后计算时间会急剧增加。如果你有上万细胞建议先做亚群聚类每个亚群抽样500-1000个细胞跑COMPASS或者只对你关心的细胞类型比如CD8 T细胞跑。2.2 代谢模型的选择与反应过滤COMPASS默认使用Recon2模型但这个模型包含7000多个反应很多反应在你的细胞类型里根本不表达。直接跑全模型不仅慢还会引入噪声。我的做法是先根据表达矩阵筛选出在至少5%细胞中检测到的代谢基因对应的反应。用compassR提供的filterModel函数把不表达的反应从模型中剔除。保留核心通路糖酵解、TCA循环、氧化磷酸化、脂肪酸氧化、氨基酸代谢等这些是免疫细胞代谢分析的重点。过滤后的模型通常只剩1500-2500个反应计算速度能提升3-5倍而且结果更干净。2.3 批次效应的预处理策略单细胞数据几乎不可能没有批次效应而COMPASS对批次效应非常敏感。因为不同批次之间即使是同一群细胞检测到的代谢基因集合可能不同导致通量推断出现系统性偏差。我的处理流程是先用Harmony或Seurat的IntegrateLayers做批次校正拿到校正后的表达矩阵再输入COMPASS。注意不要用校正后的降维结果COMPASS需要的是表达矩阵。Harmony校正后的矩阵可以直接用但要注意Harmony可能会过度校正把真实的生物学差异也抹掉。建议在校正前后各跑一次COMPASS对比结果如果差异很大说明批次效应严重需要调整校正参数。3. COMPASS运行中的参数调优与计算加速3.1 核心参数解读lambda、penalty与迭代次数COMPASS有几个关键参数直接影响结果参数含义推荐值调整逻辑lambda正则化强度0.5-1.0值越大通量推断越保守适合噪声大的数据penalty对低表达反应的惩罚0.1-0.3值越大越倾向于忽略低表达反应n_iter优化迭代次数1000-2000太少不收敛太多浪费时间n_cores并行核心数根据机器建议至少8核lambda的调整需要根据数据质量来。如果你的数据是10x Genomics的常规深度每个细胞1000-3000个基因lambda设0.5-0.7比较合适。如果是Smart-seq2的高深度数据可以降到0.3-0.5让模型更自由地推断通量。3.2 并行计算与内存管理COMPASS的计算瓶颈在于对每个细胞求解优化问题。在8核机器上1000个细胞大约需要30-60分钟。如果你有5000个细胞单机跑可能需要4-6小时。几个加速技巧分块运行把细胞分成每块500个分别跑COMPASS最后合并结果。这样即使某一块出错也不用从头再来。内存预分配COMPASS的输出矩阵是细胞×通路提前用matrix(NA, nrown_cells, ncoln_pathways)分配好内存避免运行中反复扩容。使用稀疏矩阵如果表达矩阵很稀疏单细胞数据通常如此用Matrix包的稀疏矩阵格式输入能节省大量内存。3.3 结果的后处理与归一化COMPASS跑完之后每个细胞会得到一个通路活性向量。但这个向量的绝对值没有意义需要做归一化。我通常用两种方式第一种是分位数归一化把每个通路的活性值映射到0-1之间用该通路在所有细胞中的分位数。这样不同通路之间可以比较。第二种是参考细胞归一化选一群代谢状态稳定的细胞比如初始T细胞作为参考把其他细胞的通路活性除以参考细胞的均值。这样得到的是相对活性更适合做组间比较。提示归一化方式的选择会影响后续差异分析的结果。如果你的课题是找处理组 vs 对照组的差异通路建议用参考细胞归一化因为分位数归一化会抹掉组间的整体差异。4. 差异分析从COMPASS输出中挖出代谢靶点4.1 组间差异通路的统计检验方法拿到COMPASS矩阵后差异分析和常规单细胞分析类似但有几个特殊之处。首先COMPASS的通路活性值不是正态分布很多通路的值集中在0附近因为大部分细胞在某条通路上没有通量。所以不要直接用t检验建议用非参数检验Wilcoxon秩和检验或者专门为稀疏数据设计的检验方法。我的流程是对每个通路计算处理组和对照组的活性均值差异。用Wilcoxon检验计算p值。用Benjamini-Hochberg方法做多重检验校正得到FDR。筛选FDR 0.05且活性差异倍数 1.5的通路作为候选。这里要注意COMPASS的通路活性差异倍数不能直接算因为值域是0-1很多值接近0。我通常先做log2转换加一个小的伪计数再算差异倍数。4.2 单细胞水平的代谢异质性分析组间差异只是第一步。更有价值的是看同一组内不同细胞亚群之间的代谢异质性。比如在肿瘤浸润T细胞中你可能发现耗竭T细胞和效应T细胞在脂肪酸氧化通路上有显著差异而这个差异在组间比较中被平均掉了。具体做法是先对细胞做聚类用Seurat或Leiden然后对每个亚群计算COMPASS通路的平均活性做亚群间的差异分析。我经常用热图展示亚群×通路的活性矩阵一眼就能看出哪些亚群在哪些通路上有特征性活性。4.3 代谢靶点的优先级排序与验证思路从COMPASS输出到真正的代谢靶点中间还有一段路。我的优先级排序逻辑是第一优先级组间差异显著 在关键细胞亚群中特异性高 通路有已知的小分子抑制剂或激动剂。第二优先级组间差异显著 通路与免疫治疗响应相关需要查文献。第三优先级亚群间差异显著但组间不显著可能是细胞状态差异而非处理效应。验证思路方面如果条件允许可以用Seahorse代谢分析仪做功能验证或者用代谢物检测比如乳酸、ATP做生化验证。如果只有转录组数据可以用代谢通路的限速酶表达作为间接验证。5. 实操中容易翻车的几个细节5.1 细胞过滤阈值对结果的影响COMPASS对输入细胞的质量很敏感。如果你把低质量细胞检测基因数500线粒体比例20%也放进去这些细胞的代谢基因检测率极低COMPASS会推断出大量零通量通路拉低整体信号。我的建议是在跑COMPASS之前严格过滤细胞。检测基因数至少1000线粒体比例低于10%如果是组织样本可以放宽到15%。过滤后如果细胞数不够宁可多跑几个样本也不要降低阈值。5.2 通路冗余与结果解读的陷阱COMPASS输出的通路之间有大量冗余。比如糖酵解和葡萄糖转运是两条通路但活性高度相关。如果你不做处理差异分析结果里会出现一堆高度相关的通路看起来有很多靶点实际上是一个生物学过程。处理方法是先计算通路之间的相关性矩阵把相关系数0.8的通路聚成一类每类只保留一个代表性通路通常是文献里最常报道的那个。这样能把候选通路从几十个压缩到十几个更聚焦。5.3 与常规代谢分析结果的交叉验证COMPASS的结果不要孤立地看。我通常会把它和常规的代谢基因打分比如AUCell做交叉验证。如果COMPASS推断某条通路活性高而该通路的限速酶表达也高那这个结果就比较可靠。如果两者矛盾需要仔细检查是COMPASS的模型约束有问题还是限速酶的表达被dropout影响了。另外也可以和代谢物数据如果有的话做交叉验证。比如COMPASS推断氧化磷酸化活性高而代谢物检测显示ATP水平高那就互相印证了。6. 从COMPASS结果到可验证的代谢靶点一个完整的分析链条6.1 候选通路的筛选与可视化假设你跑完COMPASS得到了处理组和对照组的通路活性矩阵。下一步是筛选候选通路。我通常分三步走第一步做组间差异分析得到差异通路列表。第二步对差异通路做聚类把功能相关的通路归为一类。第三步对每一类选一个代表通路画小提琴图或箱线图展示组间差异。可视化方面我推荐用ggplot2的geom_violin加geom_boxplot的组合既能看分布又能看中位数差异。如果通路很多可以用热图展示所有通路的组间差异倍数用颜色深浅表示差异大小。6.2 代谢靶点的功能注释与通路富集筛选出候选通路后需要做功能注释。COMPASS的通路名称通常比较专业比如脂肪酸beta氧化需要映射到KEGG或Reactome通路才能做富集分析。我的做法是把候选通路的名称手动映射到KEGG通路ID然后用clusterProfiler做富集分析看这些通路集中在哪些代谢过程中。如果富集到氧化磷酸化或糖酵解这种大通路说明你的候选靶点集中在能量代谢上可以进一步聚焦。6.3 从代谢通路到可成药靶点的转化逻辑最后一步也是最关键的一步从代谢通路找到可成药的靶点。不是所有代谢通路都有小分子抑制剂所以需要查数据库比如DrugBank、TTD看哪些通路里的酶有已知的药物。我的经验是优先关注这几类靶点糖酵解通路里的HK2、PKM2脂肪酸氧化通路里的CPT1A氨基酸代谢通路里的IDO1、GLS。这些靶点都有临床或临床前阶段的抑制剂转化潜力大。如果COMPASS推断某条通路活性高但该通路没有已知药物可以考虑两个方向一是看该通路的上下游是否有可成药靶点二是看该通路是否可以作为生物标志物而非治疗靶点。7. 一些个人经验与后续扩展方向7.1 计算资源有限时的替代方案如果你没有高性能服务器跑COMPASS确实吃力。我的替代方案是先用常规代谢基因打分做初筛找出差异最大的几条通路然后只对这些通路跑COMPASS。compassR支持指定通路子集这样计算量能减少80%以上。另一个方案是用scFEA或Compass的Python版本后者在某些场景下速度更快。但要注意不同工具的代谢模型和优化算法有差异结果不能直接混用。7.2 COMPASS与其它单细胞代谢工具的搭配使用COMPASS不是唯一的单细胞代谢分析工具。scFEA、Metabolite、scMetabolism各有侧重。我的搭配策略是用scMetabolism做快速初筛它基于KEGG通路速度快用COMPASS做精细推断它基于代谢模型更准确两者结果交叉验证。如果两者都指向同一条通路那这条通路的可信度就很高。如果只有COMPASS推断出来需要额外验证因为可能是模型约束带来的假阳性。7.3 后续可以扩展的分析维度跑完COMPASS只是起点。后续可以扩展的方向很多一是结合TCR序列数据看特定克隆型的T细胞是否有特征性代谢状态二是结合转录因子分析比如SCENIC看哪些转录因子驱动了代谢异质性三是做代谢通路的细胞间通讯分析看肿瘤细胞和免疫细胞之间是否存在代谢物交换。我个人最看好的是第一个方向。因为TCR克隆型和代谢状态的关联能直接回答为什么某些克隆型的T细胞在免疫治疗后扩增得更好这个问题临床意义很大。最后分享一个小技巧COMPASS跑完之后不要急着做差异分析。先把所有细胞的通路活性矩阵做一次PCA或UMAP看看细胞在代谢空间里是怎么分布的。很多时候代谢空间的聚类和转录组空间的聚类不一致这种不一致本身就是有意思的生物学信号。我就在一次分析中发现转录组上归为同一群的CD8 T细胞在代谢空间里分成了两群后来验证发现一群是真正具有杀伤功能的效应细胞另一群是准备进入耗竭状态的过渡态细胞。这个发现直接改变了后续实验的设计方向。