尧图网页设计 
简介这份源码资源面向从事单细胞转录组与代谢研究的科研人员及生物信息学初学者围绕scMetabolism包展开小鼠单细胞代谢激活分数分析重点解决小鼠基因名向人类基因名转换、以及适配Seurat v4/v5版本进行代谢通路打分的问题。资源包共6个文件以R脚本为主包含代谢分析主流程脚本与依赖包安装脚本另附README说明文档、HTML页面及项目配置文件压缩包约9KB体量轻便便于快速部署与二次修改。目前已有176人学习下载。读者可从中获得从基因名转换到代谢激活分数计算的完整代码示例掌握将数据导入Seurat并完成单细胞层面代谢特征解读的方法同时借助参考链接与说明文档理解分析思路适合作为单细胞代谢研究的入门模板与排错参考。1. 小鼠单细胞代谢分析源码从表达矩阵到代谢通路的可复现路径单细胞转录组测序做完之后绝大多数人停在细胞分群和标记基因注释这一步真正往代谢方向挖的人不多。原因很直接代谢分析不像差异表达那样有现成的一键流程它需要把每个细胞的表达谱映射到代谢反应网络再算通量或打分中间涉及基因ID转换、反应-基因对应关系、细胞亚群聚合等多个环节。小鼠数据又比人数据多一层麻烦——基因命名规则不同大量基因以Gm开头同源基因对应关系需要额外处理。这篇要讲的就是围绕「小鼠单细胞代谢分析源码」这条线把从原始表达矩阵到代谢通路活性打分的完整链路拆开。适合已经跑过 Seurat 或 Scanpy 基础流程、手里有小鼠单细胞数据、想往代谢方向延伸的从业者。核心工具是 scMetabolism 和 Compass 两条路线前者基于 KEGG 反应集做打分后者用约束优化算实际通量。两条路线的源码结构、参数含义、小鼠适配的坑都会落到可执行的代码层面。2. scMetabolism 源码拆解VISION 与 AUCell 两套打分引擎怎么选2.1 源码目录结构与核心函数入口scMetabolism 的源码组织比较紧凑核心逻辑集中在R/目录下。拿到源码包后先看三个文件scMetabolism.R是主函数入口utility.R负责基因ID转换和矩阵预处理get_metabolism_data.R管理 KEGG 反应集的加载。主函数sc.metabolism.Seurat()和sc.metabolism.Seurat()分别对应 Seurat 对象和普通矩阵输入。源码里最关键的设计是它不直接算代谢物浓度而是把每个 KEGG 反应关联的基因集当作一个「签名」用单细胞表达数据对这个签名打分。打分方法有两套——VISION 和 AUCell。VISION 的做法是把基因集得分投影到细胞嵌入空间适合看整体趋势AUCell 则基于排名对稀疏数据更稳健。# 加载源码包假设已 clone 到本地 devtools::load_all(/path/to/scMetabolism) # 查看主函数参数 args(sc.metabolism.Seurat) # function(obj, method AUCell, imputation FALSE, # ncores 2, metabolism.type KEGG, ...)这里method控制打分引擎imputation控制是否对 dropout 做插补metabolism.type目前支持 KEGG 和 REACTOME。ncores在 AUCell 模式下影响并行计算VISION 模式下基本用不到。2.2 小鼠基因ID转换的源码逻辑与实操小鼠数据的第一个坑就在基因ID转换。scMetabolism 内置的 KEGG 反应集是基于人类基因符号HGNC构建的直接拿小鼠的Gm基因去匹配命中率会低得离谱。源码里utility.R有一个convertHumanGeneList()函数但它是给人数据用的。小鼠数据需要先做同源转换。常见做法是用 biomaRt 或 homologene 包把小鼠基因符号转成人类同源基因符号。我一般用 homologene因为它不依赖网络速度快。library(homologene) library(Seurat) # 假设 seurat_obj 是小鼠数据 mouse_genes - rownames(seurat_obj) # homologene 转换小鼠 - 人类 human_genes - homologene(mouse_genes, inTax 10090, outTax 9606) # 去重一个小鼠基因可能对应多个人类基因取第一个 human_genes - human_genes[!duplicated(human_genes$10090), ] # 建立映射表 gene_map - setNames(human_genes$9606, human_genes$10090) # 把 Seurat 对象的基因名替换为人类符号 # 注意只保留能转换的基因 valid_genes - intersect(mouse_genes, names(gene_map)) seurat_obj - subset(seurat_obj, features valid_genes) # 重命名 new_names - gene_map[rownames(seurat_obj)] # 处理重复如果多个人类基因对应同一个小鼠基因保留表达量最高的 # 这里简化处理直接去重 seurat_obj - seurat_obj[!duplicated(new_names), ] rownames(seurat_obj) - new_names[!duplicated(new_names)]这段代码的逻辑是先建立小鼠到人类的同源映射然后只保留能映射的基因最后用人类符号替换。参数inTax 10090是小鼠的 NCBI 分类号outTax 9606是人类。去重那一步不能省否则后续打分函数会因为重复基因名报错。转换完成后命中率能从不到 30% 提升到 70% 以上。如果还是偏低检查一下数据里是不是有大量Gm基因——这些基因很多没有人类同源属于正常丢失。2.3 AUCell 打分参数调优与结果解读AUCell 的核心参数是aucMaxRank默认是基因集大小的 5%。这个参数控制排名阈值值越小越严格只取表达排名最靠前的基因。对于小鼠数据我一般会把它调到 10%因为同源转换后基因集覆盖度下降需要放宽阈值来补偿。# 运行代谢分析 seurat_obj - sc.metabolism.Seurat( obj seurat_obj, method AUCell, imputation FALSE, ncores 4, metabolism.type KEGG ) # 提取结果 metabolism_matrix - seurat_objassays$METABOLISM$score # 查看前 5 个通路在部分细胞中的得分 metabolism_matrix[1:5, 1:3]imputation FALSE是默认值对于 10x 数据插补会引入假信号不建议开。ncores根据机器配置调整AUCell 的并行效率不错4 核比单核快 2 倍左右。结果矩阵的行是 KEGG 通路列是细胞。得分范围在 0 到 1 之间越高表示该通路在细胞中越活跃。解读时不要只看绝对值要看相对差异——比如比较肿瘤细胞和正常细胞的糖酵解通路得分差异倍数比绝对得分更有意义。3. Compass 源码路线约束优化算代谢通量的落地细节3.1 Compass 的数学模型与源码依赖Compass 走的是另一条路。它不满足于「打分」而是用约束优化linear programming来估算每个细胞在代谢网络中的实际通量分布。源码核心在compass/目录下compass.py是主入口reactions.py定义反应网络solver.py封装了求解器调用。数学模型大致是给定一个细胞的基因表达谱先映射到酶活性再构建一个代谢网络最后用线性规划求解在满足稳态约束下各反应的通量。目标函数是最大化通量总和约束包括质量平衡、热力学可行性等。依赖方面Compass 需要cplex或gurobi求解器。开源方案可以用scipy.optimize.linprog但速度慢很多。源码里默认走 cplex如果没有 license需要改solver.py里的求解器配置。# compass 源码中 solver.py 的关键片段 def solve_flux(expression_matrix, model, solvercplex): if solver cplex: import cplex problem cplex.Cplex() # 设置目标函数、约束... elif solver scipy: from scipy.optimize import linprog # 用 scipy 的 linprog 替代 res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq) return res如果要用 scipy 替代需要把 cplex 的 API 调用全部重写。我一般建议直接用 cplex 的社区版学术用途免费安装也不复杂。3.2 小鼠代谢网络构建与反应-基因映射Compass 默认使用 Recon 系列代谢网络模型。Recon3D 是人类模型小鼠需要用 Recon 的小鼠版本或者做同源映射。源码里reactions.py有一个load_model()函数可以指定模型文件路径。from compass.reactions import load_model # 加载小鼠代谢模型 # 常见做法是用 Recon3D 做同源映射或者直接用 Mouse Recon model load_model(/path/to/mouse_recon.json) # 查看模型中的反应数量 print(len(model.reactions)) # 通常 10000 个反应小鼠模型的文件格式一般是 JSON 或 SBML。如果没有现成的小鼠模型可以用recon3d加同源基因映射来构建。这一步比较耗时但一次构建后可以复用。反应-基因映射是另一个关键点。Compass 需要知道每个反应由哪些基因编码的酶催化。源码里reactions.py的map_genes_to_reactions()函数负责这件事。小鼠数据同样需要先做基因符号转换逻辑和 scMetabolism 那边一样。3.3 通量结果的后处理与可视化Compass 跑完之后输出是一个通量矩阵行是反应列是细胞。这个矩阵非常稀疏直接看很难看出模式。源码里提供了一些后处理函数比如compass.visualization模块下的plot_flux()。import compass from compass.visualization import plot_flux import pandas as pd # 假设 flux_matrix 是 Compass 输出的通量矩阵 # 按细胞类型聚合 cell_types seurat_obj.obs[cell_type] flux_by_type flux_matrix.groupby(cell_types).mean() # 挑选几个关键通路 key_pathways [glycolysis, TCA, oxidative_phosphorylation] # 需要把反应映射到通路这一步可以用 KEGG 或 Reactome 的注释 plot_flux(flux_by_type, pathwayskey_pathways)后处理的核心是降维和聚合。单细胞通量矩阵维度太高直接可视化不现实。常见做法是先按细胞类型求平均再挑关键通路画热图或箱线图。参数方面groupby的粒度可以按聚类结果也可以按已知的细胞类型注释。4. 避坑与排查小鼠单细胞代谢分析里最容易翻车的五个点4.1 基因转换后命中率过低现象跑完 scMetabolism 或 Compass发现大部分通路的得分都是 0 或者接近 0检查发现基因集覆盖度不到 20%。原因小鼠基因符号没有正确转换或者转换后没有去重导致大量基因被丢弃。另一个常见原因是数据里本身就有很多Gm基因这些基因没有人类同源。解决先用 homologene 做转换检查转换率。如果低于 50%看看是不是用了错误的分类号。小鼠是 10090人类是 9606别搞反。转换后去重时保留表达量最高的那个同源基因而不是随机取一个。4.2 AUCell 打分全为 1 或全为 0现象AUCell 输出的得分矩阵里所有值都是 1 或者都是 0没有中间值。原因aucMaxRank设置不当。如果设得太小比如 1%而基因集又很大排名阈值会覆盖几乎所有基因导致得分饱和。反过来如果设得太大比如 50%得分会趋近于 0。解决把aucMaxRank调到基因集大小的 5% 到 10% 之间。具体值可以通过试跑几个通路来校准。另外检查一下输入矩阵是不是已经做了归一化AUCell 对原始 counts 和归一化数据的表现不同。4.3 Compass 求解器报错或超时现象Compass 跑到一半报错提示 solver 不可用或者跑了几小时还没结束。原因cplex 没有正确安装或 license 过期。另一个原因是细胞数量太多线性规划的计算量随细胞数线性增长。解决先确认 cplex 能正常 import。如果不行换 gurobi 或者用 scipy 的 linprog 做小规模测试。对于大规模数据建议先做细胞降采样比如每个 cluster 随机抽 200 个细胞跑完再映射回去。4.4 代谢通路得分与生物学预期不符现象明明知道某个细胞类型应该高糖酵解但得分却很低。原因可能是通路定义的问题。KEGG 的糖酵解通路包含的基因和实际糖酵解酶有出入或者小鼠的同源基因映射丢失了关键酶。解决手动检查关键基因是否在基因集里。比如糖酵解的关键酶Hk2、Pkm、Ldha看看它们有没有被正确转换和保留。如果丢了考虑用 REACTOME 通路集替代或者手动补充基因集。4.5 结果不可复现现象同样的数据两次跑出来的结果不一样。原因AUCell 的并行计算有随机性或者 Compass 的求解器有数值精度问题。另一个常见原因是基因转换时用了不同的数据库版本。解决在代码开头设set.seed(42)AUCell 的并行部分也要固定种子。Compass 那边确保求解器参数一致比如 tolerance 设成 1e-6。基因转换的数据库版本要记录在案homologene 的版本更新会导致映射结果变化。5. 进阶技巧把代谢通量映射回细胞嵌入空间做可视化验证跑完代谢分析拿到通量矩阵或得分矩阵之后怎么验证结果靠不靠谱我一般会做一件事把代谢得分映射回 UMAP 或 tSNE 嵌入空间看高得分细胞是不是聚集在特定的区域。这个操作在 scMetabolism 的源码里有现成的函数但 Compass 那边需要自己写。# scMetabolism 结果映射回 UMAP library(Seurat) library(ggplot2) # 假设 seurat_obj 已经跑完 sc.metabolism.Seurat # 提取糖酵解通路的得分 glycolysis_score - seurat_objassays$METABOLISM$score[Glycolysis, ] # 加到 metadata seurat_obj$glycolysis - glycolysis_score # 画 UMAP FeaturePlot(seurat_obj, features glycolysis, cols c(lightgrey, red), min.cutoff 0, max.cutoff 1)这段代码的逻辑很简单把通路得分作为一个 feature 加到 Seurat 对象里然后用FeaturePlot画出来。参数min.cutoff和max.cutoff控制颜色映射范围设成 0 到 1 可以避免极端值影响视觉效果。对于 Compass 的通量结果需要先做降维。常见做法是用 PCA 或 NMF 把通量矩阵降到 2 维再和 UMAP 做相关性分析。如果代谢通量的主成分和细胞类型的分布高度相关说明结果可信。from sklearn.decomposition import PCA import numpy as np # flux_matrix 是 Compass 输出 pca PCA(n_components2) flux_pca pca.fit_transform(flux_matrix.T) # 和 UMAP 嵌入做相关性 from scipy.stats import spearmanr corr, pval spearmanr(flux_pca[:, 0], umap_embedding[:, 0]) print(fCorrelation: {corr:.3f}, p-value: {pval:.2e})如果相关性显著说明代谢通量的变化和细胞状态的变化是一致的。如果完全不相关要么是代谢分析出了问题要么是这个细胞群体的代谢异质性本身就很低。我自己的习惯是每次跑完代谢分析先看 UMAP 映射图再看几个关键通路的箱线图最后和已知的生物学知识对一遍。如果糖酵解在增殖细胞里不高TCA 在静息细胞里不低那大概率是哪里出了问题。这个验证流程帮我省了很多后悔药。希望帮到你。本文还有配套的精品资源点击获取 小鼠单细胞代谢分析源码解析：从表达矩阵到代谢通路的可复现路径
简介这份源码资源面向从事单细胞转录组与代谢研究的科研人员及生物信息学初学者围绕scMetabolism包展开小鼠单细胞代谢激活分数分析重点解决小鼠基因名向人类基因名转换、以及适配Seurat v4/v5版本进行代谢通路打分的问题。资源包共6个文件以R脚本为主包含代谢分析主流程脚本与依赖包安装脚本另附README说明文档、HTML页面及项目配置文件压缩包约9KB体量轻便便于快速部署与二次修改。目前已有176人学习下载。读者可从中获得从基因名转换到代谢激活分数计算的完整代码示例掌握将数据导入Seurat并完成单细胞层面代谢特征解读的方法同时借助参考链接与说明文档理解分析思路适合作为单细胞代谢研究的入门模板与排错参考。1. 小鼠单细胞代谢分析源码从表达矩阵到代谢通路的可复现路径单细胞转录组测序做完之后绝大多数人停在细胞分群和标记基因注释这一步真正往代谢方向挖的人不多。原因很直接代谢分析不像差异表达那样有现成的一键流程它需要把每个细胞的表达谱映射到代谢反应网络再算通量或打分中间涉及基因ID转换、反应-基因对应关系、细胞亚群聚合等多个环节。小鼠数据又比人数据多一层麻烦——基因命名规则不同大量基因以Gm开头同源基因对应关系需要额外处理。这篇要讲的就是围绕「小鼠单细胞代谢分析源码」这条线把从原始表达矩阵到代谢通路活性打分的完整链路拆开。适合已经跑过 Seurat 或 Scanpy 基础流程、手里有小鼠单细胞数据、想往代谢方向延伸的从业者。核心工具是 scMetabolism 和 Compass 两条路线前者基于 KEGG 反应集做打分后者用约束优化算实际通量。两条路线的源码结构、参数含义、小鼠适配的坑都会落到可执行的代码层面。2. scMetabolism 源码拆解VISION 与 AUCell 两套打分引擎怎么选2.1 源码目录结构与核心函数入口scMetabolism 的源码组织比较紧凑核心逻辑集中在R/目录下。拿到源码包后先看三个文件scMetabolism.R是主函数入口utility.R负责基因ID转换和矩阵预处理get_metabolism_data.R管理 KEGG 反应集的加载。主函数sc.metabolism.Seurat()和sc.metabolism.Seurat()分别对应 Seurat 对象和普通矩阵输入。源码里最关键的设计是它不直接算代谢物浓度而是把每个 KEGG 反应关联的基因集当作一个「签名」用单细胞表达数据对这个签名打分。打分方法有两套——VISION 和 AUCell。VISION 的做法是把基因集得分投影到细胞嵌入空间适合看整体趋势AUCell 则基于排名对稀疏数据更稳健。# 加载源码包假设已 clone 到本地 devtools::load_all(/path/to/scMetabolism) # 查看主函数参数 args(sc.metabolism.Seurat) # function(obj, method AUCell, imputation FALSE, # ncores 2, metabolism.type KEGG, ...)这里method控制打分引擎imputation控制是否对 dropout 做插补metabolism.type目前支持 KEGG 和 REACTOME。ncores在 AUCell 模式下影响并行计算VISION 模式下基本用不到。2.2 小鼠基因ID转换的源码逻辑与实操小鼠数据的第一个坑就在基因ID转换。scMetabolism 内置的 KEGG 反应集是基于人类基因符号HGNC构建的直接拿小鼠的Gm基因去匹配命中率会低得离谱。源码里utility.R有一个convertHumanGeneList()函数但它是给人数据用的。小鼠数据需要先做同源转换。常见做法是用 biomaRt 或 homologene 包把小鼠基因符号转成人类同源基因符号。我一般用 homologene因为它不依赖网络速度快。library(homologene) library(Seurat) # 假设 seurat_obj 是小鼠数据 mouse_genes - rownames(seurat_obj) # homologene 转换小鼠 - 人类 human_genes - homologene(mouse_genes, inTax 10090, outTax 9606) # 去重一个小鼠基因可能对应多个人类基因取第一个 human_genes - human_genes[!duplicated(human_genes$10090), ] # 建立映射表 gene_map - setNames(human_genes$9606, human_genes$10090) # 把 Seurat 对象的基因名替换为人类符号 # 注意只保留能转换的基因 valid_genes - intersect(mouse_genes, names(gene_map)) seurat_obj - subset(seurat_obj, features valid_genes) # 重命名 new_names - gene_map[rownames(seurat_obj)] # 处理重复如果多个人类基因对应同一个小鼠基因保留表达量最高的 # 这里简化处理直接去重 seurat_obj - seurat_obj[!duplicated(new_names), ] rownames(seurat_obj) - new_names[!duplicated(new_names)]这段代码的逻辑是先建立小鼠到人类的同源映射然后只保留能映射的基因最后用人类符号替换。参数inTax 10090是小鼠的 NCBI 分类号outTax 9606是人类。去重那一步不能省否则后续打分函数会因为重复基因名报错。转换完成后命中率能从不到 30% 提升到 70% 以上。如果还是偏低检查一下数据里是不是有大量Gm基因——这些基因很多没有人类同源属于正常丢失。2.3 AUCell 打分参数调优与结果解读AUCell 的核心参数是aucMaxRank默认是基因集大小的 5%。这个参数控制排名阈值值越小越严格只取表达排名最靠前的基因。对于小鼠数据我一般会把它调到 10%因为同源转换后基因集覆盖度下降需要放宽阈值来补偿。# 运行代谢分析 seurat_obj - sc.metabolism.Seurat( obj seurat_obj, method AUCell, imputation FALSE, ncores 4, metabolism.type KEGG ) # 提取结果 metabolism_matrix - seurat_objassays$METABOLISM$score # 查看前 5 个通路在部分细胞中的得分 metabolism_matrix[1:5, 1:3]imputation FALSE是默认值对于 10x 数据插补会引入假信号不建议开。ncores根据机器配置调整AUCell 的并行效率不错4 核比单核快 2 倍左右。结果矩阵的行是 KEGG 通路列是细胞。得分范围在 0 到 1 之间越高表示该通路在细胞中越活跃。解读时不要只看绝对值要看相对差异——比如比较肿瘤细胞和正常细胞的糖酵解通路得分差异倍数比绝对得分更有意义。3. Compass 源码路线约束优化算代谢通量的落地细节3.1 Compass 的数学模型与源码依赖Compass 走的是另一条路。它不满足于「打分」而是用约束优化linear programming来估算每个细胞在代谢网络中的实际通量分布。源码核心在compass/目录下compass.py是主入口reactions.py定义反应网络solver.py封装了求解器调用。数学模型大致是给定一个细胞的基因表达谱先映射到酶活性再构建一个代谢网络最后用线性规划求解在满足稳态约束下各反应的通量。目标函数是最大化通量总和约束包括质量平衡、热力学可行性等。依赖方面Compass 需要cplex或gurobi求解器。开源方案可以用scipy.optimize.linprog但速度慢很多。源码里默认走 cplex如果没有 license需要改solver.py里的求解器配置。# compass 源码中 solver.py 的关键片段 def solve_flux(expression_matrix, model, solvercplex): if solver cplex: import cplex problem cplex.Cplex() # 设置目标函数、约束... elif solver scipy: from scipy.optimize import linprog # 用 scipy 的 linprog 替代 res linprog(c, A_ubA_ub, b_ubb_ub, A_eqA_eq, b_eqb_eq) return res如果要用 scipy 替代需要把 cplex 的 API 调用全部重写。我一般建议直接用 cplex 的社区版学术用途免费安装也不复杂。3.2 小鼠代谢网络构建与反应-基因映射Compass 默认使用 Recon 系列代谢网络模型。Recon3D 是人类模型小鼠需要用 Recon 的小鼠版本或者做同源映射。源码里reactions.py有一个load_model()函数可以指定模型文件路径。from compass.reactions import load_model # 加载小鼠代谢模型 # 常见做法是用 Recon3D 做同源映射或者直接用 Mouse Recon model load_model(/path/to/mouse_recon.json) # 查看模型中的反应数量 print(len(model.reactions)) # 通常 10000 个反应小鼠模型的文件格式一般是 JSON 或 SBML。如果没有现成的小鼠模型可以用recon3d加同源基因映射来构建。这一步比较耗时但一次构建后可以复用。反应-基因映射是另一个关键点。Compass 需要知道每个反应由哪些基因编码的酶催化。源码里reactions.py的map_genes_to_reactions()函数负责这件事。小鼠数据同样需要先做基因符号转换逻辑和 scMetabolism 那边一样。3.3 通量结果的后处理与可视化Compass 跑完之后输出是一个通量矩阵行是反应列是细胞。这个矩阵非常稀疏直接看很难看出模式。源码里提供了一些后处理函数比如compass.visualization模块下的plot_flux()。import compass from compass.visualization import plot_flux import pandas as pd # 假设 flux_matrix 是 Compass 输出的通量矩阵 # 按细胞类型聚合 cell_types seurat_obj.obs[cell_type] flux_by_type flux_matrix.groupby(cell_types).mean() # 挑选几个关键通路 key_pathways [glycolysis, TCA, oxidative_phosphorylation] # 需要把反应映射到通路这一步可以用 KEGG 或 Reactome 的注释 plot_flux(flux_by_type, pathwayskey_pathways)后处理的核心是降维和聚合。单细胞通量矩阵维度太高直接可视化不现实。常见做法是先按细胞类型求平均再挑关键通路画热图或箱线图。参数方面groupby的粒度可以按聚类结果也可以按已知的细胞类型注释。4. 避坑与排查小鼠单细胞代谢分析里最容易翻车的五个点4.1 基因转换后命中率过低现象跑完 scMetabolism 或 Compass发现大部分通路的得分都是 0 或者接近 0检查发现基因集覆盖度不到 20%。原因小鼠基因符号没有正确转换或者转换后没有去重导致大量基因被丢弃。另一个常见原因是数据里本身就有很多Gm基因这些基因没有人类同源。解决先用 homologene 做转换检查转换率。如果低于 50%看看是不是用了错误的分类号。小鼠是 10090人类是 9606别搞反。转换后去重时保留表达量最高的那个同源基因而不是随机取一个。4.2 AUCell 打分全为 1 或全为 0现象AUCell 输出的得分矩阵里所有值都是 1 或者都是 0没有中间值。原因aucMaxRank设置不当。如果设得太小比如 1%而基因集又很大排名阈值会覆盖几乎所有基因导致得分饱和。反过来如果设得太大比如 50%得分会趋近于 0。解决把aucMaxRank调到基因集大小的 5% 到 10% 之间。具体值可以通过试跑几个通路来校准。另外检查一下输入矩阵是不是已经做了归一化AUCell 对原始 counts 和归一化数据的表现不同。4.3 Compass 求解器报错或超时现象Compass 跑到一半报错提示 solver 不可用或者跑了几小时还没结束。原因cplex 没有正确安装或 license 过期。另一个原因是细胞数量太多线性规划的计算量随细胞数线性增长。解决先确认 cplex 能正常 import。如果不行换 gurobi 或者用 scipy 的 linprog 做小规模测试。对于大规模数据建议先做细胞降采样比如每个 cluster 随机抽 200 个细胞跑完再映射回去。4.4 代谢通路得分与生物学预期不符现象明明知道某个细胞类型应该高糖酵解但得分却很低。原因可能是通路定义的问题。KEGG 的糖酵解通路包含的基因和实际糖酵解酶有出入或者小鼠的同源基因映射丢失了关键酶。解决手动检查关键基因是否在基因集里。比如糖酵解的关键酶Hk2、Pkm、Ldha看看它们有没有被正确转换和保留。如果丢了考虑用 REACTOME 通路集替代或者手动补充基因集。4.5 结果不可复现现象同样的数据两次跑出来的结果不一样。原因AUCell 的并行计算有随机性或者 Compass 的求解器有数值精度问题。另一个常见原因是基因转换时用了不同的数据库版本。解决在代码开头设set.seed(42)AUCell 的并行部分也要固定种子。Compass 那边确保求解器参数一致比如 tolerance 设成 1e-6。基因转换的数据库版本要记录在案homologene 的版本更新会导致映射结果变化。5. 进阶技巧把代谢通量映射回细胞嵌入空间做可视化验证跑完代谢分析拿到通量矩阵或得分矩阵之后怎么验证结果靠不靠谱我一般会做一件事把代谢得分映射回 UMAP 或 tSNE 嵌入空间看高得分细胞是不是聚集在特定的区域。这个操作在 scMetabolism 的源码里有现成的函数但 Compass 那边需要自己写。# scMetabolism 结果映射回 UMAP library(Seurat) library(ggplot2) # 假设 seurat_obj 已经跑完 sc.metabolism.Seurat # 提取糖酵解通路的得分 glycolysis_score - seurat_objassays$METABOLISM$score[Glycolysis, ] # 加到 metadata seurat_obj$glycolysis - glycolysis_score # 画 UMAP FeaturePlot(seurat_obj, features glycolysis, cols c(lightgrey, red), min.cutoff 0, max.cutoff 1)这段代码的逻辑很简单把通路得分作为一个 feature 加到 Seurat 对象里然后用FeaturePlot画出来。参数min.cutoff和max.cutoff控制颜色映射范围设成 0 到 1 可以避免极端值影响视觉效果。对于 Compass 的通量结果需要先做降维。常见做法是用 PCA 或 NMF 把通量矩阵降到 2 维再和 UMAP 做相关性分析。如果代谢通量的主成分和细胞类型的分布高度相关说明结果可信。from sklearn.decomposition import PCA import numpy as np # flux_matrix 是 Compass 输出 pca PCA(n_components2) flux_pca pca.fit_transform(flux_matrix.T) # 和 UMAP 嵌入做相关性 from scipy.stats import spearmanr corr, pval spearmanr(flux_pca[:, 0], umap_embedding[:, 0]) print(fCorrelation: {corr:.3f}, p-value: {pval:.2e})如果相关性显著说明代谢通量的变化和细胞状态的变化是一致的。如果完全不相关要么是代谢分析出了问题要么是这个细胞群体的代谢异质性本身就很低。我自己的习惯是每次跑完代谢分析先看 UMAP 映射图再看几个关键通路的箱线图最后和已知的生物学知识对一遍。如果糖酵解在增殖细胞里不高TCA 在静息细胞里不低那大概率是哪里出了问题。这个验证流程帮我省了很多后悔药。希望帮到你。本文还有配套的精品资源点击获取