简介：这份源码资源面向从事单细胞转录组与代谢研究的科研人员及生物信息学初学者，围绕scMetabolism包解决小鼠单细胞代谢激活分数分析问题，重点处理小鼠基因名向人类基因名的转换，并适配Seurat v4与v5版本，帮助读者在R环境中完成从表达数据到代谢通路评分的完整流程。资源包共6个文件，以R脚本为主，包含代谢分析主程序与依赖安装脚本，另附HTML说明页、Markdown文档及项目配置文件，压缩包约9KB，结构轻量便于快速上手。目前已有176人学习下载。读者可获取可直接运行的代码示例、基因名转换与Seurat对接思路，以及参考链接指引，适合需要将小鼠单细胞数据纳入代谢维度分析、探索细胞状态与疾病机制的研究场景。
1. 小鼠单细胞代谢分析源码：从矩阵到代谢通路的落地路径
单细胞转录组做聚类、做注释、做拟时序，这些流程已经相当成熟，但一提到代谢分析，很多人就卡住了。手里明明有表达矩阵，却不知道怎么把它变成每个细胞的代谢活性评分，更不知道那些评分背后对应的是糖酵解、氧化磷酸化还是脂肪酸氧化。小鼠单细胞代谢分析源码要解决的正是这个问题：把单细胞表达数据映射到代谢通路，算出每个细胞的代谢状态，再据此做分群、做差异、做可视化。适合已经跑过 Seurat 或 Scanpy 基础流程、想往代谢方向延伸的从业者，也适合做肿瘤微环境、免疫代谢、发育代谢重编程的团队直接复用。核心思路不复杂——用代谢基因集给细胞打分，但打分方式、基因集来源、归一化策略，每一步都有讲究。
2. 代谢基因集与打分算法的选型逻辑
2.1 为什么不用普通通路富集直接套单细胞
批量 RNA 的通路富集工具，比如 GSEA 或 GSVA，默认样本是「批量」的，输入是一个基因表达矩阵加一组表型标签。单细胞数据动辄几千到几万个细胞，如果直接把每个细胞当成一个样本丢进去，会碰到两个硬伤：第一，单细胞表达矩阵极度稀疏，大量基因在单个细胞里是零，富集算法对零值敏感，结果会被 dropout 事件主导；第二，代谢通路的基因数量通常不大，十几个到几十个基因，在单细胞层面做富集，统计功效很低，容易出现假阴性。
常见做法是换一套思路：不做富集，做打分。给每个细胞算一个代谢通路活性分数，分数高低反映该通路在该细胞中的相对活跃程度。打分方法有几种，最常用的是 Seurat 的AddModuleScore，它把目标基因集的平均表达作为原始分数，再减去随机背景基因集的平均表达，得到一个校正后的分数。这个方法的优点是快、稳、对稀疏数据有一定容忍度，缺点是它假设基因之间独立，不考虑通路内部的调控关系。
另一条路是AUCell，它不直接算平均表达，而是对每个细胞的基因表达排序，看目标基因集是否富集在排序顶部。AUCell 对 dropout 更鲁棒，但计算量更大，几万个细胞跑起来需要并行。还有ssGSEA的单细胞版本，原理和 GSVA 类似，但实现上做了单细胞适配。选哪个，取决于你的数据规模和下游分析目标。如果只是做初步探索，AddModuleScore够用；如果要发文章、做精细比较，建议用AUCell或ssGSEA做交叉验证。
2.2 小鼠代谢基因集的获取与整理
人和小鼠的代谢基因集有现成的资源，但直接拿来用会踩坑。KEGG 通路里代谢相关条目很多，但 KEGG 的基因 ID 是 Entrez，而单细胞矩阵通常是 Symbol，需要做 ID 转换。Reactome 的代谢通路更细，但条目太多，直接全用会导致多重检验负担过重。MSigDB 的 Hallmark 基因集里有一组代谢相关通路，比如HALLMARK_GLYCOLYSIS、HALLMARK_OXIDATIVE_PHOSPHORYLATION、HALLMARK_FATTY_ACID_METABOLISM，数量适中，适合单细胞打分。
我一般会从 MSigDB 下载小鼠对应的基因集，或者用msigdbr包直接提取。注意，MSigDB 的小鼠基因集是通过同源映射从人转换过来的，部分基因可能丢失或一对多映射。如果做的是小鼠特有代谢过程，比如某些肝脏特有的代谢通路，建议手动补充基因列表，来源可以是 KEGG 小鼠通路或文献。
# 加载必要的包 library(Seurat) library(msigdbr) library(dplyr) # 提取小鼠 Hallmark 代谢相关基因集 m_df <- msigdbr(species = "Mus musculus", category = "H") # 筛选代谢相关通路 metabolic_pathways <- c( "HALLMARK_GLYCOLYSIS", "HALLMARK_OXIDATIVE_PHOSPHORYLATION", "HALLMARK_FATTY_ACID_METABOLISM", "HALLMARK_P53_PATHWAY" # 与代谢应激相关，可选 ) metabolic_sets <- m_df %>% filter(gs_name %in% metabolic_pathways) %>% split(x = .$gene_symbol, f = .$gs_name) # 查看每个通路的基因数 sapply(metabolic_sets, length)这段代码先加载msigdbr，指定物种为小鼠、类别为 Hallmark，然后筛选出四个代谢相关通路。split函数把数据框按通路名拆成列表，每个元素是一个通路的基因 Symbol 向量。最后一行查看每个通路的基因数量，一般糖酵解和氧化磷酸化各有 200 个左右基因，脂肪酸代谢约 150 个。如果某个通路基因数少于 30，说明映射过程中丢失太多，需要手动补充。
参数说明：species必须写"Mus musculus"，写"mouse"会报错；category = "H"表示 Hallmark 集合，如果要更细的代谢通路，可以换成category = "C2"配合subcategory = "CP:KEGG"，但基因集数量会大幅增加，后续要做多重检验校正。
2.3 用 AddModuleScore 给每个细胞打代谢分
拿到基因集后，下一步是给 Seurat 对象里的每个细胞打分。AddModuleScore是 Seurat 内置函数，用法简单，但有几个参数必须调对。
# 假设 seurat_obj 已经完成标准化，且细胞类型注释已完成 # 给每个代谢通路打分 for (pathway in names(metabolic_sets)) { seurat_obj <- AddModuleScore( object = seurat_obj, features = list(metabolic_sets[[pathway]]), name = paste0(pathway, "_Score"), nbin = 24, # 背景基因分箱数 ctrl = 100, # 每个细胞选取的背景基因数 seed = 42 # 随机种子，保证可重复 ) } # 查看打分结果列名 grep("_Score", colnames(seurat_obj@meta.data), value = TRUE)AddModuleScore的核心逻辑是：对目标基因集，计算每个细胞的平均表达值；然后从表达量相近的基因中随机抽取ctrl个背景基因，计算背景平均表达；两者相减得到校正分数。nbin = 24表示把基因按平均表达量分成 24 个箱，从同一箱里抽背景基因，这样背景基因的表达分布和目标基因更接近，校正更合理。ctrl = 100是每个细胞抽 100 个背景基因，这个值不能太小，否则背景噪声大；也不能太大，否则计算慢。seed固定后结果可重复，这在做差异分析时很重要。
打分完成后，seurat_obj@meta.data里会多出几列，列名是HALLMARK_GLYCOLYSIS_Score1这样的格式。注意，Seurat 会自动在名字后面加数字，如果多次运行同一个名字，数字会递增。建议在循环里用paste0拼一个唯一名字，避免混淆。
2.4 打分结果的归一化和可视化
打分出来不能直接用，因为不同通路的分数范围不同，有的通路分数在 -0.5 到 0.5 之间，有的在 -1 到 1 之间。做跨通路比较时，需要先归一化。我一般用 z-score 归一化，把每个通路的分数转成标准正态分布，这样不同通路之间可以横向比较。
# 提取打分列 score_cols <- grep("_Score", colnames(seurat_obj@meta.data), value = TRUE) # z-score 归一化 seurat_obj@meta.data[score_cols] <- scale(seurat_obj@meta.data[score_cols]) # 可视化：用 FeaturePlot 看糖酵解分数在 UMAP 上的分布 FeaturePlot( seurat_obj, features = "HALLMARK_GLYCOLYSIS_Score1", cols = c("lightgrey", "red"), min.cutoff = -1, max.cutoff = 1 ) + ggtitle("Glycolysis Score")scale函数默认对每列做 z-score，结果是一个矩阵，赋值回meta.data时要注意列名对齐。FeaturePlot的min.cutoff和max.cutoff用来截断极端值，避免个别细胞分数过高导致颜色映射失真。一般取 -1 到 1 或 -2 到 2，根据实际分布调整。
如果想看不同细胞类型的代谢分数差异，可以用VlnPlot或DotPlot。DotPlot更适合展示多个通路在多个细胞类型中的平均分数，点的大小表示表达比例，颜色表示平均分数。
# 按细胞类型展示代谢分数 Idents(seurat_obj) <- "cell_type" # 假设已有细胞类型注释 DotPlot( seurat_obj, features = score_cols, group.by = "cell_type", cols = c("blue", "white", "red") ) + RotatedAxis()DotPlot的features传入打分列名，group.by指定细胞类型列。颜色映射用蓝-白-红，蓝色表示低分，红色表示高分。RotatedAxis把 x 轴标签旋转 45 度，避免重叠。
3. 从代谢分数到生物学结论：差异分析与通路关联
3.1 代谢分数的差异比较与统计检验
拿到每个细胞的代谢分数后，下一步是比较不同组别或不同细胞类型之间的差异。比如比较肿瘤细胞和正常细胞的糖酵解分数，或者比较不同亚群的氧化磷酸化水平。这里要注意，单细胞数据的统计检验不能用普通的 t 检验，因为细胞之间不独立，同一患者的细胞有批次效应。
常见做法是先用FindMarkers做差异表达，但FindMarkers默认是对基因表达做检验，不是对代谢分数。要对代谢分数做检验，可以手动提取分数列，用wilcox.test或limma做。如果样本有多个生物学重复，建议用混合效应模型或 pseudobulk 方法，把同一患者的细胞聚合成一个样本，再做组间比较。
# 提取代谢分数和分组信息 score_data <- seurat_obj@meta.data[, c(score_cols, "group", "sample_id")] # 按样本聚合，取平均分 pseudobulk <- score_data %>% group_by(sample_id, group) %>% summarise(across(all_of(score_cols), mean)) # 用 limma 做差异分析 library(limma) design <- model.matrix(~ group, data = pseudobulk) fit <- lmFit(t(pseudobulk[, score_cols]), design) fit <- eBayes(fit) topTable(fit, coef = 2, adjust.method = "BH")这段代码先把细胞水平的分数按样本聚合，每个样本每个通路得到一个平均分。然后用limma做线性模型拟合，design矩阵里group是分组变量。topTable输出差异分析结果，coef = 2表示比较组 vs 对照组的系数，adjust.method = "BH"做 Benjamini-Hochberg 多重检验校正。如果pseudobulk里样本数少于 3，limma的结果不稳定，建议改用非参数检验或增加样本量。
3.2 代谢通路之间的相关性分析
代谢通路不是孤立的，糖酵解和氧化磷酸化之间往往有代偿关系，脂肪酸氧化和糖酵解也可能此消彼长。分析通路之间的相关性，可以发现代谢重编程的模式。我一般会计算通路分数的 Spearman 相关系数，然后做聚类热图。
# 计算通路之间的 Spearman 相关性 cor_mat <- cor(seurat_obj@meta.data[, score_cols], method = "spearman") # 可视化相关性热图 library(pheatmap) pheatmap( cor_mat, cluster_rows = TRUE, cluster_cols = TRUE, display_numbers = TRUE, number_format = "%.2f", color = colorRampPalette(c("blue", "white", "red"))(100), main = "Metabolic Pathway Correlation" )cor函数计算列之间的相关系数，method = "spearman"对非正态分布更稳健。pheatmap做聚类热图，display_numbers在格子里显示相关系数，number_format控制小数位数。如果某些通路之间相关系数绝对值大于 0.6，说明它们在该数据集中高度协同或拮抗，值得进一步做基因层面的机制分析。
3.3 代谢分数与基因表达的联合分析
代谢分数只是表型，背后是基因表达的变化。找到与代谢分数高度相关的基因，可以揭示调控代谢的关键基因。做法是计算每个基因与代谢分数的 Spearman 相关系数，然后做排序，取 top 基因做富集分析。
# 计算基因与糖酵解分数的相关性 glycolysis_score <- seurat_obj@meta.data$HALLMARK_GLYCOLYSIS_Score1 expr_matrix <- GetAssayData(seurat_obj, slot = "data") # 对每个基因计算 Spearman 相关系数 cor_results <- apply(expr_matrix, 1, function(x) { if (sum(x > 0) < 10) return(c(0, 1)) # 表达细胞太少，跳过 ct <- cor.test(x, glycolysis_score, method = "spearman") return(c(ct$estimate, ct$p.value)) }) cor_df <- data.frame( gene = rownames(expr_matrix), rho = cor_results[1, ], pval = cor_results[2, ] ) cor_df$padj <- p.adjust(cor_df$pval, method = "BH") cor_df <- cor_df[order(abs(cor_df$rho), decreasing = TRUE), ] head(cor_df, 20)这段代码对每个基因做 Spearman 相关检验，apply遍历所有基因。sum(x > 0) < 10过滤掉表达细胞数少于 10 的基因，避免噪声。cor.test返回相关系数和 p 值，p.adjust做多重检验校正。最后按相关系数绝对值排序，取 top 20 基因。这些基因可能是代谢通路的直接成员，也可能是调控因子，需要结合注释判断。
4. 避坑与排查：小鼠单细胞代谢分析源码的五个血泪教训
4.1 基因 ID 不匹配导致打分全为零
现象：跑完AddModuleScore后，所有细胞的代谢分数都是零或接近零，FeaturePlot上没有任何颜色变化。
原因：基因集里的基因 Symbol 和 Seurat 对象里的基因名不一致。比如基因集里是Gapdh，但矩阵里是GAPDH或ENSMUSG00000057666。小鼠基因 Symbol 通常首字母大写、其余小写，但不同来源的数据可能用全大写或全小写。
解决：先检查基因集和矩阵的基因名交集。用intersect(names(metabolic_sets[[1]]), rownames(seurat_obj))看交集大小。如果交集很小，用toupper或tolower统一大小写，或者用bitr做 ID 转换。转换后重新跑打分。
4.2 背景基因数设置不当导致分数失真
现象：代谢分数在细胞类型之间没有差异，或者差异方向与预期相反。
原因：AddModuleScore的ctrl参数默认是 100，但如果目标基因集很大（比如 200 个基因），背景基因数太少会导致校正不充分。另外，nbin设置太小会让背景基因的表达分布和目标基因不匹配。
解决：把ctrl提高到 200 或 300，nbin保持在 24 到 30 之间。如果目标基因集超过 300 个基因，考虑拆分成子集分别打分，或者改用AUCell。跑完后用VlnPlot检查分数分布，正常应该是近似正态，如果出现双峰或长尾，说明校正有问题。
4.3 稀疏矩阵导致 AUCell 运行内存爆炸
现象：用AUCell给几万个细胞打分时，R 会话内存占用飙升，最后报cannot allocate vector of size错误。
原因：AUCell需要对每个细胞的基因表达排序，如果矩阵是稠密矩阵，内存占用是稀疏矩阵的几十倍。Seurat 默认的scale.data是稠密矩阵，直接传给AUCell会爆内存。
解决：用GetAssayData(seurat_obj, slot = "counts")提取稀疏矩阵，传给AUCell_buildRankings。如果还是不够，用AUCell的splitIntoBatches参数分批计算，每批 5000 个细胞。另外，提前用DietSeurat精简对象，去掉不需要的 assay 和降维结果。
4.4 批次效应未校正导致代谢分数假阳性
现象：不同样本的同一细胞类型代谢分数差异很大，但生物学分组之间没有差异。
原因：单细胞数据通常有批次效应，不同样本的测序深度、细胞活性、建库质量不同，导致代谢分数被批次主导。如果直接做组间比较，会把批次差异当成生物学差异。
解决：在打分前先做批次校正，用Harmony或Seurat的IntegrateLayers整合数据。打分后，用pseudobulk方法把细胞聚合成样本，再做组间比较。如果批次效应仍然明显，在limma模型里加入批次作为协变量。
4.5 代谢分数与细胞周期混淆
现象：增殖期细胞的糖酵解分数普遍偏高，导致分群时增殖细胞单独聚成一类。
原因：增殖细胞代谢活跃，糖酵解和氧化磷酸化相关基因表达上调，这是真实生物学现象，但会干扰细胞类型注释。如果研究目标不是增殖代谢，需要把细胞周期影响回归掉。
解决：用CellCycleScoring给每个细胞打细胞周期分数，然后在AddModuleScore之后，用scale函数对代谢分数做回归，把细胞周期分数作为协变量。或者，在差异分析时把细胞周期作为协变量纳入模型。如果增殖细胞是研究重点，则不需要回归，反而要单独分析。
5. 进阶技巧：用代谢分数做细胞亚群细分与轨迹推断
代谢分数不仅能做差异比较，还能用来细分细胞亚群。比如在肿瘤微环境里，同样注释为巨噬细胞的群体，糖酵解分数高的可能是 M1 样促炎表型，氧化磷酸化分数高的可能是 M2 样抑炎表型。做法是把代谢分数作为特征，和原来的基因表达矩阵一起做降维聚类。
# 把代谢分数加入降维特征 seurat_obj[["MetabolicScore"]] <- CreateAssayObject( data = t(seurat_obj@meta.data[, score_cols]) ) # 用代谢分数做 PCA seurat_obj <- RunPCA( seurat_obj, assay = "MetabolicScore", features = score_cols, reduction.name = "metabolic_pca", npcs = 5 ) # 用代谢 PCA 做 UMAP seurat_obj <- RunUMAP( seurat_obj, reduction = "metabolic_pca", dims = 1:5, reduction.name = "metabolic_umap" ) # 可视化 DimPlot(seurat_obj, reduction = "metabolic_umap", group.by = "cell_type")这段代码把代谢分数矩阵转成一个新的 Assay，然后用RunPCA做降维，reduction.name指定为metabolic_pca，避免和原来的 PCA 混淆。npcs = 5是因为代谢通路数量少，5 个主成分通常够用。RunUMAP基于代谢 PCA 做 UMAP，reduction.name指定为metabolic_umap。最后用DimPlot按细胞类型着色，看代谢分数能否把某些亚群分开。
如果代谢 UMAP 上出现了明显的亚群分离，可以进一步做轨迹推断。用monocle3或slingshot，把代谢 UMAP 的降维坐标作为输入，推断细胞在代谢状态之间的转换轨迹。比如从高糖酵解状态到高氧化磷酸化状态的转变，可能对应巨噬细胞的极化过程。
# 用 monocle3 做轨迹推断 library(monocle3) cds <- new_cell_data_set( expression_data = GetAssayData(seurat_obj, slot = "counts"), cell_metadata = seurat_obj@meta.data, gene_metadata = data.frame(gene_short_name = rownames(seurat_obj)) ) cds <- preprocess_cds(cds, method = "PCA", num_dim = 5) cds <- reduce_dimension(cds, reduction_method = "UMAP") cds <- cluster_cells(cds) cds <- learn_graph(cds) cds <- order_cells(cds) plot_cells(cds, color_cells_by = "HALLMARK_GLYCOLYSIS_Score1")new_cell_data_set创建 monocle3 对象，preprocess_cds做 PCA，reduce_dimension做 UMAP，cluster_cells聚类，learn_graph学习轨迹图，order_cells排序细胞。最后用plot_cells按糖酵解分数着色，看轨迹上的分数变化。如果轨迹起点是高糖酵解、终点是高氧化磷酸化，说明代谢状态转换方向明确。
我自己的习惯是，每次跑完代谢分析，先检查基因集交集，再跑打分，然后做 pseudobulk 差异分析，最后用代谢 UMAP 验证分群。这套流程跑下来，基本能避开大部分坑。代谢分析不像基因表达分析那么直接，但一旦跑通，能挖出很多常规分析看不到的生物学信息。希望帮到你。
本文还有配套的精品资源，点击获取
