简介面向单细胞转录组研究人员与生物信息学入门者这份源码包提供基于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 展开。维度AUCellssGSEA打分单位单个细胞单个细胞输入要求原始表达矩阵归一化矩阵受 dropout 影响弱基于排名中等交叉验证价值主方法备选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(objassays$RNAdata, plotStats TRUE)逻辑说明输入是归一化后的 data 矩阵不是 count也不是 scale.data。plotStats TRUE会输出每个细胞检测基因数的分布图方便你确认细胞质量。这一步计算量不小三万细胞大概要几分钟如果内存有限可以先随机抽两万细胞做测试确认通路打分逻辑没问题后再全量跑。rankings 对象是后面所有计算的中间产物建议用saveRDS存下来避免重复构建。3.2 AUCell_calcAUCaucMaxRank 是该调的第一个参数按基因集算 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(objassays$RNAdata) 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 把分数写回 SeuratAddMetaData 与 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.05UMAP 上没有任何梯度。原因最常见的是小鼠数据配了人类基因集。小鼠基因符号首字母大写后面小写如“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.9UMAP 上三张图几乎一模一样。原因基因集之间有重叠或者打分方法被总表达量主导。高表达基因多的细胞任何基因集排名都靠前本质上是细胞大小或测序深度效应。解决检查通路基因重叠长度如果重叠超过 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(objassays$RNAdata) 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% 区间之外才算站得住。这个办法对样本数少的数据尤其重要——三个样本组里有一个异常样本就足以制造假阳性。我自己吃过一次亏肝病模型里高代谢亚群与病理分期显著相关结果发现是测序深度差异导致的假象置换检验第一遍就跑出 p0.6白高兴两周。后来把置换检验写进标准化流程当成和质控图同级的固定环节不再跳过。做单细胞代谢分析分数只是起点能经受住随机背景拷问的差异才是值得写进结论的东西。希望这个习惯对你也有用。本文还有配套的精品资源点击获取