1. 为什么单细胞代谢分析值得你花时间折腾单细胞测序技术这几年从科研小众变成了免疫学、肿瘤学、发育生物学领域的标配工具。大家最开始关注的都是细胞分群、marker基因、拟时序这些常规操作但做着做着就会发现一个尴尬的问题你分出来的亚群除了表面标记物不同它们在功能层面到底差在哪里尤其是免疫治疗背景下为什么有的T细胞亚群对PD-1阻断响应好有的却完全没反应这个问题的答案很多时候藏在代谢层面。传统的单细胞转录组分析能告诉你某个亚群高表达什么基因但基因表达和代谢通量之间并不是一一对应的关系。一个代谢通路里的酶基因表达上调了不代表这条通路真的在高速运转因为代谢反应还受到底物浓度、产物反馈抑制、酶活性的翻译后调控等多重因素影响。这就是为什么我们需要COMPASS这类专门做代谢通量估计的算法。COMPASS的全称是ComputationalOptimization ofMetabolicActivity fromSingle-cellSequencing翻译过来就是“从单细胞测序数据计算优化代谢活性”。它的核心思路是把单细胞转录组数据映射到一个基因组尺度的代谢网络模型上然后通过线性规划求解每个细胞在给定营养条件下各条代谢反应的通量分布。说白了它不是在看你表达了多少代谢酶而是在算你这些酶组合起来能跑出多大的代谢流量。这个思路对于免疫治疗研究特别有价值。比如你想知道肿瘤浸润T细胞的代谢状态是不是导致免疫治疗耐药的原因用COMPASS跑一遍就能看到不同T细胞亚群在糖酵解、氧化磷酸化、脂肪酸氧化、氨基酸代谢等通路上的通量差异。这些差异往往比单纯的基因表达差异更能解释功能表型。适合读这篇内容的人已经做过基础的单细胞分析、想往代谢方向拓展的科研人员正在做免疫治疗相关课题、需要找新靶点的研究者以及任何对单细胞数据挖掘有需求、想了解COMPASS这个工具的人。我会从算法原理讲到实操步骤再到差异分析和靶点筛选尽量把踩过的坑都标出来。2. COMPASS算法到底在算什么2.1 从基因表达到代谢通量的映射逻辑COMPASS的核心是一个基于约束的代谢模型。你可以把它想象成一个城市的交通网络道路就是代谢反应路口就是代谢物车流量就是通量。每条道路的通行能力是有限的对应酶的最大催化速率而且进入路口的车必须等于离开路口的车对应代谢物的质量平衡。具体来说COMPASS使用的是一个经过整理的基因组尺度代谢模型通常是人类代谢网络Recon系列模型。这个模型包含了数千个代谢反应和代谢物覆盖了中心碳代谢、氨基酸代谢、脂质代谢、核苷酸代谢等主要通路。算法的输入是单细胞的基因表达矩阵。它首先根据每个基因的表达水平来设定对应代谢反应的容量约束。这里有一个关键步骤COMPASS不是简单地用表达量乘以一个系数而是通过一个基因-反应映射规则GPRGene-Protein-Reaction association来确定。比如一个反应需要酶A和酶B同时存在才能催化那么它的容量就取决于A和B表达量的较小值如果酶A或酶B任意一个存在就能催化那就取较大值。这个映射过程决定了COMPASS对输入数据质量比较敏感。如果你的单细胞数据dropout率很高很多代谢酶基因检测不到那COMPASS估计出来的通量就会偏低甚至为零。所以实际操作中数据预处理和基因覆盖度检查是第一步要做的。2.2 线性规划求解与通量分布设定好容量约束之后COMPASS会为每个细胞求解一个线性规划问题。目标函数通常是最大化某个生物量反应的通量或者最大化整体代谢网络的活跃程度。约束条件包括每个代谢物的净通量为零稳态假设每条反应的通量不超过其容量上限热力学不可逆反应的通量方向固定求解这个线性规划问题就能得到每个细胞在代谢网络上的一个通量分布向量。这个向量的维度等于代谢反应的个数通常有几千维。直接拿这个高维向量做下游分析不太现实所以COMPASS还会计算一个代谢通路活性评分把通量按通路聚合得到每个细胞在各条代谢通路上的活性值。这里有一个实操中容易忽略的点COMPASS的求解是随机规划的也就是说同一个细胞跑两次可能得到略微不同的通量分布。这是因为线性规划的最优解可能不唯一COMPASS会在最优解空间里随机采样。所以如果你要做差异分析建议设置足够的采样次数通常1000次以上然后取平均值或者中位数作为该细胞的代谢活性估计。2.3 COMPASS与其他代谢分析工具的区别市面上做单细胞代谢分析的工具不止COMPASS一个还有scFEA、Compass注意大小写不同是另一个工具、MetabolicActivity等。它们之间的核心区别在于工具核心方法输出适用场景COMPASS约束优化线性规划反应水平通量分布需要精确通量估计、通路间交互分析scFEA神经网络代谢网络通路活性评分大规模数据、快速筛选MetabolicActivity基因集评分通路评分快速探索、初步筛选Compass通量平衡分析通量值微生物代谢、特定模型COMPASS的优势在于它的机制可解释性。因为它是基于真实的代谢网络约束来求解的所以得到的通量值有明确的生化意义。你可以追踪某个代谢物是从哪条通路来的、到哪条通路去的。这对于找靶点特别重要因为你要找的不是一个统计上显著的差异而是一个有生物学机制支撑的干预点。但COMPASS的劣势也很明显计算量大。一个包含几千个细胞的数据集跑COMPASS可能需要几个小时甚至更久取决于你设置的采样次数和并行核数。所以如果你的数据量很大建议先做细胞亚群抽样每个亚群取几百个细胞代表跑完COMPASS再映射回全体。3. 从原始数据到COMPASS输入完整实操流程3.1 数据准备与质控要点假设你已经有了一个经过基础处理的单细胞Seurat对象或者Scanpy的AnnData对象。在跑COMPASS之前有几件事必须确认第一基因注释版本要统一。COMPASS使用的代谢模型是基于特定版本的基因注释通常是Entrez ID或Ensembl ID。如果你的单细胞数据用的是基因symbol需要先做ID转换。我建议统一转成Entrez ID因为Recon模型就是用Entrez ID索引的。转换工具可以用org.Hs.eg.db人类或org.Mm.eg.db小鼠注意处理多对一和一对多的映射关系。第二检查代谢基因的覆盖度。人类代谢模型涉及大约2000-3000个代谢基因。你可以先算一下你的数据里检测到了多少个。如果覆盖率低于60%COMPASS的结果可靠性会打折扣。提高覆盖度的方法包括使用更宽松的质控阈值比如min.cells3而不是10、考虑使用imputation方法补dropout但要谨慎imputation可能引入假信号。第三营养条件设定。COMPASS需要你指定培养基的营养成分包括葡萄糖、谷氨酰胺、氧气等的浓度。这个参数直接影响通量分布。如果你做的是肿瘤微环境研究可以设置低葡萄糖、低氧的条件来模拟肿瘤内环境如果是体外培养实验就用标准培养基的浓度。这个参数没有绝对正确的值但一定要在文章方法里写清楚。# 检查代谢基因覆盖度示例R library(org.Hs.eg.db) metabolic_genes - read.csv(recon_metabolic_genes.csv) # 从模型文件提取 entrez_ids - mapIds(org.Hs.eg.db, keys rownames(seurat_obj), column ENTREZID, keytype SYMBOL) detected - sum(entrez_ids %in% metabolic_genes$EntrezID) coverage - detected / length(metabolic_genes$EntrezID) cat(代谢基因覆盖率:, round(coverage * 100, 1), %\n)3.2 构建COMPASS可识别的输入对象COMPASS的R包COMPASS接受的是一个特定格式的列表对象包含表达矩阵、基因ID映射、营养条件等。构建这个对象的过程有点繁琐但每一步都有明确的逻辑。首先你需要从单细胞对象中提取原始计数矩阵不是归一化后的数据。COMPASS内部会自己做归一化如果你提前归一化了反而会干扰它的容量估计。提取的时候注意保持基因和细胞的对应关系。然后你需要准备一个基因-反应映射表。这个表可以从Recon模型文件里提取包含每个反应对应的基因ID和GPR规则。COMPASS包通常会自带这个映射表但如果你用的是自定义模型就需要自己构建。# 构建COMPASS输入对象R伪代码 library(COMPASS) # 提取原始计数 counts_matrix - GetAssayData(seurat_obj, slot counts) # 准备反应-基因映射 reaction_gene_map - prepare_reaction_gene_map(model recon2) # 设定营养条件 nutrient_conditions - list( glucose 5, # mM glutamine 2, # mM oxygen 0.1, # 相对浓度低氧模拟肿瘤 serine 0.5, glycine 0.5 ) # 构建COMPASS对象 compass_input - COMPASSContainer( data counts_matrix, metadata seurat_objmeta.data, reaction_gene_map reaction_gene_map, nutrient_conditions nutrient_conditions )这里有一个容易踩的坑COMPASS对输入矩阵的稀疏性很敏感。如果你的数据里有很多零值COMPASS可能会把某些反应的容量设为零导致整条通路不通。解决办法是在构建对象之前先对表达矩阵做一次轻度的平滑处理比如用Seurat::SmoothData()或者简单的kNN平滑。但注意不要过度平滑否则会抹掉真实的生物学差异。3.3 运行COMPASS与参数调优COMPASS的核心函数是COMPASS()主要参数包括n_samples每个细胞的采样次数默认1000。建议至少1000如果计算资源允许可以设到5000。n_cores并行核数。COMPASS支持多核并行建议设为可用核数的80%。timeout单个细胞求解的超时时间。如果某个细胞特别难求解可以设一个上限避免卡死。verbose是否输出进度信息。第一次跑建议设为TRUE方便监控。# 运行COMPASS compass_result - COMPASS( data compass_input, n_samples 2000, n_cores 8, timeout 300, verbose TRUE )跑完之后你会得到一个包含每个细胞代谢通量分布的列表。这个结果对象比较大建议及时保存为RDS文件避免重复计算。saveRDS(compass_result, file compass_result_full.rds)参数调优的经验如果你发现很多细胞的求解状态是infeasible不可行说明你的营养条件设得太苛刻了或者基因覆盖度太低。可以尝试放宽营养条件比如提高葡萄糖浓度或者检查一下是不是有大量代谢基因没检测到。另一个常见问题是求解时间过长这时候可以适当降低n_samples或者先用一个子集测试参数确认没问题再跑全量。4. 差异分析与代谢靶点筛选实战4.1 从通量分布到通路活性评分COMPASS跑完之后你拿到的是每个细胞在每个代谢反应上的通量值。这个数据维度太高直接做差异分析容易假阳性。所以第一步是降维把反应水平的通量按代谢通路聚合得到通路水平的活性评分。聚合的方法有几种可以简单求和也可以取平均还可以用通量变异系数来加权。我通常用的是通路内反应通量的中位数因为中位数对异常值更稳健。聚合的映射关系可以从代谢模型文件里提取每个反应属于哪个子系统subsystem都有标注。# 通路活性聚合R伪代码 pathway_activity - aggregate_pathway_activity( compass_result compass_result, model recon2, method median ) # 结果是一个细胞×通路的矩阵 dim(pathway_activity) # [1] 5000 80 假设5000个细胞80条通路得到通路活性矩阵之后就可以做常规的差异分析了。如果你有分组信息比如响应组vs非响应组可以用Wilcoxon秩和检验或者limma来做差异通路分析。注意做多重检验校正FDR 0.05作为阈值。4.2 差异代谢通路的可视化与解读差异分析做完之后怎么展示结果也很关键。我常用的几种可视化方式第一种通路活性热图。把差异显著的代谢通路挑出来画一个细胞×通路的热图按分组排序。这样能直观看到哪些通路在哪个组里更活跃。第二种火山图。横轴是通路活性的差异倍数log2FC纵轴是显著性-log10 FDR。可以标注出最显著的那几个通路。第三种代谢网络图。把差异通路映射到代谢网络图上用颜色表示上调或下调。这种图能帮你看到通路之间的关联比如糖酵解上调的同时乳酸分泌通路是不是也上调了。第四种单细胞水平的通量分布图。比如你想看糖酵解通路在T细胞亚群间的差异可以画一个小提琴图或者箱线图每个点代表一个细胞的糖酵解活性。这里有一个解读上的注意事项COMPASS估计的是相对通量不是绝对通量。也就是说你只能比较不同细胞或不同组之间的相对高低不能说某个细胞的糖酵解通量是多少毫摩尔每秒。所以在写文章的时候措辞要准确用higher glycolytic flux而不是increased glycolytic rate by X-fold。4.3 从差异通路到可干预靶点的筛选逻辑找到差异代谢通路只是第一步真正的目标是找到可干预的靶点。这里的逻辑链条是差异通路 → 关键酶 → 可成药性 → 实验验证。具体来说你可以按以下步骤筛选第一步锁定差异通路中的关键酶。比如你发现氧化磷酸化通路在耐药组中显著下调那就去看这条通路里哪些酶的基因表达也同步下调。COMPASS的结果里其实包含了反应水平的通量你可以直接提取氧化磷酸化相关反应的通量值看哪些反应差异最大。第二步评估可成药性。不是所有代谢酶都有现成的小分子抑制剂。你可以查一下DrugBank、TTD等数据库看目标酶有没有已知的抑制剂或者正在临床试验的药物。优先选择那些已有工具化合物的靶点这样后续验证会快很多。第三步检查靶点的表达特异性。理想情况下你希望靶点在肿瘤细胞或特定免疫细胞亚群里高表达而在正常组织里低表达。可以用GTEx或者Human Protein Atlas的数据来交叉验证。第四步设计体外验证实验。常见的验证思路包括用siRNA或CRISPR敲低目标酶看是否逆转代谢表型用抑制剂处理细胞看是否影响免疫功能比如T细胞的杀伤能力、细胞因子分泌等做代谢流实验如Seahorse确认通量变化。我自己的经验是不要一上来就盯着最显著的那条通路。有时候最显著的通路是下游效应不是驱动因素。反而是一些中等显著但处于网络枢纽位置的通路更值得关注。你可以用通路间的相关性分析来识别枢纽通路如果一条通路和很多其他差异通路都高度相关那它可能是上游调控点。5. 常见问题与排查技巧实录5.1 COMPASS运行报错与解决方案问题一求解器报错infeasible problem。这是最常见的问题通常是因为营养条件设得太严格或者某些必需反应的基因表达为零。解决办法先检查营养条件把葡萄糖和谷氨酰胺浓度适当调高然后检查基因覆盖度如果某些关键代谢基因缺失可以考虑用同源基因替代或者放宽GPR规则。问题二运行时间过长。COMPASS的计算复杂度是O(n_cells × n_reactions × n_samples)。如果你的细胞数超过1万建议先做细胞抽样每个亚群取200-500个代表细胞。或者降低n_samples到500虽然精度会下降但对于初步探索够用了。问题三结果中大量细胞的通量为零。这通常是因为输入矩阵太稀疏。可以尝试用ALRA或SAVER做imputation但要注意imputation可能引入假阳性。另一个办法是聚合细胞把同一亚群的细胞表达谱平均一下当成一个元细胞来跑COMPASS这样能提高信噪比。问题四不同批次的数据结果不可比。COMPASS的结果受批次效应影响很大因为不同批次的测序深度和基因检测率不同。解决办法在跑COMPASS之前先做批次校正如Harmony、Seurat Integration或者把批次作为协变量纳入下游差异分析模型。5.2 差异分析中的统计陷阱陷阱一伪重复问题。如果你的每个处理组只有一个样本但有很多细胞直接把细胞当独立样本做检验会严重高估显著性。正确的做法是用混合效应模型把样本作为随机效应或者用伪批量pseudobulk方法先把细胞聚合成样本水平的均值再做检验。陷阱二多重检验校正方法选择。COMPASS输出的通路有几十到上百条做多重检验校正是必须的。但不同校正方法Bonferroni、BH、Storey的结果可能差异很大。我通常用BH方法因为它控制FDR的同时保留了较好的检验效能。如果通路数很少20可以用Bonferroni更严格。陷阱三忽略代谢通路的层级结构。代谢通路不是独立的糖酵解和氧化磷酸化通过丙酮酸脱氢酶连接氨基酸代谢和核苷酸代谢通过一碳单位连接。做差异分析时如果只盯着单条通路可能会错过通路间的协同变化。建议做一次通路共表达网络分析看看哪些通路模块在组间差异最大。5.3 从COMPASS结果到实验验证的衔接COMPASS给的是计算预测最终还是要实验验证。这里分享几个衔接上的经验经验一优先验证通量差异最大的反应而不是通路。通路水平的聚合会平滑掉一些细节有时候通路整体不显著但通路里的某条关键反应显著。你可以直接从COMPASS结果里提取反应水平的通量做差异分析往往能发现更有意思的靶点。经验二用代谢物检测来交叉验证。如果COMPASS预测糖酵解通量上调你可以用质谱检测乳酸分泌量来验证如果预测脂肪酸氧化下调可以检测酰基肉碱水平。代谢物数据比转录组更接近表型验证价值更高。经验三注意COMPASS的预测方向可能和基因表达相反。这听起来很奇怪但确实会发生。比如某个代谢酶的基因表达上调了但COMPASS预测的通量反而下降可能是因为底物供应不足或者产物反馈抑制。这种情况反而更有意思说明存在翻译后调控或者代谢物介导的调控值得深入挖掘。经验四不要忽视低通量但高特异性的通路。有些代谢通路比如色氨酸-犬尿氨酸通路在整体通量中占比很小但在免疫调节中作用关键。COMPASS的通量估计对这些低丰度通路可能不够敏感你可以单独提取相关反应的通量值做更细致的分析。6. 一些实操中的个人体会COMPASS这个工具我从两年前开始用踩过的坑不算少。最大的一个体会是不要把它当成一个黑箱。很多人跑完COMPASS拿到通路活性矩阵就直接扔进差异分析流程然后对着结果硬解释。这样很容易得出似是而非的结论。我的做法是每次跑完COMPASS都会先做一轮质量检查看看通量分布是否合理比如糖酵解通量在活跃增殖的细胞里应该偏高、看看有没有大量零值、看看不同批次的分布是否一致。这些检查花不了多少时间但能避免很多后续的麻烦。另一个体会是COMPASS的结果要和基因表达数据结合看。如果COMPASS预测某条通路通量上调但通路里的酶基因表达没变化那可能是酶活性调控如果酶基因表达也上调那就是转录调控。这两种情况的干预策略完全不同前者可能需要找激酶或者变构调节剂后者可以直接用转录抑制剂或者RNAi。还有一点营养条件的设定对结果影响很大。我试过同一批数据用高糖和低糖条件跑COMPASS得到的差异通路列表重叠度只有60%左右。所以如果你的研究涉及肿瘤微环境或者免疫代谢一定要根据实际生理条件来设定营养参数不要用默认值。最后关于计算资源如果你的数据量很大建议用服务器跑至少32GB内存起步。COMPASS的R包对内存管理不是特别优化细胞数超过2万的时候可能会爆内存。一个变通的办法是分批次跑每次跑5000个细胞然后把结果合并。这个方向后续还可以往空间代谢组结合的方向扩展。现在有空间转录组数据了如果能结合COMPASS做空间代谢通量估计就能看到代谢异质性在组织空间上的分布这对于理解肿瘤微环境里的免疫代谢互作会很有帮助。不过这是另一个大话题了有机会再展开聊。
单细胞代谢分析实战：COMPASS算法原理、实操与免疫治疗靶点筛选
1. 为什么单细胞代谢分析值得你花时间折腾单细胞测序技术这几年从科研小众变成了免疫学、肿瘤学、发育生物学领域的标配工具。大家最开始关注的都是细胞分群、marker基因、拟时序这些常规操作但做着做着就会发现一个尴尬的问题你分出来的亚群除了表面标记物不同它们在功能层面到底差在哪里尤其是免疫治疗背景下为什么有的T细胞亚群对PD-1阻断响应好有的却完全没反应这个问题的答案很多时候藏在代谢层面。传统的单细胞转录组分析能告诉你某个亚群高表达什么基因但基因表达和代谢通量之间并不是一一对应的关系。一个代谢通路里的酶基因表达上调了不代表这条通路真的在高速运转因为代谢反应还受到底物浓度、产物反馈抑制、酶活性的翻译后调控等多重因素影响。这就是为什么我们需要COMPASS这类专门做代谢通量估计的算法。COMPASS的全称是ComputationalOptimization ofMetabolicActivity fromSingle-cellSequencing翻译过来就是“从单细胞测序数据计算优化代谢活性”。它的核心思路是把单细胞转录组数据映射到一个基因组尺度的代谢网络模型上然后通过线性规划求解每个细胞在给定营养条件下各条代谢反应的通量分布。说白了它不是在看你表达了多少代谢酶而是在算你这些酶组合起来能跑出多大的代谢流量。这个思路对于免疫治疗研究特别有价值。比如你想知道肿瘤浸润T细胞的代谢状态是不是导致免疫治疗耐药的原因用COMPASS跑一遍就能看到不同T细胞亚群在糖酵解、氧化磷酸化、脂肪酸氧化、氨基酸代谢等通路上的通量差异。这些差异往往比单纯的基因表达差异更能解释功能表型。适合读这篇内容的人已经做过基础的单细胞分析、想往代谢方向拓展的科研人员正在做免疫治疗相关课题、需要找新靶点的研究者以及任何对单细胞数据挖掘有需求、想了解COMPASS这个工具的人。我会从算法原理讲到实操步骤再到差异分析和靶点筛选尽量把踩过的坑都标出来。2. COMPASS算法到底在算什么2.1 从基因表达到代谢通量的映射逻辑COMPASS的核心是一个基于约束的代谢模型。你可以把它想象成一个城市的交通网络道路就是代谢反应路口就是代谢物车流量就是通量。每条道路的通行能力是有限的对应酶的最大催化速率而且进入路口的车必须等于离开路口的车对应代谢物的质量平衡。具体来说COMPASS使用的是一个经过整理的基因组尺度代谢模型通常是人类代谢网络Recon系列模型。这个模型包含了数千个代谢反应和代谢物覆盖了中心碳代谢、氨基酸代谢、脂质代谢、核苷酸代谢等主要通路。算法的输入是单细胞的基因表达矩阵。它首先根据每个基因的表达水平来设定对应代谢反应的容量约束。这里有一个关键步骤COMPASS不是简单地用表达量乘以一个系数而是通过一个基因-反应映射规则GPRGene-Protein-Reaction association来确定。比如一个反应需要酶A和酶B同时存在才能催化那么它的容量就取决于A和B表达量的较小值如果酶A或酶B任意一个存在就能催化那就取较大值。这个映射过程决定了COMPASS对输入数据质量比较敏感。如果你的单细胞数据dropout率很高很多代谢酶基因检测不到那COMPASS估计出来的通量就会偏低甚至为零。所以实际操作中数据预处理和基因覆盖度检查是第一步要做的。2.2 线性规划求解与通量分布设定好容量约束之后COMPASS会为每个细胞求解一个线性规划问题。目标函数通常是最大化某个生物量反应的通量或者最大化整体代谢网络的活跃程度。约束条件包括每个代谢物的净通量为零稳态假设每条反应的通量不超过其容量上限热力学不可逆反应的通量方向固定求解这个线性规划问题就能得到每个细胞在代谢网络上的一个通量分布向量。这个向量的维度等于代谢反应的个数通常有几千维。直接拿这个高维向量做下游分析不太现实所以COMPASS还会计算一个代谢通路活性评分把通量按通路聚合得到每个细胞在各条代谢通路上的活性值。这里有一个实操中容易忽略的点COMPASS的求解是随机规划的也就是说同一个细胞跑两次可能得到略微不同的通量分布。这是因为线性规划的最优解可能不唯一COMPASS会在最优解空间里随机采样。所以如果你要做差异分析建议设置足够的采样次数通常1000次以上然后取平均值或者中位数作为该细胞的代谢活性估计。2.3 COMPASS与其他代谢分析工具的区别市面上做单细胞代谢分析的工具不止COMPASS一个还有scFEA、Compass注意大小写不同是另一个工具、MetabolicActivity等。它们之间的核心区别在于工具核心方法输出适用场景COMPASS约束优化线性规划反应水平通量分布需要精确通量估计、通路间交互分析scFEA神经网络代谢网络通路活性评分大规模数据、快速筛选MetabolicActivity基因集评分通路评分快速探索、初步筛选Compass通量平衡分析通量值微生物代谢、特定模型COMPASS的优势在于它的机制可解释性。因为它是基于真实的代谢网络约束来求解的所以得到的通量值有明确的生化意义。你可以追踪某个代谢物是从哪条通路来的、到哪条通路去的。这对于找靶点特别重要因为你要找的不是一个统计上显著的差异而是一个有生物学机制支撑的干预点。但COMPASS的劣势也很明显计算量大。一个包含几千个细胞的数据集跑COMPASS可能需要几个小时甚至更久取决于你设置的采样次数和并行核数。所以如果你的数据量很大建议先做细胞亚群抽样每个亚群取几百个细胞代表跑完COMPASS再映射回全体。3. 从原始数据到COMPASS输入完整实操流程3.1 数据准备与质控要点假设你已经有了一个经过基础处理的单细胞Seurat对象或者Scanpy的AnnData对象。在跑COMPASS之前有几件事必须确认第一基因注释版本要统一。COMPASS使用的代谢模型是基于特定版本的基因注释通常是Entrez ID或Ensembl ID。如果你的单细胞数据用的是基因symbol需要先做ID转换。我建议统一转成Entrez ID因为Recon模型就是用Entrez ID索引的。转换工具可以用org.Hs.eg.db人类或org.Mm.eg.db小鼠注意处理多对一和一对多的映射关系。第二检查代谢基因的覆盖度。人类代谢模型涉及大约2000-3000个代谢基因。你可以先算一下你的数据里检测到了多少个。如果覆盖率低于60%COMPASS的结果可靠性会打折扣。提高覆盖度的方法包括使用更宽松的质控阈值比如min.cells3而不是10、考虑使用imputation方法补dropout但要谨慎imputation可能引入假信号。第三营养条件设定。COMPASS需要你指定培养基的营养成分包括葡萄糖、谷氨酰胺、氧气等的浓度。这个参数直接影响通量分布。如果你做的是肿瘤微环境研究可以设置低葡萄糖、低氧的条件来模拟肿瘤内环境如果是体外培养实验就用标准培养基的浓度。这个参数没有绝对正确的值但一定要在文章方法里写清楚。# 检查代谢基因覆盖度示例R library(org.Hs.eg.db) metabolic_genes - read.csv(recon_metabolic_genes.csv) # 从模型文件提取 entrez_ids - mapIds(org.Hs.eg.db, keys rownames(seurat_obj), column ENTREZID, keytype SYMBOL) detected - sum(entrez_ids %in% metabolic_genes$EntrezID) coverage - detected / length(metabolic_genes$EntrezID) cat(代谢基因覆盖率:, round(coverage * 100, 1), %\n)3.2 构建COMPASS可识别的输入对象COMPASS的R包COMPASS接受的是一个特定格式的列表对象包含表达矩阵、基因ID映射、营养条件等。构建这个对象的过程有点繁琐但每一步都有明确的逻辑。首先你需要从单细胞对象中提取原始计数矩阵不是归一化后的数据。COMPASS内部会自己做归一化如果你提前归一化了反而会干扰它的容量估计。提取的时候注意保持基因和细胞的对应关系。然后你需要准备一个基因-反应映射表。这个表可以从Recon模型文件里提取包含每个反应对应的基因ID和GPR规则。COMPASS包通常会自带这个映射表但如果你用的是自定义模型就需要自己构建。# 构建COMPASS输入对象R伪代码 library(COMPASS) # 提取原始计数 counts_matrix - GetAssayData(seurat_obj, slot counts) # 准备反应-基因映射 reaction_gene_map - prepare_reaction_gene_map(model recon2) # 设定营养条件 nutrient_conditions - list( glucose 5, # mM glutamine 2, # mM oxygen 0.1, # 相对浓度低氧模拟肿瘤 serine 0.5, glycine 0.5 ) # 构建COMPASS对象 compass_input - COMPASSContainer( data counts_matrix, metadata seurat_objmeta.data, reaction_gene_map reaction_gene_map, nutrient_conditions nutrient_conditions )这里有一个容易踩的坑COMPASS对输入矩阵的稀疏性很敏感。如果你的数据里有很多零值COMPASS可能会把某些反应的容量设为零导致整条通路不通。解决办法是在构建对象之前先对表达矩阵做一次轻度的平滑处理比如用Seurat::SmoothData()或者简单的kNN平滑。但注意不要过度平滑否则会抹掉真实的生物学差异。3.3 运行COMPASS与参数调优COMPASS的核心函数是COMPASS()主要参数包括n_samples每个细胞的采样次数默认1000。建议至少1000如果计算资源允许可以设到5000。n_cores并行核数。COMPASS支持多核并行建议设为可用核数的80%。timeout单个细胞求解的超时时间。如果某个细胞特别难求解可以设一个上限避免卡死。verbose是否输出进度信息。第一次跑建议设为TRUE方便监控。# 运行COMPASS compass_result - COMPASS( data compass_input, n_samples 2000, n_cores 8, timeout 300, verbose TRUE )跑完之后你会得到一个包含每个细胞代谢通量分布的列表。这个结果对象比较大建议及时保存为RDS文件避免重复计算。saveRDS(compass_result, file compass_result_full.rds)参数调优的经验如果你发现很多细胞的求解状态是infeasible不可行说明你的营养条件设得太苛刻了或者基因覆盖度太低。可以尝试放宽营养条件比如提高葡萄糖浓度或者检查一下是不是有大量代谢基因没检测到。另一个常见问题是求解时间过长这时候可以适当降低n_samples或者先用一个子集测试参数确认没问题再跑全量。4. 差异分析与代谢靶点筛选实战4.1 从通量分布到通路活性评分COMPASS跑完之后你拿到的是每个细胞在每个代谢反应上的通量值。这个数据维度太高直接做差异分析容易假阳性。所以第一步是降维把反应水平的通量按代谢通路聚合得到通路水平的活性评分。聚合的方法有几种可以简单求和也可以取平均还可以用通量变异系数来加权。我通常用的是通路内反应通量的中位数因为中位数对异常值更稳健。聚合的映射关系可以从代谢模型文件里提取每个反应属于哪个子系统subsystem都有标注。# 通路活性聚合R伪代码 pathway_activity - aggregate_pathway_activity( compass_result compass_result, model recon2, method median ) # 结果是一个细胞×通路的矩阵 dim(pathway_activity) # [1] 5000 80 假设5000个细胞80条通路得到通路活性矩阵之后就可以做常规的差异分析了。如果你有分组信息比如响应组vs非响应组可以用Wilcoxon秩和检验或者limma来做差异通路分析。注意做多重检验校正FDR 0.05作为阈值。4.2 差异代谢通路的可视化与解读差异分析做完之后怎么展示结果也很关键。我常用的几种可视化方式第一种通路活性热图。把差异显著的代谢通路挑出来画一个细胞×通路的热图按分组排序。这样能直观看到哪些通路在哪个组里更活跃。第二种火山图。横轴是通路活性的差异倍数log2FC纵轴是显著性-log10 FDR。可以标注出最显著的那几个通路。第三种代谢网络图。把差异通路映射到代谢网络图上用颜色表示上调或下调。这种图能帮你看到通路之间的关联比如糖酵解上调的同时乳酸分泌通路是不是也上调了。第四种单细胞水平的通量分布图。比如你想看糖酵解通路在T细胞亚群间的差异可以画一个小提琴图或者箱线图每个点代表一个细胞的糖酵解活性。这里有一个解读上的注意事项COMPASS估计的是相对通量不是绝对通量。也就是说你只能比较不同细胞或不同组之间的相对高低不能说某个细胞的糖酵解通量是多少毫摩尔每秒。所以在写文章的时候措辞要准确用higher glycolytic flux而不是increased glycolytic rate by X-fold。4.3 从差异通路到可干预靶点的筛选逻辑找到差异代谢通路只是第一步真正的目标是找到可干预的靶点。这里的逻辑链条是差异通路 → 关键酶 → 可成药性 → 实验验证。具体来说你可以按以下步骤筛选第一步锁定差异通路中的关键酶。比如你发现氧化磷酸化通路在耐药组中显著下调那就去看这条通路里哪些酶的基因表达也同步下调。COMPASS的结果里其实包含了反应水平的通量你可以直接提取氧化磷酸化相关反应的通量值看哪些反应差异最大。第二步评估可成药性。不是所有代谢酶都有现成的小分子抑制剂。你可以查一下DrugBank、TTD等数据库看目标酶有没有已知的抑制剂或者正在临床试验的药物。优先选择那些已有工具化合物的靶点这样后续验证会快很多。第三步检查靶点的表达特异性。理想情况下你希望靶点在肿瘤细胞或特定免疫细胞亚群里高表达而在正常组织里低表达。可以用GTEx或者Human Protein Atlas的数据来交叉验证。第四步设计体外验证实验。常见的验证思路包括用siRNA或CRISPR敲低目标酶看是否逆转代谢表型用抑制剂处理细胞看是否影响免疫功能比如T细胞的杀伤能力、细胞因子分泌等做代谢流实验如Seahorse确认通量变化。我自己的经验是不要一上来就盯着最显著的那条通路。有时候最显著的通路是下游效应不是驱动因素。反而是一些中等显著但处于网络枢纽位置的通路更值得关注。你可以用通路间的相关性分析来识别枢纽通路如果一条通路和很多其他差异通路都高度相关那它可能是上游调控点。5. 常见问题与排查技巧实录5.1 COMPASS运行报错与解决方案问题一求解器报错infeasible problem。这是最常见的问题通常是因为营养条件设得太严格或者某些必需反应的基因表达为零。解决办法先检查营养条件把葡萄糖和谷氨酰胺浓度适当调高然后检查基因覆盖度如果某些关键代谢基因缺失可以考虑用同源基因替代或者放宽GPR规则。问题二运行时间过长。COMPASS的计算复杂度是O(n_cells × n_reactions × n_samples)。如果你的细胞数超过1万建议先做细胞抽样每个亚群取200-500个代表细胞。或者降低n_samples到500虽然精度会下降但对于初步探索够用了。问题三结果中大量细胞的通量为零。这通常是因为输入矩阵太稀疏。可以尝试用ALRA或SAVER做imputation但要注意imputation可能引入假阳性。另一个办法是聚合细胞把同一亚群的细胞表达谱平均一下当成一个元细胞来跑COMPASS这样能提高信噪比。问题四不同批次的数据结果不可比。COMPASS的结果受批次效应影响很大因为不同批次的测序深度和基因检测率不同。解决办法在跑COMPASS之前先做批次校正如Harmony、Seurat Integration或者把批次作为协变量纳入下游差异分析模型。5.2 差异分析中的统计陷阱陷阱一伪重复问题。如果你的每个处理组只有一个样本但有很多细胞直接把细胞当独立样本做检验会严重高估显著性。正确的做法是用混合效应模型把样本作为随机效应或者用伪批量pseudobulk方法先把细胞聚合成样本水平的均值再做检验。陷阱二多重检验校正方法选择。COMPASS输出的通路有几十到上百条做多重检验校正是必须的。但不同校正方法Bonferroni、BH、Storey的结果可能差异很大。我通常用BH方法因为它控制FDR的同时保留了较好的检验效能。如果通路数很少20可以用Bonferroni更严格。陷阱三忽略代谢通路的层级结构。代谢通路不是独立的糖酵解和氧化磷酸化通过丙酮酸脱氢酶连接氨基酸代谢和核苷酸代谢通过一碳单位连接。做差异分析时如果只盯着单条通路可能会错过通路间的协同变化。建议做一次通路共表达网络分析看看哪些通路模块在组间差异最大。5.3 从COMPASS结果到实验验证的衔接COMPASS给的是计算预测最终还是要实验验证。这里分享几个衔接上的经验经验一优先验证通量差异最大的反应而不是通路。通路水平的聚合会平滑掉一些细节有时候通路整体不显著但通路里的某条关键反应显著。你可以直接从COMPASS结果里提取反应水平的通量做差异分析往往能发现更有意思的靶点。经验二用代谢物检测来交叉验证。如果COMPASS预测糖酵解通量上调你可以用质谱检测乳酸分泌量来验证如果预测脂肪酸氧化下调可以检测酰基肉碱水平。代谢物数据比转录组更接近表型验证价值更高。经验三注意COMPASS的预测方向可能和基因表达相反。这听起来很奇怪但确实会发生。比如某个代谢酶的基因表达上调了但COMPASS预测的通量反而下降可能是因为底物供应不足或者产物反馈抑制。这种情况反而更有意思说明存在翻译后调控或者代谢物介导的调控值得深入挖掘。经验四不要忽视低通量但高特异性的通路。有些代谢通路比如色氨酸-犬尿氨酸通路在整体通量中占比很小但在免疫调节中作用关键。COMPASS的通量估计对这些低丰度通路可能不够敏感你可以单独提取相关反应的通量值做更细致的分析。6. 一些实操中的个人体会COMPASS这个工具我从两年前开始用踩过的坑不算少。最大的一个体会是不要把它当成一个黑箱。很多人跑完COMPASS拿到通路活性矩阵就直接扔进差异分析流程然后对着结果硬解释。这样很容易得出似是而非的结论。我的做法是每次跑完COMPASS都会先做一轮质量检查看看通量分布是否合理比如糖酵解通量在活跃增殖的细胞里应该偏高、看看有没有大量零值、看看不同批次的分布是否一致。这些检查花不了多少时间但能避免很多后续的麻烦。另一个体会是COMPASS的结果要和基因表达数据结合看。如果COMPASS预测某条通路通量上调但通路里的酶基因表达没变化那可能是酶活性调控如果酶基因表达也上调那就是转录调控。这两种情况的干预策略完全不同前者可能需要找激酶或者变构调节剂后者可以直接用转录抑制剂或者RNAi。还有一点营养条件的设定对结果影响很大。我试过同一批数据用高糖和低糖条件跑COMPASS得到的差异通路列表重叠度只有60%左右。所以如果你的研究涉及肿瘤微环境或者免疫代谢一定要根据实际生理条件来设定营养参数不要用默认值。最后关于计算资源如果你的数据量很大建议用服务器跑至少32GB内存起步。COMPASS的R包对内存管理不是特别优化细胞数超过2万的时候可能会爆内存。一个变通的办法是分批次跑每次跑5000个细胞然后把结果合并。这个方向后续还可以往空间代谢组结合的方向扩展。现在有空间转录组数据了如果能结合COMPASS做空间代谢通量估计就能看到代谢异质性在组织空间上的分布这对于理解肿瘤微环境里的免疫代谢互作会很有帮助。不过这是另一个大话题了有机会再展开聊。
