简介面向从事单细胞转录组与代谢研究的科研人员及生物信息学学习者这套源码包聚焦解决scMetabolism包应用中的关键痛点——小鼠基因名向人类基因名的准确转换以及在不同Seuratv4/v5版本下的平滑集成。压缩包共含6个文件以两个R脚本搭建核心分析流程同时附有HTML可视化页面、Markdown说明文档与安装依赖脚本文件类型清晰、分工明确整个包仅9KB属于轻量级代码包便于快速下载与二次修改。透过源码可掌握从基因名转换、代谢激活分数计算到结果展示的完整步骤并配套安装配置说明及参考链接降低实际操作中的排错成本。目前已有171人学习使用适合希望以示例驱动方式快速上手单细胞代谢分析的研究者。1. 为什么单细胞转录组里要单独做代谢分析先聊一个很多人问过我的问题单细胞转录组本身就能看到基因表达为什么还要单独做代谢分析这一层原因是这样的。常规的差异表达、细胞分群、拟时序分析回答的是细胞是什么、从哪里来、要往哪去这一类问题。但代谢重编程这件事在单细胞层面是跨界的——代谢酶的转录水平不一定和代谢通路的通量直接挂钩而且代谢物本身是动态变化的纯看单个基因的表达量很容易漏掉真实的代谢状态变化。举个例子你在UMAP上看到一群细胞它们和相邻亚群在经典marker上其实差别不大但代谢偏好可能完全不一样一群偏糖酵解另一群偏脂肪酸氧化。这种差异在常规转录组层面会被淹没因为单个代谢酶的表达波动太小达不到差异阈值。小鼠单细胞代谢分析解决的正是这个问题把几十条甚至上百条代谢通路当作特征来打分然后基于这些通路活性对细胞进行重新分群、轨迹分析、组间比较。这也是代谢异质性这个概念的由来——代谢状态本身就可以作为定义细胞身份和功能状态的一个维度。这篇博客我会把我跑通的一套源码项目完整拆开讲从环境准备、核心分析流程、通路打分实现到结果可视化和我踩过的几个大坑。适用对象是已经会用Seurat做基本单细胞分析、但想在代谢维度上做深一层挖掘的朋友。如果你还没跑过单细胞分析建议先补一下标准的降维聚类流程再回来看这篇。2. 分析框架设计从表达矩阵到代谢通路活性说实话单细胞代谢分析没有官方唯一指定流程现在主流是两条路线一条是先做细胞分群注释再看各群代谢差异另一条是直接基于通路活性做无监督分群。我这次项目的做法是两条都走了一遍但主体以先分群、再打分、后比较为主线。整体框架大致是这样原始表达矩阵raw count → QC 归一化 降维 聚类 → 细胞类型注释经典marker → 代谢通路基因集的收集与整理 → 计算每个细胞在各通路上的活性分数 → 基于通路分数做差异分析 / 轨迹分析 / 组间比较 → 可视化气泡图、小提琴图、UMAP叠加分数、代谢轨迹这个设计思路的核心价值在于分层先锁定细胞身份再叠加代谢状态。如果跳过聚类注释直接做通路打分群与群之间的差异除了代谢状态可能还混入细胞类型本身的表达偏好结论很容易被污染。数据方面我用的是小鼠肝脏的单细胞数据Smart-seq2和10x都有测试因为肝脏代谢通路非常丰富糖代谢、脂代谢、氨基酸代谢都活跃用来验证流程最有说服力。当然了代码本身不绑定组织类型放到脑、肿瘤、肾脏数据上都可以跑只需要把基因集的物种对应关系改对就行。3. 核心源码拆解通路打分到底是怎么实现的3.1 基因集的准备与统一通路打分的第一步是拿到靠谱的代谢通路基因集。我优先用的就是KEGG里和代谢相关的通路配合Reactome的补充。如果你没有精力自己整理有个好消息——现在有不少现成的R包比如scMetabolism、VISION都是内置了代谢通路基因集的。但你千万别懒到完全不看基因集内容。我这次就踩过坑KEGG的小鼠基因注释用的是mmu前缀比如mmu00010糖酵解通路而有些包默认给你的是人类的hsa前缀。直接混用的话基因名匹配率低得感人。我的处理方式是统一做一次基因名转换以Ensembl ID或Symbol为准保证所有通路基因集和表达矩阵的基因名在同一种命名系统下。# 从msigdbr获取小鼠KEGG代谢通路基因集 library(msigdbr) library(dplyr) # msigdbr可以用species参数指定小鼠 mm_gs - msigdbr(species Mus musculus, category C2, subcategory CP:KEGG) # 筛代谢相关通路关键词匹配根据实际数据可以调整 metabolic_keywords - c(GLYCOLYSIS, OXIDATIVE_PHOSPHORYLATION, FATTY_ACID, PENTOSE, CITRATE, AMINO, PURINE, PYRIMIDINE, GLUTATHIONE, PORPHYRIN, TERPENOID) metabolic_gs - mm_gs %% filter(gs_name %in% grep(paste(metabolic_keywords, collapse |), gs_name, value TRUE)) %% select(gs_name, gene_symbol) %% distinct()这里有个细节要注意msigdbr的species Mus musculus会把基因映射成小鼠的gene symbol如果你拿到的表达矩阵里是Gapdh那这里也是Gapdh大小写一定要统一。我一般会再做一个toupper()以防止大小写不一致导致的匹配缺失。3.2 打分方法选型AddModuleScore vs AUCell vs ssGSEA实现了基因集之后剩下的核心问题就是怎么给每个细胞打出一个通路活性分数我在这个项目里对比了三种常用方法方法核心思路优点缺点Seurat AddModuleScore目标基因集平均表达减去背景基因平均表达快内置在Seurat方便背景基因抽样的随机性影响结果稳定性AUCell基于AUC评估基因集在细胞中的富集程度鲁棒性好适合做二元判断运行较慢需要设阈值ssGSEA单样本GSEA排序后统计富集单细胞上表现稳定适合比较对内存要求高实现稍复杂我的结论是如果只是做探索性分析直接用Seurat的AddModuleScore就够了因为它在速度和结果可解释性上最平衡。但如果你要做正式发表级别的比较建议同时用AUCell跑一版互相验证结论。我的建议是别只依赖一种方法。因为AddModuleScore对背景基因的选择比较敏感默认是随机抽取目标基因数量相同的基因作为背景遇到高表达基因占比大的通路比如氧化磷酸化分数会整体偏高。这时候你再用AUCell跑一遍如果两个方法的结果趋势一致才说明你的发现是稳健的。AUCell的大致用法library(AUCell) cells_rankings - AUCell_buildRankings(sce_obj, plotStats FALSE) cells_AUC - AUCell_calcAUC(metabolic_gs_list, cells_rankings) auc_matrix - getAUC(cells_AUC)然后你可以把这个auc_matrix直接塞回Seurat对象里当作一个新的assay来用后续的可视化、差异分析都是常规操作。3.3 通路活性分数的完整计算流程我跑通的一套核心代码大概是这样的逻辑这里只展示最关键的一段完整的源码项目里会有封装好的函数。# 假设seu是经过标准流程处理的Seurat对象metabolic_gs_list是通路基因集 # 方案一Seurat AddModuleScore seu - AddModuleScore(seu, features metabolic_gs_list, name MetaScore_) # 每一个通路的分数会存为MetaScore_1, MetaScore_2, ... # 你可以在meta.data里看到这些列这里我要特别强调一个操作细节AddModuleScore的features参数要求是一个list每个list元素是一组基因。如果你传入的是一个数据框而不是list它会报错或者只算第一个通路。另外基因名必须百分之百匹配表达矩阵的行名一个基因对不上整个通路分数都会有偏差。算完之后务必要做一个QC步骤看每个通路分数的分布。有些通路可能因为基因数太少少于5个打出来的分数方差极大这种通路要么删掉要么把基因数下限设到10以上再纳入分析。基因数太少的话分数基本没有统计说服力。4. 细胞代谢状态的比较差异通路筛选与组间分析4.1 差异代谢通路怎么筛拿到每个细胞在各条通路上的分数以后最常见的问题就变成了A组和B组的代谢状态到底有什么差异做法一般有两种一种是把细胞按分组做wilcox.test直接对通路分数做差异检验另一种是基于通路分数做置换检验控制混淆因素。第一种快第二种更严谨但计算量比较扛不住。我一般先用第一种筛出候选通路再用第二种验证关键结论。例如比较肿瘤组和对照组的糖酵解通路分数差异# 提取通路分数 score_col - grep(MetaScore_, colnames(seumeta.data), value TRUE) # 按分组做wilcox检验 df - data.frame( group seu$group, score seu$MetaScore_1 # 假设MetaScore_1是糖酵解通路 ) wilcox.test(score ~ group, data df)你需要特别注意多重检验校正。通路数多了以后P值很容易出现假阳性建议跑完所有通路之后统一做一次BH校正p.adjust只看校正后仍然显著的条目。4.2 代谢轨迹分析在拟时序里看代谢转变如果你的数据里有发育过程、疾病进展或者细胞状态转变强烈建议把代谢通路分数叠加到拟时序轨迹上看变化。这个操作其实非常直观——先跑Monocle3或者slingshot得到拟时序再把通路分数作为feature映射到轨迹上。我这次项目里做了一个比较有意思的结果肝脏的内皮细胞在沿着静息-活化轨迹变化的时候糖酵解相关通路是先升后降的而氧化磷酸化是持续上升的。单独看基因表达你很难发现这种代谢切换的规律但通路分数在轨迹上的分布一目了然。关键代码如下简化版# 假设已用monocle3构建好cds对象 # 把通路分数加入到cds的meta.data里 cdscolData$glycolysis_score - seu$MetaScore_1[colnames(cds)] # 用plot_cells可视化 plot_cells(cds, color_cells_by glycolysis_score, label_groups_by_cluster FALSE)这里有个容易栽的坑Monocle3的colData和Seurat的meta.data行名顺序不一定一致直接赋值可能导致错位。赋值之前务必先做行列名的匹配校验最好是用match()函数显式对齐而不是直接下标。这个坑我中过一次出来的图看似正常实际上分数全串位了。4.3 组间比较的可视化方案组间比较的可视化我推荐两件套小提琴图加气泡图。小提琴图用来展示单条通路在两个或多个组之间的分布差异适合展示具体通路的变化趋势气泡图用来展示多通路多组间的总体模式每个气泡大小代表显著比例或平均分数颜色代表方向上调/下调。气泡图参考实现library(ggplot2) # p_vals是通路差异检验P值avg_diff是平均分数差 plot_data - data.frame( pathway rownames(p_vals), pval -log10(p_vals$padj), diff avg_diff ) ggplot(plot_data, aes(x group_a, y pathway, size pval, color diff)) geom_point() scale_color_gradient2(low blue, mid white, high red) theme_minimal()我这里只是演示了绘图思路实际项目里我还会按通路大类糖代谢、脂代谢、氨基酸代谢、核苷酸代谢等分组加上分面这样从整体到局部都能看得很清楚。5. 小鼠数据特有的坑与避坑经验既然标题点明是小鼠这部分我单独拿出来说都是我自己实际踩过、重新跑过才算明白的经验。5.1 基因名大小写与版本问题小鼠基因symbol是首字母大写、其余小写比如Gapdh、Apoe。而很多重注释工具默认输入是人类的全部大写比如GAPDH。一旦你混用了这两套命名基因匹配率可能从90%直接掉到50%以下整个通路打分结果就是错的。我的建议是在分析的最早期统一基因名不要拖到打分阶段才考虑这个问题。用toupper()把表达矩阵行名和基因集都转成大写一次性消除大小写不一致如果你要保留小鼠原格式那就两边都保持原格式同时做好人工抽查。还有一个版本问题容易被忽略同一个基因在不同版本的小鼠基因组注释里symbol可能会变。比如以前叫GmXXXX的基因后面可能会被正式命名。如果你的表达矩阵来自较老的注释版本而基因集来自最新数据库就会出现部分基因匹配不上。解决办法是尽量用Ensembl ID做匹配最后再映射回symbol但这会增加工作量和踩坑概率新手建议直接保证两者来源版本不要差太远。5.2 线粒体基因的过滤阈值要重新跑线粒体基因比例过滤在人类数据里通常卡在20%但在小鼠组织里这个值不一定适用。比如小鼠肝脏组织的线粒体基因比例整体就偏高有些文献里也报到25%-30%以上如果机械地卡20%会过滤掉一大批真实的肝细胞。我这次项目里做QC时就发现肝细胞群的线粒体基因比例中位数就在22%左右如果按常规20%过滤肝细胞直接被削掉一大半后面所有分析都失真了。所以建议先画出线粒体基因比例的分布图再根据分布形状选拐点作为阈值而不是直接用别的文献里的数字。5.3 代谢基因集的物种匹配前面提到过msigdbr可以指定species Mus musculus来获取小鼠基因集。但如果你用的工具不直接支持物种切换比如某个Python包内置了人类KEGG基因集你要么手动做同源基因转换利用biomaRt的homolog映射要么干脆只保留那些在小鼠和人类之间同源性好、名称一致的基因。我个人的建议是优先选择支持物种参数的R包比如msigdbr、escape包在读取通路时会做物种匹配比你自己后期做同源映射要省心太多。5.4 组织特异性代谢特征会影响分群最后讲一个分析思路上的坑。小鼠不同组织的代谢基线差异极大肝脏以脂肪酸氧化和糖异生为主脑以葡萄糖氧化为主肌肉偏好脂肪酸和糖酵解混合。当你做代谢通路打分的时候同一个糖酵解通路在不同组织里的绝对分数没有可比性。所以如果你做的是多组织比较不要直接拿通路分数做跨组织的数值大小比较而是应该关注各自组织内部的组间差异方向和趋势。这次项目里我最核心的一个发现就是肝脏中不显著的代谢通路差异放到肿瘤相关成纤维细胞里却极其显著。这进一步说明代谢状态必须结合具体细胞类型和组织微环境来解读单纯跑了个高分意义有限。6. 可视化细节与图形表达怎么让审稿人一眼看到重点可视化本质上是对你分析结论的翻译。代谢通路数据量大、指标多如果一张图既要塞20个通路又要塞8个细胞类型必然糊成一团。我在多次迭代后形成了两个基本策略第一先宏观后微观。第一张图给全貌——UMAP或者分面气泡图把所有细胞类型和所有通路大类的关系展示出来第二张图聚焦关键的少数通路用小提琴图或者堆叠箱线图展示组间或群间差异第三张图给轨迹或者通路流图展示代谢状态的变化路径。这三层递进读者跟着你的图走一遍自然就能还原你的分析逻辑。第二颜色编码要有意义。用在代谢通路分数展示上的连续渐变颜色我建议用blue-white-red这类发散色阶中间色对应0分或中位数这样哪些通路偏高哪些偏低一眼可见。避免用那些从暗到亮不连续的自定义色板会让读者把注意力放在颜色差异上而不是数据模式上。我这次项目里最后敲定了两个核心展示方案一个是用ggplot2画的多组气泡图另一个是把代谢分数映射到Monocle3轨迹上的图。每一张图的caption里都写清楚了用的是哪种打分方法AddModuleScore还是AUCell以及用了哪个版本的基因集。这一点在发文章或者写技术报告时非常重要很多细节影响可重复性。7. 源码项目的目录结构与管理方式最后说下源码项目的组织方式因为很多人拿到代码后跑不通一半原因出在源码结构不清楚而不是代码本身有问题。我的项目目录是这样的mouse_single_cell_metabolism/ ├── 01_qc_and_clustering.R # QC、降维、聚类 ├── 02_cell_type_annotation.R # 细胞类型注释 ├── 03_metabolic_scoring.R # 通路打分核心代码 ├── 04_differential_metabolism.R # 组间差异分析 ├── 05_trajectory_analysis.R # 拟时序与代谢轨迹 ├── 06_visualization.R # 所有图和表输出 ├── data/ │ ├── raw/ # 原始表达矩阵 │ └── processed/ # 分析中间产物 ├── results/ │ ├── tables/ # 差异表、打分表 │ └── figures/ # 输出图 └── README.md每个脚本我都在开始部分写好输入是什么、输出是什么、依赖哪些上游文件的产物这样可以保证你按1到6的顺序执行就能完整复现。实际运行中很容易因为中间产物路径不对导致下游脚本报错这也是我踩过的最浪费时间的坑。另外强烈建议用renv或conda锁定环境版本。单细胞分析工具链的版本兼容性极差Seurat 4和Seurat 5之间有些API直接变了Monocle3又依赖特定的R版本。我这次把sessionInfo()记录到了results/session_info.txt里方便后续回来看清依赖版本也方便别人复现时对齐环境。花点时间把源码组织好短期看是增加了工作量长期看给你节约的时间是数倍的尤其是当你三四个月后要回来看自己写的代码时。本文还有配套的精品资源点击获取
小鼠单细胞代谢分析实战：通路打分与可视化全流程解析
简介面向从事单细胞转录组与代谢研究的科研人员及生物信息学学习者这套源码包聚焦解决scMetabolism包应用中的关键痛点——小鼠基因名向人类基因名的准确转换以及在不同Seuratv4/v5版本下的平滑集成。压缩包共含6个文件以两个R脚本搭建核心分析流程同时附有HTML可视化页面、Markdown说明文档与安装依赖脚本文件类型清晰、分工明确整个包仅9KB属于轻量级代码包便于快速下载与二次修改。透过源码可掌握从基因名转换、代谢激活分数计算到结果展示的完整步骤并配套安装配置说明及参考链接降低实际操作中的排错成本。目前已有171人学习使用适合希望以示例驱动方式快速上手单细胞代谢分析的研究者。1. 为什么单细胞转录组里要单独做代谢分析先聊一个很多人问过我的问题单细胞转录组本身就能看到基因表达为什么还要单独做代谢分析这一层原因是这样的。常规的差异表达、细胞分群、拟时序分析回答的是细胞是什么、从哪里来、要往哪去这一类问题。但代谢重编程这件事在单细胞层面是跨界的——代谢酶的转录水平不一定和代谢通路的通量直接挂钩而且代谢物本身是动态变化的纯看单个基因的表达量很容易漏掉真实的代谢状态变化。举个例子你在UMAP上看到一群细胞它们和相邻亚群在经典marker上其实差别不大但代谢偏好可能完全不一样一群偏糖酵解另一群偏脂肪酸氧化。这种差异在常规转录组层面会被淹没因为单个代谢酶的表达波动太小达不到差异阈值。小鼠单细胞代谢分析解决的正是这个问题把几十条甚至上百条代谢通路当作特征来打分然后基于这些通路活性对细胞进行重新分群、轨迹分析、组间比较。这也是代谢异质性这个概念的由来——代谢状态本身就可以作为定义细胞身份和功能状态的一个维度。这篇博客我会把我跑通的一套源码项目完整拆开讲从环境准备、核心分析流程、通路打分实现到结果可视化和我踩过的几个大坑。适用对象是已经会用Seurat做基本单细胞分析、但想在代谢维度上做深一层挖掘的朋友。如果你还没跑过单细胞分析建议先补一下标准的降维聚类流程再回来看这篇。2. 分析框架设计从表达矩阵到代谢通路活性说实话单细胞代谢分析没有官方唯一指定流程现在主流是两条路线一条是先做细胞分群注释再看各群代谢差异另一条是直接基于通路活性做无监督分群。我这次项目的做法是两条都走了一遍但主体以先分群、再打分、后比较为主线。整体框架大致是这样原始表达矩阵raw count → QC 归一化 降维 聚类 → 细胞类型注释经典marker → 代谢通路基因集的收集与整理 → 计算每个细胞在各通路上的活性分数 → 基于通路分数做差异分析 / 轨迹分析 / 组间比较 → 可视化气泡图、小提琴图、UMAP叠加分数、代谢轨迹这个设计思路的核心价值在于分层先锁定细胞身份再叠加代谢状态。如果跳过聚类注释直接做通路打分群与群之间的差异除了代谢状态可能还混入细胞类型本身的表达偏好结论很容易被污染。数据方面我用的是小鼠肝脏的单细胞数据Smart-seq2和10x都有测试因为肝脏代谢通路非常丰富糖代谢、脂代谢、氨基酸代谢都活跃用来验证流程最有说服力。当然了代码本身不绑定组织类型放到脑、肿瘤、肾脏数据上都可以跑只需要把基因集的物种对应关系改对就行。3. 核心源码拆解通路打分到底是怎么实现的3.1 基因集的准备与统一通路打分的第一步是拿到靠谱的代谢通路基因集。我优先用的就是KEGG里和代谢相关的通路配合Reactome的补充。如果你没有精力自己整理有个好消息——现在有不少现成的R包比如scMetabolism、VISION都是内置了代谢通路基因集的。但你千万别懒到完全不看基因集内容。我这次就踩过坑KEGG的小鼠基因注释用的是mmu前缀比如mmu00010糖酵解通路而有些包默认给你的是人类的hsa前缀。直接混用的话基因名匹配率低得感人。我的处理方式是统一做一次基因名转换以Ensembl ID或Symbol为准保证所有通路基因集和表达矩阵的基因名在同一种命名系统下。# 从msigdbr获取小鼠KEGG代谢通路基因集 library(msigdbr) library(dplyr) # msigdbr可以用species参数指定小鼠 mm_gs - msigdbr(species Mus musculus, category C2, subcategory CP:KEGG) # 筛代谢相关通路关键词匹配根据实际数据可以调整 metabolic_keywords - c(GLYCOLYSIS, OXIDATIVE_PHOSPHORYLATION, FATTY_ACID, PENTOSE, CITRATE, AMINO, PURINE, PYRIMIDINE, GLUTATHIONE, PORPHYRIN, TERPENOID) metabolic_gs - mm_gs %% filter(gs_name %in% grep(paste(metabolic_keywords, collapse |), gs_name, value TRUE)) %% select(gs_name, gene_symbol) %% distinct()这里有个细节要注意msigdbr的species Mus musculus会把基因映射成小鼠的gene symbol如果你拿到的表达矩阵里是Gapdh那这里也是Gapdh大小写一定要统一。我一般会再做一个toupper()以防止大小写不一致导致的匹配缺失。3.2 打分方法选型AddModuleScore vs AUCell vs ssGSEA实现了基因集之后剩下的核心问题就是怎么给每个细胞打出一个通路活性分数我在这个项目里对比了三种常用方法方法核心思路优点缺点Seurat AddModuleScore目标基因集平均表达减去背景基因平均表达快内置在Seurat方便背景基因抽样的随机性影响结果稳定性AUCell基于AUC评估基因集在细胞中的富集程度鲁棒性好适合做二元判断运行较慢需要设阈值ssGSEA单样本GSEA排序后统计富集单细胞上表现稳定适合比较对内存要求高实现稍复杂我的结论是如果只是做探索性分析直接用Seurat的AddModuleScore就够了因为它在速度和结果可解释性上最平衡。但如果你要做正式发表级别的比较建议同时用AUCell跑一版互相验证结论。我的建议是别只依赖一种方法。因为AddModuleScore对背景基因的选择比较敏感默认是随机抽取目标基因数量相同的基因作为背景遇到高表达基因占比大的通路比如氧化磷酸化分数会整体偏高。这时候你再用AUCell跑一遍如果两个方法的结果趋势一致才说明你的发现是稳健的。AUCell的大致用法library(AUCell) cells_rankings - AUCell_buildRankings(sce_obj, plotStats FALSE) cells_AUC - AUCell_calcAUC(metabolic_gs_list, cells_rankings) auc_matrix - getAUC(cells_AUC)然后你可以把这个auc_matrix直接塞回Seurat对象里当作一个新的assay来用后续的可视化、差异分析都是常规操作。3.3 通路活性分数的完整计算流程我跑通的一套核心代码大概是这样的逻辑这里只展示最关键的一段完整的源码项目里会有封装好的函数。# 假设seu是经过标准流程处理的Seurat对象metabolic_gs_list是通路基因集 # 方案一Seurat AddModuleScore seu - AddModuleScore(seu, features metabolic_gs_list, name MetaScore_) # 每一个通路的分数会存为MetaScore_1, MetaScore_2, ... # 你可以在meta.data里看到这些列这里我要特别强调一个操作细节AddModuleScore的features参数要求是一个list每个list元素是一组基因。如果你传入的是一个数据框而不是list它会报错或者只算第一个通路。另外基因名必须百分之百匹配表达矩阵的行名一个基因对不上整个通路分数都会有偏差。算完之后务必要做一个QC步骤看每个通路分数的分布。有些通路可能因为基因数太少少于5个打出来的分数方差极大这种通路要么删掉要么把基因数下限设到10以上再纳入分析。基因数太少的话分数基本没有统计说服力。4. 细胞代谢状态的比较差异通路筛选与组间分析4.1 差异代谢通路怎么筛拿到每个细胞在各条通路上的分数以后最常见的问题就变成了A组和B组的代谢状态到底有什么差异做法一般有两种一种是把细胞按分组做wilcox.test直接对通路分数做差异检验另一种是基于通路分数做置换检验控制混淆因素。第一种快第二种更严谨但计算量比较扛不住。我一般先用第一种筛出候选通路再用第二种验证关键结论。例如比较肿瘤组和对照组的糖酵解通路分数差异# 提取通路分数 score_col - grep(MetaScore_, colnames(seumeta.data), value TRUE) # 按分组做wilcox检验 df - data.frame( group seu$group, score seu$MetaScore_1 # 假设MetaScore_1是糖酵解通路 ) wilcox.test(score ~ group, data df)你需要特别注意多重检验校正。通路数多了以后P值很容易出现假阳性建议跑完所有通路之后统一做一次BH校正p.adjust只看校正后仍然显著的条目。4.2 代谢轨迹分析在拟时序里看代谢转变如果你的数据里有发育过程、疾病进展或者细胞状态转变强烈建议把代谢通路分数叠加到拟时序轨迹上看变化。这个操作其实非常直观——先跑Monocle3或者slingshot得到拟时序再把通路分数作为feature映射到轨迹上。我这次项目里做了一个比较有意思的结果肝脏的内皮细胞在沿着静息-活化轨迹变化的时候糖酵解相关通路是先升后降的而氧化磷酸化是持续上升的。单独看基因表达你很难发现这种代谢切换的规律但通路分数在轨迹上的分布一目了然。关键代码如下简化版# 假设已用monocle3构建好cds对象 # 把通路分数加入到cds的meta.data里 cdscolData$glycolysis_score - seu$MetaScore_1[colnames(cds)] # 用plot_cells可视化 plot_cells(cds, color_cells_by glycolysis_score, label_groups_by_cluster FALSE)这里有个容易栽的坑Monocle3的colData和Seurat的meta.data行名顺序不一定一致直接赋值可能导致错位。赋值之前务必先做行列名的匹配校验最好是用match()函数显式对齐而不是直接下标。这个坑我中过一次出来的图看似正常实际上分数全串位了。4.3 组间比较的可视化方案组间比较的可视化我推荐两件套小提琴图加气泡图。小提琴图用来展示单条通路在两个或多个组之间的分布差异适合展示具体通路的变化趋势气泡图用来展示多通路多组间的总体模式每个气泡大小代表显著比例或平均分数颜色代表方向上调/下调。气泡图参考实现library(ggplot2) # p_vals是通路差异检验P值avg_diff是平均分数差 plot_data - data.frame( pathway rownames(p_vals), pval -log10(p_vals$padj), diff avg_diff ) ggplot(plot_data, aes(x group_a, y pathway, size pval, color diff)) geom_point() scale_color_gradient2(low blue, mid white, high red) theme_minimal()我这里只是演示了绘图思路实际项目里我还会按通路大类糖代谢、脂代谢、氨基酸代谢、核苷酸代谢等分组加上分面这样从整体到局部都能看得很清楚。5. 小鼠数据特有的坑与避坑经验既然标题点明是小鼠这部分我单独拿出来说都是我自己实际踩过、重新跑过才算明白的经验。5.1 基因名大小写与版本问题小鼠基因symbol是首字母大写、其余小写比如Gapdh、Apoe。而很多重注释工具默认输入是人类的全部大写比如GAPDH。一旦你混用了这两套命名基因匹配率可能从90%直接掉到50%以下整个通路打分结果就是错的。我的建议是在分析的最早期统一基因名不要拖到打分阶段才考虑这个问题。用toupper()把表达矩阵行名和基因集都转成大写一次性消除大小写不一致如果你要保留小鼠原格式那就两边都保持原格式同时做好人工抽查。还有一个版本问题容易被忽略同一个基因在不同版本的小鼠基因组注释里symbol可能会变。比如以前叫GmXXXX的基因后面可能会被正式命名。如果你的表达矩阵来自较老的注释版本而基因集来自最新数据库就会出现部分基因匹配不上。解决办法是尽量用Ensembl ID做匹配最后再映射回symbol但这会增加工作量和踩坑概率新手建议直接保证两者来源版本不要差太远。5.2 线粒体基因的过滤阈值要重新跑线粒体基因比例过滤在人类数据里通常卡在20%但在小鼠组织里这个值不一定适用。比如小鼠肝脏组织的线粒体基因比例整体就偏高有些文献里也报到25%-30%以上如果机械地卡20%会过滤掉一大批真实的肝细胞。我这次项目里做QC时就发现肝细胞群的线粒体基因比例中位数就在22%左右如果按常规20%过滤肝细胞直接被削掉一大半后面所有分析都失真了。所以建议先画出线粒体基因比例的分布图再根据分布形状选拐点作为阈值而不是直接用别的文献里的数字。5.3 代谢基因集的物种匹配前面提到过msigdbr可以指定species Mus musculus来获取小鼠基因集。但如果你用的工具不直接支持物种切换比如某个Python包内置了人类KEGG基因集你要么手动做同源基因转换利用biomaRt的homolog映射要么干脆只保留那些在小鼠和人类之间同源性好、名称一致的基因。我个人的建议是优先选择支持物种参数的R包比如msigdbr、escape包在读取通路时会做物种匹配比你自己后期做同源映射要省心太多。5.4 组织特异性代谢特征会影响分群最后讲一个分析思路上的坑。小鼠不同组织的代谢基线差异极大肝脏以脂肪酸氧化和糖异生为主脑以葡萄糖氧化为主肌肉偏好脂肪酸和糖酵解混合。当你做代谢通路打分的时候同一个糖酵解通路在不同组织里的绝对分数没有可比性。所以如果你做的是多组织比较不要直接拿通路分数做跨组织的数值大小比较而是应该关注各自组织内部的组间差异方向和趋势。这次项目里我最核心的一个发现就是肝脏中不显著的代谢通路差异放到肿瘤相关成纤维细胞里却极其显著。这进一步说明代谢状态必须结合具体细胞类型和组织微环境来解读单纯跑了个高分意义有限。6. 可视化细节与图形表达怎么让审稿人一眼看到重点可视化本质上是对你分析结论的翻译。代谢通路数据量大、指标多如果一张图既要塞20个通路又要塞8个细胞类型必然糊成一团。我在多次迭代后形成了两个基本策略第一先宏观后微观。第一张图给全貌——UMAP或者分面气泡图把所有细胞类型和所有通路大类的关系展示出来第二张图聚焦关键的少数通路用小提琴图或者堆叠箱线图展示组间或群间差异第三张图给轨迹或者通路流图展示代谢状态的变化路径。这三层递进读者跟着你的图走一遍自然就能还原你的分析逻辑。第二颜色编码要有意义。用在代谢通路分数展示上的连续渐变颜色我建议用blue-white-red这类发散色阶中间色对应0分或中位数这样哪些通路偏高哪些偏低一眼可见。避免用那些从暗到亮不连续的自定义色板会让读者把注意力放在颜色差异上而不是数据模式上。我这次项目里最后敲定了两个核心展示方案一个是用ggplot2画的多组气泡图另一个是把代谢分数映射到Monocle3轨迹上的图。每一张图的caption里都写清楚了用的是哪种打分方法AddModuleScore还是AUCell以及用了哪个版本的基因集。这一点在发文章或者写技术报告时非常重要很多细节影响可重复性。7. 源码项目的目录结构与管理方式最后说下源码项目的组织方式因为很多人拿到代码后跑不通一半原因出在源码结构不清楚而不是代码本身有问题。我的项目目录是这样的mouse_single_cell_metabolism/ ├── 01_qc_and_clustering.R # QC、降维、聚类 ├── 02_cell_type_annotation.R # 细胞类型注释 ├── 03_metabolic_scoring.R # 通路打分核心代码 ├── 04_differential_metabolism.R # 组间差异分析 ├── 05_trajectory_analysis.R # 拟时序与代谢轨迹 ├── 06_visualization.R # 所有图和表输出 ├── data/ │ ├── raw/ # 原始表达矩阵 │ └── processed/ # 分析中间产物 ├── results/ │ ├── tables/ # 差异表、打分表 │ └── figures/ # 输出图 └── README.md每个脚本我都在开始部分写好输入是什么、输出是什么、依赖哪些上游文件的产物这样可以保证你按1到6的顺序执行就能完整复现。实际运行中很容易因为中间产物路径不对导致下游脚本报错这也是我踩过的最浪费时间的坑。另外强烈建议用renv或conda锁定环境版本。单细胞分析工具链的版本兼容性极差Seurat 4和Seurat 5之间有些API直接变了Monocle3又依赖特定的R版本。我这次把sessionInfo()记录到了results/session_info.txt里方便后续回来看清依赖版本也方便别人复现时对齐环境。花点时间把源码组织好短期看是增加了工作量长期看给你节约的时间是数倍的尤其是当你三四个月后要回来看自己写的代码时。本文还有配套的精品资源点击获取
