简介：这份源码资源面向从事单细胞转录组与代谢研究的科研人员及生物信息学初学者，围绕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_obj@assays$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, solver='cplex'): if solver == 'cplex': import cplex problem = cplex.Cplex() # 设置目标函数、约束... elif solver == 'scipy': from scipy.optimize import linprog # 用 scipy 的 linprog 替代 res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_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, pathways=key_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_obj@assays$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_components=2) flux_pca = pca.fit_transform(flux_matrix.T) # 和 UMAP 嵌入做相关性 from scipy.stats import spearmanr corr, pval = spearmanr(flux_pca[:, 0], umap_embedding[:, 0]) print(f"Correlation: {corr:.3f}, p-value: {pval:.2e}")如果相关性显著，说明代谢通量的变化和细胞状态的变化是一致的。如果完全不相关，要么是代谢分析出了问题，要么是这个细胞群体的代谢异质性本身就很低。
我自己的习惯是：每次跑完代谢分析，先看 UMAP 映射图，再看几个关键通路的箱线图，最后和已知的生物学知识对一遍。如果糖酵解在增殖细胞里不高，TCA 在静息细胞里不低，那大概率是哪里出了问题。这个验证流程帮我省了很多后悔药。希望帮到你。
本文还有配套的精品资源，点击获取
