小鼠单细胞代谢分析从原始矩阵到代谢活性：AUCell与ssGSEA完整源码路线

简介：面向单细胞转录组研究人员与生物信息学入门者，这份源码包提供基于scMetabolism包的小鼠单细胞代谢激活分数分析实现，重点解决小鼠基因名向人类基因名转换、以及适配Seurat v4/v5版本等实际操作问题，能帮助用户在单细胞层面量化代谢通路活性。包体非常轻量，压缩后仅9KB，共包含6个文件，其中以R脚本和Markdown说明文档为主，另有HTML预览页、配置文件和依赖安装脚本，便于快速搭建运行环境、对照学习与直接复用。目前已有180人学习下载，适合作为代谢分析入门参考。源码覆盖了从基因表达数据读取、基因名转换到代谢激活分数计算与Seurat对象集成的完整流程，并配有代码注释和说明文档；读者可据此理解scMetabolism包的数据处理逻辑，并迁移至自身的小鼠单细胞数据中，也可通过参考链接进一步探索代谢通路背景知识。整体来看，这是一份精简实用的代谢分析工具包，能有效降低单细胞代谢研究的代码门槛，尤其适合想将代谢分析纳入现有Seurat工作流的科研人员。

1. 小鼠单细胞代谢分析怎么做：从原始矩阵到代谢活性的完整源码路线

小鼠单细胞代谢分析在做一件听起来很直接、实际很别扭的事：从单细胞转录组里推断每个细胞的代谢活性。直接看代谢酶的表达量会翻车，酶有翻译后修饰、有底物浓度影响，表达高不代表通路就开着。源码路线一般落在 AUCell 或 ssGSEA 打分上，把代谢基因集的表达排名转成一个分数，再贴回聚类、拟时序和病理分组里解读。适合做肿瘤微环境、发育、肝病或代谢疾病模型的科研团队。这篇笔记按数据准备、打分、可视化、验证的顺序拆开讲，附能直接跑的 R 代码，也把小鼠数据和人类数据的差异单独拎出来说。

2. 数据准备与打分原理：为什么代谢分析不能只盯表达量

2.1 表达量不是活性：选 AUCell 还是 ssGSEA

单细胞数据落到手，最常见的第一反应是看基因表达：代谢酶高就是代谢旺。这个推理在代谢分析里站不住。代谢通路的通量由酶活性、辅因子浓度、底物供应共同决定，转录水平只是其中一环。比如糖酵解里 PFK 的转录变化很小，但活性被 ATP 和柠檬酸强抑制；脂肪酸氧化相关基因在静息 T 细胞里也表达，可细胞根本不在氧化脂肪。所以代谢分析实际做的是“活性推断”，不是“表达汇总”。

目前从业人员用得最多的两类方法是 AUCell 和 ssGSEA。ssGSEA 的思路是把基因集在全体细胞的表达排名上做富集，给出一个连续分；AUCell 的思路更偏单细胞，它先给每个细胞内部所有基因按表达量排座次，再看目标基因集在排名中的位置是否显著靠前。两者的差别决定了适用场景：ssGSEA 对通路整体状态敏感，AUCell 对细胞个体更稳。我一般建议以 AUCell 为主、ssGSEA 作交叉验证，理由在 2.4 展开。

2.2 代谢基因集准备：小鼠符号版本别拿人的直接用

基因集是整个代谢分析的地基。常见做法是去 MSigDB 下载小鼠符号版（文件名类似 mh.all.v2023.Mm.symbols.gmt），注意不要误下人源注释文件。如果手里只有人类基因集，就得先做同源转换，这一步常常是后面分数全零的元凶。

我用 GSEABase 读 GMT 的代码很固定，直接放在这里。基因集文件放在本地 data 目录，不依赖网络。

library(GSEABase) gsc <- getGmt("data/mh.all.v2023.Mm.symbols.gmt") metab_pathways <- gsc[grep("KEGG|REACTOME", names(gsc))] metab_genes <- unique(unlist(geneIds(metab_pathways))) length(metab_genes)

逻辑说明：getGmt把 GMT 读成 GeneSet 集合，grep("KEGG|REACTOME", ...)只挑代谢相关通路，避免把免疫、增殖的基因混进来。最后用geneIds()提取基因符号，unlist后去重。参数调整上，如果只想看糖酵解和氧化磷酸化，可以进一步按通路名过滤，比如加grep("GLYCOLYSIS|OXIDATIVE_PHOSPHORYLATION", ...)。这里最容易踩的坑是geneIds()返回的是列表，忘了unlist会导致后面矩阵匹配全空。

如果你从 KEGG 官网手工整理基因列表，强烈建议整理成三列 CSV：通路名、基因符号、物种。至少留一个文件备份，不然日后审稿人问“基因集哪来的”，你拿不出一份可复现的清单，会很被动。

2.3 Seurat 标准预处理：小鼠线粒体基因前缀是 mt- 不是 MT-

拿到 10X 的 filtered_feature_bc_matrix 后，先用 Seurat 建对象。这里给的是最小可用代码，后面所有代谢分析都基于这个对象跑，不另起炉灶。

library(Seurat) counts <- Read10X("data/filtered_feature_bc_matrix") obj <- CreateSeuratObject(counts = counts, project = "mouse_metab", min.cells = 3, min.features = 200) obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^mt-") obj <- subset(obj, subset = nFeature_RNA > 500 & nFeature_RNA < 6000 & percent.mt < 20) obj <- NormalizeData(obj) %>% FindVariableFeatures(nfeatures = 3000) %>% ScaleData() %>% RunPCA(npcs = 30) %>% RunUMAP(dims = 1:20)

参数说明：min.cells = 3要求基因至少在 3 个细胞中出现，过滤掉文库污染带来的假基因；min.features = 200保证每个细胞有足够多的检测基因，空液滴通常低于这个值。pattern = "^mt-"针对小鼠线粒体基因，前缀是小写 mt-，人类才是 MT-，写错就彻底算不出线粒体比例。nFeature_RNA 上限 6000 和 percent.mt 上限 20 是常规过滤阈值，如果组织是肝或心肌，线粒体比例本身偏高，上限可以放宽到 30，按实际分布图调整。

一个容易被忽略的点：NormalizeData默认是 LogNormalize，这一步之后的 data 矩阵才是 AUCell 和 ssGSEA 的输入。不要在 scale.data 上打分，scale 后的残差会把代谢基因的低幅波动放大成伪信号。

2.4 为什么不先跑差异基因再打分

很多新手习惯先 FindMarkers 挑差异基因，再拿差异基因去做通路富集。这在代谢分析里是绕远路。代谢相关基因大多是持家基因，表达变化幅度小，差异检验往往不显著，直接筛选就把信息丢光了。AUCell 和 ssGSEA 的共性优点是它们用“基因集整体在排名中的位置”打分，不需要单个基因显著差异。所以基因集的完整度比差异倍数重要得多。

换句话说，代谢分析的正确姿势是：先准备好完整的代谢基因集，再在归一化矩阵上打分，最后把分数拿去和聚类、拟时序、分组信息做关联。顺序反了，结论就会偏。这个逻辑后面每一章都在围绕它展开。

3. 用 AUCell 给每个细胞打代谢分：源码与参数说明

3.1 AUCell_buildRankings：先给每个细胞排座次

AUCell 的第一步不是直接打分，而是构建 rankings。它对每个细胞的所有基因按表达量从高到低排序，得到一个基因排名矩阵。一个细胞表达量再稀疏，排名矩阵也是完整的，这就是 AUCell 在单细胞 dropout 严重时依然稳定的原因。

library(AUCell) cells_rankings <- AUCell_buildRankings(obj@assays$RNA@data, plotStats = TRUE)

逻辑说明：输入是归一化后的 data 矩阵，不是 count，也不是 scale.data。plotStats = TRUE会输出每个细胞检测基因数的分布图，方便你确认细胞质量。这一步计算量不小，三万细胞大概要几分钟，如果内存有限，可以先随机抽两万细胞做测试，确认通路打分逻辑没问题后再全量跑。rankings 对象是后面所有计算的中间产物，建议用saveRDS存下来，避免重复构建。

3.2 AUCell_calcAUC：aucMaxRank 是该调的第一个参数

按基因集算 AUC 这一步，核心参数是aucMaxRank，默认值是排名前 5% 的基因数量。默认值在大部分通路场景够用，但代谢基因集往往很大，包含几百个基因，如果aucMaxRank设得太小，基因集里有一半基因排在 5% 之外，截断后打分区分度会变差。

metab_genes <- unique(unlist(geneIds(metab_pathways))) cells_AUC <- AUCell_calcAUC(list(metabolism = metab_genes), cells_rankings, aucMaxRank = nrow(cells_rankings) * 0.05)

参数说明：list(metabolism = metab_genes)这里给基因集起名字，后面取分数时就用这个名字。aucMaxRank = nrow(cells_rankings) * 0.05就是把默认阈值显式写出来，方便后续调整。我一般会先跑默认 5%，再试 10% 和 20%，看三个阈值下分数在 UMAP 上的分布是否稳定。如果阈值变化后结论完全反转，说明基因集本身信号弱，不是参数问题。

AUCell_calcAUC的返回值是 AUCell 对象，里面包含每个细胞在每个基因集上的 AUC 分数。分数越高，说明该细胞中代谢基因集的整体表达排名越靠前，代谢活性推断越强。

3.3 ssGSEA 交叉验证：别只信单一方法

单靠 AUCell 不够。AUCell 对基因集大小敏感，ssGSEA 对基因集组成的敏感方向不同，两个方法结论一致时，可信度才够。用 GSVA 包跑 ssGSEA 是最常见的做法。

library(GSVA) expr_mat <- as.matrix(obj@assays$RNA@data) gset_list <- list(metabolism = metab_genes) ssgsea_res <- gsva(expr_mat, gset_list, method = "ssgsea")

逻辑说明：gsva返回一个基因集乘以细胞的矩阵，每一列就是一个细胞的 ssGSEA 分数。method = "ssgsea"用的是单样本 GSEA 思路，不依赖分组标签，适合单细胞这种无先验分组的数据。这里必须用归一化矩阵转成普通矩阵，GSVA 不接受稀疏矩阵。一个隐藏问题：gsva默认会做基因注释匹配，如果你的基因名是“Gapdh”这种小鼠写法，而基因集里混了人类写法“GAPDH”，匹配率会骤降，分数全往零缩。匹配率低于 70% 就要回头查基因集。

3.4 把分数写回 Seurat：AddMetaData 与 FeaturePlot

打分是中间状态，最终要让代谢分数和细胞身份一起看。把 AUCell 分数写回 Seurat 元数据，再用 FeaturePlot 做第一眼判断。

obj$AUC_metab <- as.numeric(getAUC(cells_AUC)["metabolism", ]) obj$ssGSEA_metab <- as.numeric(ssgsea_res["metabolism", ]) cor(obj$AUC_metab, obj$ssGSEA_metab, method = "spearman") FeaturePlot(obj, features = "AUC_metab", order = TRUE)

逻辑说明：getAUC(cells_AUC)返回基因集 × 细胞矩阵，用["metabolism", ]取目标基因集那一行，转成数值后直接存入obj$AUC_metab。cor计算 AUCell 和 ssGSEA 的 Spearman 相关，相关系数低于 0.3 就要警惕：要么基因集有问题，要么数据质量问题，两个方法意见完全不同时别急着下结论。FeaturePlot的order = TRUE让高分细胞最后画、显示在最上层，避免低分细胞把画面盖住。

到这里，每个细胞已经有一个代谢分数。下一步不是直接信这个分数，而是要把它拆到通路级别去解读，否则你会看到一张“全图一个颜色”的 UMAP。

4. 代谢异质性的解读：聚类、轨迹与通路级三层视角

4.1 通路级打分比全局打分更有信息量

把几百个代谢基因揉成一个大集合，打出来的分数往往在 UMAP 上是连续渐变，看不出结构。原因很简单：糖酵解、氧化磷酸化、脂肪酸氧化、氨基酸代谢是四条不同逻辑的通路，有的细胞靠糖酵解，有的靠脂肪酸氧化，合在一起平均，信号互相抵消。

常见做法是按 KEGG 二级通路拆成多个基因子集，分别打分。糖酵解和 TCA 循环、脂肪酸氧化、氨基酸代谢各算一个分数，再去看哪个通路驱动了聚类差异。

glycolysis_genes <- unlist(geneIds(gsc[["KEGG_GLYCOLYSIS_GLUCONEOGENESIS"]])) oxphos_genes <- unlist(geneIds(gsc[["KEGG_OXIDATIVE_PHOSPHORYLATION"]])) fa_ox_genes <- unlist(geneIds(gsc[["KEGG_FATTY_ACID_OXIDATION"]])) pathway_list <- list(glycolysis = glycolysis_genes, oxphos = oxphos_genes, fa_ox = fa_ox_genes) pathway_AUC <- AUCell_calcAUC(pathway_list, cells_rankings)

逻辑说明：gsc[["KEGG_GLYCOLYSIS_GLUCONEOGENESIS"]]直接从 GeneSet 集合里按名字取子集，geneIds再提取基因。这里注意 KEGG 里脂肪酸氧化通路名在不同版本里叫法不一致，有的版本叫KEGG_FATTY_ACID_METABOLISM，需要先names(gsc)确认。算出三条通路的 AUC 后，逐一写回 Seurat 元数据。

我在实际项目里发现，糖酵解和氧化磷酸化两条通路往往呈现互斥分布，这在肿瘤微环境里尤其明显，CD8 T 细胞耗竭亚群通常糖酵解分数偏低而线粒体氧化分数偏高。单看全局代谢分数永远看不出这个规律。

4.2 UMAP 与分组小提琴图：代谢分数的分布形态比颜色更重要

FeaturePlot 适合找感觉，不适合下结论。代谢分数在细胞群体内往往不是正态分布，而是右偏的，少数细胞分数极高，多数细胞贴底。这种分布形态用小提琴图比散点图清晰得多。

p1 <- VlnPlot(obj, features = "AUC_glycolysis", group.by = "cell_type") + geom_boxplot(width = 0.1, outlier.shape = NA) p2 <- VlnPlot(obj, features = "AUC_oxphos", group.by = "cell_type") + geom_boxplot(width = 0.1, outlier.shape = NA)

逻辑说明：group.by换成你自己的细胞类型注释列，通常是harmony_clusters或cell_type。加一层geom_boxplot是因为 VlnPlot 默认的密度曲线容易让人忽略中位数差异，箱线图把中位数和四分位间距暴露出来。比较两组的代谢分数时，别只看均值，用 Wilcoxon 检验补一个 p 值，但要注意细胞数多时 p 值必然显著，效应量（两组中位数差）比 p 值更有意义。

另外提示一句：FeaturePlot 的颜色刻度默认是每个面板独立缩放的，两张子图拼一起时，务必确认两个分数面板的颜色范围一致，否则视觉上的差异是假的。

4.3 拟时序上的代谢活性：monocle3 里看动态变化

代谢状态不是静止的。从初始状态到效应状态、从幼稚到耗竭，代谢重编程是连续的。拟时序分析能把这种连续变化投射到一维轨迹上，再和代谢分数做关联。

library(monocle3) cds <- as.cell_data_set(obj) cds <- cluster_cells(cds) cds <- learn_graph(cds, use_partition = TRUE) cds <- order_cells(cds) obj$pseudotime <- pseudotime(cds) plot(obj$pseudotime, obj$AUC_oxphos, pch = 16, cex = 0.3, col = "grey40") lines(lowess(obj$pseudotime, obj$AUC_oxphos), col = "red", lwd = 2)

逻辑说明：as.cell_data_set把 Seurat 对象转成 monocle3 格式，learn_graph构建轨迹图，order_cells需要在交互窗口选根细胞。根的选取很重要，我一般选对照组中 marker 基因表达典型的细胞群作为起点。lowess是局部加权回归，专门用来无视噪声看趋势；如果曲线呈现先升后降的倒 U 形，说明代谢状态是过渡态而非终态，这时候再去比较终态两端反而会漏掉中间的高代谢细胞。

细胞量大的时候，三万细胞全画成散点会糊成一团。常规做法是先按分数分箱，再对每箱计算中位数和置信区间，用误差条替代散点。

4.4 代谢亚型划分：什么时候按分数二分才合理

把连续分数变成离散分组，本身是一种降维操作，有信息损失。只有当分数分布明显双峰，或者分数和某个功能表型强相关时，才建议用阈值分组。

thresh_high <- quantile(obj$AUC_oxphos, 0.75) thresh_low <- quantile(obj$AUC_oxphos, 0.25) obj$oxphos_group <- ifelse(obj$AUC_oxphos >= thresh_high, "high", ifelse(obj$AUC_oxphos <= thresh_low, "low", NA))

参数说明：用上下四分位数取 Top 25% 和 Bottom 25%，中间 50% 标成 NA，不参与后续差异分析。这样做的逻辑是中间态本身噪声大，区分度不足，强行二分会把真实信号稀释掉。分组后可以做代谢基因差异表达、转录因子富集、以及和其他功能标记的相关分析。注意不要用中位数二分，中位数切出来的两组差异通常经不起置换检验。

5. 小鼠单细胞代谢分析的五个常见坑与排查记录

5.1 基因集匹配率极低，分数全部趋近于零

现象：AUCell 打分跑完，getAUC出来几乎所有细胞分数都集中在 0 到 0.05，UMAP 上没有任何梯度。

原因：最常见的是小鼠数据配了人类基因集。小鼠基因符号首字母大写后面小写，如“Gapdh”，人类是全大写“GAPDH”。还有一种情况是 GMT 文件本身只带了 Entrez ID，代码里没有转换，geneIds提取到的是一串数字而不是符号。

解决：验证基因与表达矩阵 rownames 的匹配率，代码写mean(metab_genes %in% rownames(obj))。低于 0.7 直接换小鼠符号版基因集。如果是 Entrez ID，用bitr或同源转换表映射到 symbol 后再跑。

5.2 aucMaxRank 调小后通路分数全零

现象：同一份数据，aucMaxRank设为默认值时有分数，改成特别低的值后一个大群分数变成零。

原因：代谢基因集合通常有几百个基因，aucMaxRank 过窄意味着只允许排名前 1% 基因参与计算，代谢基因大部分排在阈值之外，AUC 自然趋近于零。

解决：将aucMaxRank设为排名前 10% 或前 20%，跑一遍看分布；再用两个阈值的结果做 Spearman 相关。相关大于 0.8 说明阈值不敏感，结论可靠；小于 0.5 说明基因集信号弱，要回头审查基因集质量，而不是继续调参数。

5.3 三个不同通路分数高度相关，区分不出功能差异

现象：糖酵解、氧化磷酸化、脂肪酸氧化分数两两相关系数超过 0.9，UMAP 上三张图几乎一模一样。

原因：基因集之间有重叠，或者打分方法被总表达量主导。高表达基因多的细胞，任何基因集排名都靠前，本质上是细胞大小或测序深度效应。

解决：检查通路基因重叠长度；如果重叠超过 30%，先做减法去重。另外把代谢分数对 nFeature_RNA 和 nCount_RNA 做回归，取残差再作图。常见做法是确认残差后的分数仍然与分组相关，说明信号不是文库大小带来的。

5.4 代谢分数与批次完全绑定

现象：两个不同批次的小鼠样本合并分析，UMAP 上按批次分堆，代谢分数也随批次明显不同。

原因：未做批次整合。代谢基因表达受实验批次影响大，尤其是组织消化时间不同，线粒体基因比例与代谢基因表达联动变化。

解决：先对 obj 做整合。Seurat 的整合流程或者 Harmony 都可以，这里给出常用代码：obj <- RunHarmony(obj, group.by.vars = "batch")，之后 PCA 和 UMAP 一律基于 harmony 降维结果。整合后再重新提取代谢分数，对比整合前后分数分布，如果差异消失，说明先前的代谢差异是批次假象。

5.5 ssGSEA 在 scale.data 上跑，分数偏离真实表达

现象：AUCell 结果和 ssGSEA 结果完全相反，同一细胞一个分数高一个分数低。

原因：ssGSEA 输入用了 scale 后的数据，负值残差和零膨胀打乱了基因排名，ssGSEA 对表达值的符号非常敏感。

解决：ssGSEA 输入始终用GetAssayData(obj, slot = "data")，确保是归一化非缩放数据。AUCell 输出与 ssGSEA 相关低于 0.3 时，先做这一项检查，再决定重新跑哪一步。

6. 让代谢结论站得住：一个置换检验的验证习惯

代谢分析分值本身没有绝对意义，0.3 和 0.6 相差多少，不同基因集定义下不可比。所以真正重要的不是分数大小，而是分数差异能否在随机背景下复现。我习惯做置换检验：随机抽取与目标基因集同等大小的基因集，重复计算 AUC 或 ssGSEA 分数多次，得到零分布，再看真实的代谢分数差异落在零分布的什么位置。

set.seed(42) n_perm <- 100 null_auc <- matrix(nrow = n_perm, ncol = ncol(cells_AUC)) bg_genes <- rownames(obj@assays$RNA@data) for (i in seq_len(n_perm)) { set.seed(100 + i) random_genes <- sample(bg_genes, length(metab_genes)) rand_auc <- AUCell_calcAUC(list(random = random_genes), cells_rankings, aucMaxRank = ncol(cells_rankings) * 0.05) null_auc[i, ] <- getAUC(rand_auc)["random", ] }

逻辑说明：循环里每次抽一组与代谢基因集等长的随机基因，跑同样的 AUCell 流程，得到 100 行背景分数矩阵，每一行都是一个随机基因集的细胞分数向量。然后计算真实代谢分数在每列的位置，比如统计真实分数高于 95% 随机分布的比例，这个比例高于 0.9 就说明信号不是尺寸效应，低于 0.6 就说明可以收工换基因集了。set.seed是保证这份结果在别人手里能复现，这在审稿阶段特别重要。

如果对照组设计明确，也可以做样本标签置换：把高分和低分细胞的标签随机打乱，重复 1000 次，观察组间中位数差的分布，真实差异落在 95% 区间之外才算站得住。这个办法对样本数少的数据尤其重要——三个样本组里有一个异常样本就足以制造假阳性。

我自己吃过一次亏，肝病模型里高代谢亚群与病理分期显著相关，结果发现是测序深度差异导致的假象，置换检验第一遍就跑出 p=0.6，白高兴两周。后来把置换检验写进标准化流程，当成和质控图同级的固定环节，不再跳过。做单细胞代谢分析，分数只是起点，能经受住随机背景拷问的差异，才是值得写进结论的东西。希望这个习惯对你也有用。

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