简介：面向从事单细胞转录组和代谢研究的科研人员、生物信息学初学者，这份源码包聚焦小鼠单细胞代谢激活分数分析。其核心思路基于scMetabolism工具包，涵盖从小鼠基因名到人类基因名的准确转换，同时适配Seurat的第四、第五版本，方便融入既有单细胞数据处理流程。压缩包共六个文件，包含两个R脚本（主分析脚本与依赖安装脚本）、一个HTML说明页、一个Markdown文档及配套配置文件，整体体积约九KB，结构简明、便于按需查阅。目前已有180人学习，适合希望快速掌握单细胞代谢评分实现细节、节省前期摸索时间的读者。借助包内代码与说明文档，可理解从表达矩阵到代谢特征评分的完整流程，并迁移应用到自身数据集，对深入认识不同细胞类型的代谢差异具有直接参考价值。

1. 单细胞代谢分析源码包：把表达矩阵变成代谢通量图

拿到一份小鼠组织的 10x 单细胞表达矩阵时，很多人第一反应是看几个代谢酶基因（比如 GAPDH、LDHA）的表达量。但单细胞代谢分析源码包做的事情要深一层：它把每个细胞的转录组信号映射到代谢通路上，算通路活性，再用代谢网络模型推断通量方向，最后把"细胞状态"和"代谢偏好"绑定。这份源码包适合已经会用 Seurat 或 Scanpy 处理基础数据的从业者，以及要做肿瘤微环境、免疫代谢、发育代谢这类课题的研究者。它解决的核心问题就一个：表达矩阵里有几千个代谢基因，怎么把它们整理成能解释生物学差异的代谢图景。

2. 核心原理与选型：代谢分析不能只看基因表达量

2.1 代谢分析的三层逻辑：表达量、通路活性、通量方向

常见误区是拿单个基因表达高低当代谢状态。代谢系统是网络结构，一个通路里几十个酶共同决定代谢产出，单基因的波动会被上游底物、下游产物和反馈抑制抵消。所以源码包把分析拆成三层：

第一层是表达矩阵本身。这一步只做质控和标准化，不涉及代谢判断。这个源码包以 Scanpy 为主流程，输入是 10x 的 h5 文件或 mtx 目录。

第二层是通路活性打分。做法是把 KEGG 或 Reactome 的代谢通路基因集拿出来，看每个细胞在这些基因上的整体表达趋势。常用方法有 GSVA、AUCell，这个包里也兼容 scMetabolism 的打分逻辑。打分结果是一个细胞 × 通路的矩阵，每个数值代表这条通路在该细胞里的"活跃程度"。

第三层是代谢通量推断。这是重头戏，源码包默认用 Compass 做这件事。Compass 基于小鼠的全基因组代谢网络模型，结合每个细胞的表达数据，用线性规划去推断每个反应在细胞里的通量强度。它输出的不是"这条通路高不高"，而是"某个反应往哪个方向推进、推进多少"。

用表格归纳一下：

2.2 工具链选型：Compass 为主、Scanpy 为辅的原因

为什么这个源码包没有全部用 R 生态？虽然 Seurat 的用户基数大，但 Compass 本身是 Python 包，输入输出都围绕 AnnData 设计，硬要跨语言调用会多一层数据转换。源码包选择 Scanpy 做上游处理、Compass 做核心推断，一体化程度最高。

Compass 的原理值得一提。它把每个细胞当作一个独立的代谢网络优化问题来做——给定了这个细胞的基因表达水平，就约束了反应的最大速率，然后通过通量平衡分析框架去找最优通量分布。这意味着 Compass 的计算量不小，默认要跑几十分钟到几小时，具体取决于细胞数量和反应数量。

这个源码包做了个关键设计：把细胞聚类后的平均表达作为 Compass 输入，而不是每个细胞单独跑。这样做的代价是丢失单细胞内的通量异质性，但换来的是大多数普通工作站能跑完的可行性。包内同时保留了 per-cell 模式，机器好的课题组可以直接全量跑。

我一般会建议：先跑聚类平均模式得到全局代谢图，再挑感兴趣的亚群单独跑 per-cell 模式，两者结合既省时间又不牺牲关键分辨率。

提示：不要在一开始就把所有细胞丢进 Compass，默认配置下 10 万细胞的全量计算很可能会把内存吃满。先聚类、后分层跑的思路是大多数已发表工作的做法。

3. 源码包实战：从原始矩阵到代谢状态

3.1 环境搭建与依赖安装

源码包的核心依赖是 Scanpy、Compass 和用于批次整合的 Harmony，另有一些附属库用于基因名映射和绘图。我自己的习惯是用 conda 单独建一个环境，避免与已有项目互相污染：

conda create -n metasc python=3.9 -y conda activate metasc pip install scanpy harmonypy anndata pip install compass

安装时可以留意几点：Compass 在 PyPI 上的包名就叫 compass，不需要额外编译；harmonypy 负责批次整合，装的版本太新有时会跟 scanpy 依赖的 numpy 冲突，建议安装时指定numpy<1.24这一档。如果机器上有 GPU，Compass 也可以把线性规划的部分交给 GPU 加速，但我实际跑下来的感受是，多数情况 CPU 多线程更省心，GPU 版本反而容易遇到驱动兼容问题。

3.2 数据预处理与代谢基因匹配

装好环境后，第一步是把原始表达矩阵读进 AnnData，并确认基因名是小鼠的官方符号（Symbol）。很多公开数据集的基因名是 Ensembl ID，或者人类基因名，这两种情况直接跑代谢分析都会导致基因命中率暴跌，必须做映射。

import scanpy as sc import pandas as pd # 读取10x的h5文件，产出未过滤的AnnData adata = sc.read_10x_h5("mouse_data.h5") adata.var_names_make_unique() # 只保留蛋白编码基因，避免lncRNA和假基因干扰代谢打分 adata = adata[:, adata.var["gene_types"] == "protein_coding"].copy() # 归一化与对数化，Compass对输入表达值有依赖 sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata)

这段代码的注意点有两个：第一，sc.read_10x_h5读出来的gene_types列不一定存在，看数据来源，如果缺失就跳过这一行过滤；第二，target_sum=1e4是 10x 流程里最常见的测序深度归一化目标，改成1e6也不会影响 Compass 的相对关系，但会影响后续与参考数据的对比。

基因名映射我会放在这个位置，因为 Compass 自动做的符号匹配不够稳：

# 用源数据自带的映射表，把Ensembl转成Symbol ensembl_to_symbol = pd.read_csv("mart_export_mouse.txt", sep="\t") id_map = dict(zip(ensembl_to_symbol["Gene stable ID"], ensembl_to_symbol["Gene name"])) adata.var["symbol"] = adata.var.index.map(id_map) adata = adata[:, adata.var["symbol"].notna()].copy() adata.var_names = adata.var["symbol"]

参数说明：mart_export_mouse.txt是自己在 Ensembl Biomart 导出的小鼠基因映射表，只需要两列，一列是稳定 ID，一列是基因名。映射后要确认一下重复基因名的情况，比如 HIST1 家族有多条记录映射到同一个 Symbol，这时用adata = adata[~adata.var_names.duplicated()].copy()去掉重复项即可。

3.3 运行 Compass：通量推断的两种模式

这一步是源码包里最耗时也最关键的部分。前面分两种模式讲过，这里给出命令行写法。Compass 的常见参数不多，但每个都影响最终结果：

# 聚类平均模式，推荐先用这个 compass --data clustered_counts.csv \ --num-processes 8 \ --species mus_musculus \ --which-model all \ --output-dir results/ # per-cell模式，挑出感兴趣的亚群单独跑 compass --data subset_counts.csv \ --num-processes 8 \ --species mus_musculus \ --which-model all \ --output-dir results_cell/

clustered_counts.csv的格式是行为基因、列为样本，这里的样本在平均模式下是每个聚类簇，per-cell 模式下则是一个个细胞。数据需预先从 AnnData 里导出，并且基因名要和 Compass 内置的模型基因能对上。--which-model all表示跑全部模型，换成--which-model central_carbon可以只看中心碳代谢，速度快很多，适合先验证流程能跑通。

跑完后输出目录里最重要的文件是reactions.csv，它记录每个样本（聚类簇或细胞）在每个代谢反应上的通量值，行是反应，列是样本。这个文件就是后面所有可视化和统计分析的原料。

4. 避坑与排查：代谢分析里最常见的翻车现场

4.1 基因名映射后代谢基因命中率太低

现象：Compass 跑完日志里显示基因匹配数量只有几百个，而正常小鼠数据应该匹配到 3000 以上。

原因：最常见是两个，一是 Ensembl ID 没转换成 Symbol，Compass 的词表里只认标准基因名；二是数据本身来自人类样本，却指定了mus_musculus，这种匹配率基本在 20% 以下。

解决：在预处理阶段强制做一次基因名标准化，映射完成后再跑len(adata.var_names)对比一下。如果命中率还是低，把 Compass 日志里加载的基因表打出来，人工核对几个关键代谢基因（如 GAPDH、PKM、LDHA）是否在表里。

4.2 稀疏矩阵导致的内存爆炸

现象：per-cell 模式下十几个进程同时加载多个细胞的表达矩阵，内存一度顶到 60GB 以上，普通工作站直接 OOM。

原因：Compass 默认把每个细胞的表达向量转成 Dense 矩阵，10 万基因 × 数千细胞的展开量非常大。

解决：先用聚类平均模式把矩阵压到几十行；非要用 per-cell 模式就控制子集规模，一次只跑 2000 个细胞。另外--num-processes 4比8更稳，多进程带来的加速在 I/O 瓶颈前不明显。

4.3 批次效应把代谢差异完全掩盖

现象：两个实验批次（比如不同月份测的 10x 数据）合并后，代谢打分聚类出来不是按细胞类型分，而是按批次分。

原因：代谢基因本身对测序深度敏感，批次间的 library size 差异会直接传导到通路活性和通量结果上。

解决：预处理阶段加 Harmony 批次整合，把batch_key指向样本批次列，然后再做代谢打分。这一步要在 Compass 之前做，整合后重新导出表达矩阵，Compass 的输入必须用整合后的数据。

4.4 通量结果和已知生物学结论相悖

现象：文献里普遍认为肿瘤相关巨噬细胞偏向脂肪酸氧化，但跑出来的通量图显示糖酵解占主导。

原因：聚类簇划分太粗，或者误把多个亚群混进了一个簇。代谢信号很细，不同功能状态的细胞平均后会把特征抵消。

解决：回到聚类参数，把n_neighbors从 15 调到 30、resolution从 0.5 调到 1.0，看是否可以分出更细的亚群；再用 AUCell 对单个细胞做脂肪酸氧化基因集打分，与 Compass 的簇级结果互相印证。

4.5 物种参数写错导致全盘白跑

现象：所有流程走完，发现反应列表里出现大量人类才有的反应（比如 CYP 酶系），和小鼠生物学不一致。

原因：--species参数写成了homo_sapiens，或者用了人类模型谱。

解决：这一步在源码包的解压目录里有一个config.yaml可以写默认值，我一般会在里面固定mus_musculus，避免每次命令行都手敲出错。基因名处理时也可以顺带加一个物种名列，在映射阶段就做一次交叉校验。

注意：Compass 的--which-model all会同时跑多个代谢子模型，输出文件数量很多，磁盘空间不够也是容易忽略的坑，提前留出至少 20GB 比较稳妥。

5. 进阶技巧：代谢状态可视化与下游验证

5.1 把通量矩阵映射回 UMAP

Compass 输出的reactions.csv行是反应、列是样本（聚类簇），要把它和单细胞 UMAP 对齐，需要先把通量数据降维，再把每个簇的通量值回填到该簇所有细胞上。

import scanpy as sc import pandas as pd import numpy as np reactions = pd.read_csv("results/reactions.csv", index_col=0) # 用PCA把高维反应空间压到20维，方便后续与UMAP叠加 pca = sc.tl.pca(reactions.T, n_comps=20) # 将每个簇的通量PCA坐标赋值给该簇内的所有细胞 adata.obs["metab_pca_1"] = adata.obs["leiden"].map( dict(zip(pca.index, pca[:, 0])) ) adata.obs["metab_pca_2"] = adata.obs["leiden"].map( dict(zip(pca.index, pca[:, 1])) ) sc.pl.umap(adata, color=["metab_pca_1", "metab_pca_2"], cmap="RdBu_r")

逻辑说明：pca计算时列是反应描述、行是聚类簇，所以需要转置成簇 × 反应矩阵再丢给sc.tl.pca。回填时用leiden聚类标签做键，把每个簇的 PCA 坐标广播回所有细胞，最终形成一张代谢视角的 UMAP。RdBu_r色板在高值和低值之间对比明显，适合代谢通量这种连续变量。

参数说明：n_comps=20不是硬编码，反应数量在几千个时取 20 维足够解释大部分方差；如果只关注单一代谢模块（比如氨基酸代谢），可以直接用那部分反应做子集再降维，比全量 PCA 更有生物学针对性。

5.2 代谢特征与差异基因的交叉验证

通量结果说到底是一个模拟结果，需要用实验或已知文献去校验。这个源码包里附带了一个相关性分析模块，功能是把代谢通量和差异表达基因做关联，用来确认代谢变化在表达层面有支撑。

# 计算每个反应与目标基因表达量的Spearman相关 target_gene = "Cpt1a" # 脂肪酸氧化关键酶 reaction_scores = reactions.T # 簇 × 反应 corr = reaction_scores.corrwith( adata.obs.groupby("leiden")[target_gene].mean() ) print(corr.sort_values(ascending=False).head(10))

这一步的思路是：如果 Compass 推断出脂肪酸氧化通量升高，而Cpt1a的表达在同一批簇里也同步上调，两者相互印证，结论可信度大幅提升。反过来如果通量高但关键酶表达没有变化，就要警惕是不是网络模型对某些反应的约束过松导致误报，回到参数检查。

5.3 运行时间与参数的折中方案

单次 full 流程（2 万细胞，平均模式）在我日常用的 8 核 32GB 机器上大约要跑 3 到 4 小时。第一次调试肯定会跑崩，我强烈建议先用--which-model central_carbon跑一个几百细胞的子集，把链路全走通再放全量。全量跑的时候--num-processes设置为 CPU 核数减一，预留出一个核给系统和其他任务，--lambda 0.01是这套源码包里默认的平衡参数，控制模型对表达数据的信任程度，不建议新手动它。

我自己的习惯是把每次 Compass 运行的参数记录在config.yaml里，基因命中率、运行时长、反应数都写进去。从那以后每次代谢分析我都强制走一遍这个流程：先子集试跑、再全量、最后用Cpt1a这类关键酶做交叉验证。这步能筛掉九成以上的流程错误，希望帮到你。

本文还有配套的精品资源，点击获取