简介这份源码资源面向从事单细胞转录组与代谢研究的科研人员及生物信息学初学者围绕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_objmeta.data), value TRUE)AddModuleScore的核心逻辑是对目标基因集计算每个细胞的平均表达值然后从表达量相近的基因中随机抽取ctrl个背景基因计算背景平均表达两者相减得到校正分数。nbin 24表示把基因按平均表达量分成 24 个箱从同一箱里抽背景基因这样背景基因的表达分布和目标基因更接近校正更合理。ctrl 100是每个细胞抽 100 个背景基因这个值不能太小否则背景噪声大也不能太大否则计算慢。seed固定后结果可重复这在做差异分析时很重要。打分完成后seurat_objmeta.data里会多出几列列名是HALLMARK_GLYCOLYSIS_Score1这样的格式。注意Seurat 会自动在名字后面加数字如果多次运行同一个名字数字会递增。建议在循环里用paste0拼一个唯一名字避免混淆。2.4 打分结果的归一化和可视化打分出来不能直接用因为不同通路的分数范围不同有的通路分数在 -0.5 到 0.5 之间有的在 -1 到 1 之间。做跨通路比较时需要先归一化。我一般用 z-score 归一化把每个通路的分数转成标准正态分布这样不同通路之间可以横向比较。# 提取打分列 score_cols - grep(_Score, colnames(seurat_objmeta.data), value TRUE) # z-score 归一化 seurat_objmeta.data[score_cols] - scale(seurat_objmeta.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_objmeta.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里样本数少于 3limma的结果不稳定建议改用非参数检验或增加样本量。3.2 代谢通路之间的相关性分析代谢通路不是孤立的糖酵解和氧化磷酸化之间往往有代偿关系脂肪酸氧化和糖酵解也可能此消彼长。分析通路之间的相关性可以发现代谢重编程的模式。我一般会计算通路分数的 Spearman 相关系数然后做聚类热图。# 计算通路之间的 Spearman 相关性 cor_mat - cor(seurat_objmeta.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_objmeta.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 或 300nbin保持在 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_objmeta.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 做 UMAPreduction.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_objmeta.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做 PCAreduce_dimension做 UMAPcluster_cells聚类learn_graph学习轨迹图order_cells排序细胞。最后用plot_cells按糖酵解分数着色看轨迹上的分数变化。如果轨迹起点是高糖酵解、终点是高氧化磷酸化说明代谢状态转换方向明确。我自己的习惯是每次跑完代谢分析先检查基因集交集再跑打分然后做 pseudobulk 差异分析最后用代谢 UMAP 验证分群。这套流程跑下来基本能避开大部分坑。代谢分析不像基因表达分析那么直接但一旦跑通能挖出很多常规分析看不到的生物学信息。希望帮到你。本文还有配套的精品资源点击获取