简介：这份源码资源面向从事单细胞转录组与代谢研究的生信分析人员及R语言学习者，聚焦小鼠单细胞代谢激活分数分析这一具体场景，解决从基因表达数据出发、借助scMetabolism包完成代谢通路打分并适配Seurat v4/v5版本的实际问题。资源包共6个文件，以R脚本为主要类型，包含代谢分析主流程脚本与依赖包安装脚本，另附README说明文档、HTML页面及配置文件，压缩包整体约9KB，体量轻便、结构清晰，便于直接运行与二次修改。目前已有176人学习下载。读者可从中获得一套可复用的分析代码，掌握小鼠基因名向人类基因名转换以对接人类代谢通路数据库的关键处理思路，并了解如何将代谢打分结果导入Seurat进行后续降维聚类与可视化，适合作为单细胞代谢研究的入门参考与排错对照。
1. 小鼠单细胞代谢分析源码：从表达矩阵到代谢通路的完整拆解
单细胞测序做到代谢层面，很多人第一反应是“表达矩阵都拿到了，代谢还能难到哪去”。真上手才发现，代谢基因本身表达量低、稀疏严重，常规的 Seurat 流程跑完聚类，代谢通路那一步基本是空的。这份小鼠单细胞代谢分析源码，解决的正是这个断层——它把单细胞表达矩阵到代谢通路活性评分之间的链路补全了，包含数据预处理、代谢基因集匹配、通路活性打分、细胞亚群代谢异质性比较几个模块。适合已经跑通过单细胞基础流程、想往代谢方向延伸的从业者，也适合做肿瘤微环境、免疫代谢、肝脏代谢这类课题的研究生。源码是 Python 写的，依赖 scanpy 和 scipy，不绑定特定平台，拿到矩阵就能跑。
2. 代谢分析源码的环境搭建与数据准备：scanpy 版本和矩阵格式的硬约束
2.1 为什么选 scanpy 而不是 Seurat 做代谢打分
单细胞代谢分析的核心操作是“按基因集给每个细胞打分”，这件事在 R 里用 Seurat 的AddModuleScore也能做，但代谢基因集动辄几百个基因，跨物种同源基因映射在 R 里要额外维护一个 ortholog 表，容易出错。scanpy 的score_genes直接吃 Python 字典格式的基因集，配合mygene做小鼠-人类基因符号转换，链路更短。源码里用的是 KEGG 代谢通路基因集，也预留了 Reactome 和 GO 代谢相关条目的接口。另一个实际原因是内存——小鼠单细胞数据动辄几万个细胞，scanpy 的稀疏矩阵处理比 Seurat 在同等内存下能多扛约 30% 的细胞数，这个数字是我在 16G 内存机器上反复跑出来的经验值，不是官方 benchmark。
2.2 环境依赖与安装步骤
源码根目录下有一个requirements.txt，但直接pip install -r大概率会在 scanpy 版本上翻车。我建议按下面的顺序手动装，版本号是源码注释里标明的兼容区间：
# 先建独立环境，避免和已有 scanpy 冲突 conda create -n sc_metab python=3.9 -y conda activate sc_metab # scanpy 锁在 1.9.x，1.10 之后 score_genes 的默认参数变了 pip install scanpy==1.9.6 pip install mygene==3.2.2 pip install scipy==1.10.1 pip install pandas==1.5.3 pip install matplotlib==3.7.1 pip install seaborn==0.12.2 # 验证核心依赖 python -c "import scanpy as sc; print(sc.__version__)"逻辑说明：scanpy 1.9.6 的score_genes在ctrl_size参数默认值上和 1.10 不同，源码里的打分阈值是基于 1.9.x 调的，换版本会导致同一份数据打分结果偏移。mygene 用来做基因符号转换，小鼠基因符号首字母大写、人类全大写，不转换的话 KEGG 基因集匹配率会掉到 40% 以下。scipy 锁 1.10.1 是因为 1.11 改了稀疏矩阵的eliminate_zeros行为，会影响代谢基因稀疏过滤那一步。
2.3 输入矩阵的格式要求与预处理
源码接受三种输入：10x Genomics 的filtered_feature_bc_matrix目录、h5ad 文件、以及 CSV 格式的稀疏矩阵三件套（matrix.mtx + genes.tsv + barcodes.tsv）。我一般直接用 h5ad，因为前一步的质控和归一化已经在 scanpy 里做完了。如果你的矩阵还带着原始 counts，源码里的preprocess.py会走一遍标准流程：
import scanpy as sc import numpy as np # 读入原始矩阵 adata = sc.read_10x_mtx("filtered_feature_bc_matrix/", var_names="gene_symbols") # 基础质控：线粒体基因比例过滤 adata.var["mt"] = adata.var_names.str.startswith("mt-") sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True) adata = adata[adata.obs["pct_counts_mt"] < 20, :].copy() adata = adata[adata.obs["n_genes_by_counts"] > 200, :].copy() # 归一化与对数化 sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) # 高变基因与降维（代谢分析前必须做，否则聚类不准） sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor="seurat") adata = adata[:, adata.var["highly_variable"]].copy() sc.pp.scale(adata, max_value=10) sc.tl.pca(adata, n_comps=30) sc.pp.neighbors(adata, n_neighbors=15, n_pcs=30) sc.tl.leiden(adata, resolution=0.8)参数说明：pct_counts_mt < 20是小鼠数据的经验阈值，人类数据一般卡 10，小鼠线粒体基因少、比例天然偏高，卡太严会丢掉真实细胞。n_top_genes=2000是 scanpy 教程的默认值，但代谢分析里我建议提到 3000，因为代谢基因本身表达低，2000 个高变基因里可能只覆盖到几十个代谢基因。resolution=0.8是聚类分辨率，做代谢异质性比较时如果亚群分得太粗，代谢差异会被平均掉，我一般会试 0.6 到 1.2 三个值，看哪个分辨率下代谢通路打分在亚群间有显著差异。
注意：如果你的数据来自 Smart-seq2 而不是 10x，
normalize_total的target_sum要改成 1e6，因为 Smart-seq2 的测序深度和 10x 不是一个量级，用 1e4 会把低表达代谢基因压没。
3. 代谢基因集匹配与通路活性打分：KEGG 基因集怎么落到小鼠基因符号上
3.1 基因符号转换的坑与 mygene 批量查询
KEGG 的代谢通路基因集默认是人类基因符号，小鼠基因符号只有首字母大写这一个区别，听起来简单，但实际数据里混着 Ensembl ID、旧版符号、别名，直接大写转换会漏掉一批。源码里用 mygene 做批量查询，把小鼠符号映射到人类 Entrez ID 再取人类符号，匹配率能从 60% 提到 85% 以上。这一步的代码在gene_mapping.py里：
import mygene import pandas as pd mg = mygene.MyGeneInfo() def map_mouse_to_human(mouse_genes): """把小鼠基因符号批量映射到人类符号""" # 去掉线粒体基因和未命名基因 mouse_genes = [g for g in mouse_genes if not g.startswith("mt-") and g != ""] # 批量查询，scopes 指定小鼠符号，fields 取人类同源符号 results = mg.querymany( mouse_genes, scopes="symbol", species="mouse", fields="homologene, symbol", as_dataframe=True ) # 取 homologene 里的人类同源符号 mapping = {} for idx, row in results.iterrows(): if isinstance(row.get("homologene"), dict): human_sym = row["homologene"].get("human", {}).get("symbol") if human_sym: mapping[idx] = human_sym return mapping逻辑说明：querymany的scopes="symbol"表示输入是小鼠基因符号，species="mouse"限定物种，fields="homologene"取同源基因信息。返回的 DataFrame 里homologene字段是一个嵌套字典，人类同源符号在human.symbol路径下。这个查询走的是 NCBI 的接口，一次查几千个基因大概要 2 到 3 分钟，源码里加了缓存机制，第二次跑同一批基因直接读本地 pickle，不用重复请求。
参数说明：as_dataframe=True让返回结果直接是 DataFrame，方便后续遍历。如果查询的基因里有别名或旧符号，mygene 会自动做模糊匹配，但会在notfound字段里标记出来，源码里对notfound的基因做了二次查询，用scopes="alias"再试一次，能再捞回来 5% 左右。
3.2 KEGG 代谢通路基因集的加载与过滤
源码的pathways/目录下有一个kegg_metabolism.json，里面是手动整理过的 KEGG 代谢通路基因集，去掉了疾病通路和信号通路，只保留碳水化合物代谢、脂质代谢、氨基酸代谢、能量代谢这几大类。加载和过滤的代码：
import json def load_metabolic_pathways(path="pathways/kegg_metabolism.json", min_genes=10): """加载代谢通路基因集，过滤掉基因数太少的通路""" with open(path, "r") as f: pathways = json.load(f) # 过滤：通路基因数少于 min_genes 的丢掉 filtered = {} for name, genes in pathways.items(): if len(genes) >= min_genes: filtered[name] = genes print(f"原始通路数: {len(pathways)}, 过滤后: {len(filtered)}") return filtered逻辑说明：min_genes=10是个经验值，KEGG 里有些通路只有三五个基因，打分结果噪声极大，不如直接丢掉。源码里默认保留约 80 条代谢通路，覆盖糖酵解、TCA 循环、氧化磷酸化、脂肪酸氧化这些核心代谢过程。
参数说明：如果你做的是特定方向，比如只关注脂质代谢，可以把kegg_metabolism.json里其他类别的通路删掉，减少多重检验校正的负担。min_genes可以调到 15，但再高就会把一些有生物学意义的短通路也滤掉，我试过 20，糖酵解通路只剩 8 个基因，打分结果和预期完全对不上。
3.3 用 score_genes 做通路活性打分
打分这一步是整份源码的核心。scanpy 的score_genes原理是：对每个细胞，计算通路基因集的平均表达量，再减去随机背景基因集的平均表达量，得到相对活性分数。源码里对每个通路循环打分，结果存到adata.obs里：
import scanpy as sc def score_pathways(adata, pathways, ctrl_size=50): """对每个代谢通路打分，结果写入 adata.obs""" for name, genes in pathways.items(): # 只保留在 adata 里存在的基因 valid_genes = [g for g in genes if g in adata.var_names] if len(valid_genes) < 5: continue score_name = f"metab_{name}" sc.tl.score_genes( adata, gene_list=valid_genes, ctrl_size=ctrl_size, score_name=score_name, random_state=42 ) return adata逻辑说明：valid_genes过滤掉数据里不存在的基因，如果通路里有效基因少于 5 个就跳过，避免噪声。ctrl_size=50是背景基因集的基因数，源码默认 50，我试过 100，打分结果更平滑但亚群间差异也被抹平了，50 是个平衡点。random_state=42固定随机种子，保证每次跑结果一致，这个在做多次分析对比时很重要。
参数说明：score_name前缀metab_是为了和后续其他打分区分开。打分结果是一个 Z-score 形式的相对值，不是绝对活性，所以跨数据集比较时要注意——不同数据集的背景基因集不同，同一个通路的打分值不能直接比大小，只能在同一数据集内部比较亚群间的相对高低。
4. 细胞亚群代谢异质性比较与可视化：从打分矩阵到差异通路
4.1 亚群间代谢通路差异的统计检验
打分跑完后，adata.obs里每个细胞每个通路都有一个分数。下一步是看哪些通路在不同亚群间有显著差异。源码里用 Mann-Whitney U 检验做两两比较，再用 Benjamini-Hochberg 做多重检验校正：
from scipy.stats import mannwhitneyu from statsmodels.stats.multitest import multipletests import pandas as pd def compare_pathways(adata, group_key="leiden", score_prefix="metab_"): """比较各亚群间代谢通路打分的差异""" score_cols = [c for c in adata.obs.columns if c.startswith(score_prefix)] groups = adata.obs[group_key].unique() results = [] for pathway in score_cols: for i, g1 in enumerate(groups): for g2 in groups[i+1:]: s1 = adata.obs.loc[adata.obs[group_key] == g1, pathway] s2 = adata.obs.loc[adata.obs[group_key] == g2, pathway] stat, pval = mannwhitneyu(s1, s2, alternative="two-sided") results.append({ "pathway": pathway.replace(score_prefix, ""), "group1": g1, "group2": g2, "pval": pval, "median_diff": s1.median() - s2.median() }) df = pd.DataFrame(results) # BH 校正 df["padj"] = multipletests(df["pval"], method="fdr_bh")[1] return df[df["padj"] < 0.05].sort_values("padj")逻辑说明：Mann-Whitney U 检验不假设正态分布，适合单细胞打分这种偏态数据。alternative="two-sided"做双尾检验，因为我们不确定哪个亚群代谢活性更高。median_diff记录中位数差值，正负号表示方向。BH 校正用fdr_bh方法，这是单细胞分析里的标准做法，控制假发现率在 5% 以下。
参数说明：如果亚群数超过 10 个，两两比较的组合数会爆炸，multipletests的校正会非常严格，很多真实差异会被滤掉。这时候我建议先做整体 Kruskal-Wallis 检验，只对整体显著的 pathway 做两两比较，源码里预留了这个开关，把compare_pathways的pre_filter参数设为True就会先跑 Kruskal-Wallis。
4.2 代谢通路热图与气泡图的绘制
差异结果出来后，可视化是给合作者看的关键。源码里提供了两个函数：plot_pathway_heatmap画亚群-通路热图，plot_pathway_dotplot画气泡图。热图用 seaborn 的 clustermap，气泡图用 matplotlib 手搓：
import seaborn as sns import matplotlib.pyplot as plt import numpy as np def plot_pathway_heatmap(adata, pathways, group_key="leiden", figsize=(12, 8)): """画亚群-代谢通路热图""" # 按亚群取每个通路的平均打分 score_cols = [f"metab_{p}" for p in pathways if f"metab_{p}" in adata.obs.columns] mean_scores = adata.obs.groupby(group_key)[score_cols].mean() # Z-score 标准化，让通路间可比 mean_scores = (mean_scores - mean_scores.mean()) / mean_scores.std() # 聚类热图 g = sns.clustermap( mean_scores.T, cmap="RdBu_r", center=0, figsize=figsize, dendrogram_ratio=0.15, cbar_pos=(0.02, 0.8, 0.03, 0.15) ) g.ax_heatmap.set_xlabel("细胞亚群") g.ax_heatmap.set_ylabel("代谢通路") plt.savefig("pathway_heatmap.pdf", bbox_inches="tight") return g逻辑说明：mean_scores先按亚群取平均，得到“亚群 × 通路”矩阵。Z-score 标准化是按通路做的（axis=0方向），让每个通路的打分在亚群间可比，不然高表达通路的绝对值会压过低表达通路。cmap="RdBu_r"是红蓝配色，center=0让零值对应白色，正负差异一目了然。dendrogram_ratio=0.15控制聚类树占图的比例，太大热图区域会被压缩。
参数说明：figsize根据通路数调，80 条通路建议至少 12 英寸高，不然通路名会挤成一团。cbar_pos是 colorbar 的位置，seaborn 的 clustermap 默认 colorbar 位置经常和热图重叠，手动调一下。保存用 PDF 格式，矢量图放大不糊，投稿时直接能用。
4.3 代谢通路的富集方向与生物学解释
打分和差异分析跑完，最后一步是解释。源码里有一个interpret.py，把差异通路的基因集和方向性输出成表格，方便写文章时直接引用。比如某个亚群的糖酵解通路打分显著高，同时氧化磷酸化打分低，这通常提示该亚群偏向糖酵解代谢模式，在肿瘤细胞里是 Warburg 效应的典型表现。源码会把每个差异通路的 leading edge 基因（对打分贡献最大的基因）列出来，这些基因就是后续做实验验证的候选。
注意：代谢通路打分是相对值，不能直接说“这个亚群糖酵解活性是另一个的 2 倍”，只能说“显著高于”。绝对活性需要代谢组学数据来验证，单细胞转录组只能给方向性提示。
5. 避坑与排查：代谢分析里最容易翻车的五个地方
5.1 打分结果全是 NaN 或零
现象：跑完score_pathways，adata.obs里新增的列全是 NaN 或者零。原因通常是基因符号没转换，KEGG 基因集里的人类符号和 adata 里的小鼠符号对不上，valid_genes过滤后一个不剩。解决：在score_pathways之前先跑map_mouse_to_human，把 adata 的var_names替换成人类符号，或者把 KEGG 基因集转成小鼠符号。我一般选前者，因为 KEGG 更新时直接下人类版本就行，不用每次重新映射。
5.2 亚群间代谢差异不显著
现象：差异检验跑完，padj < 0.05的通路只有个位数。原因可能是聚类分辨率太低，亚群分得太粗，代谢异质性被平均掉了。解决：把sc.tl.leiden的resolution从 0.8 提到 1.2 甚至 1.5，重新聚类再打分。另一个原因是ctrl_size设得太大，背景基因集把信号稀释了，调到 30 试试。还有一个容易被忽略的点：如果数据没有做批次校正，批次效应会掩盖真实的代谢差异，先跑sc.external.pp.harmony_integrate再做后续。
5.3 mygene 查询超时或返回空
现象：map_mouse_to_human跑一半报连接超时，或者返回的 mapping 字典是空的。原因：mygene 走的是 NCBI 接口，网络不稳定时会超时。解决：源码里加了重试机制，querymany外面套一个for attempt in range(3)循环，每次失败等 5 秒再试。如果还是不行，可以先把小鼠基因列表存成文件，用 NCBI 的官方datasets工具离线做同源映射，再读进来。返回空的情况通常是基因符号格式不对，检查一下有没有混入 Ensembl ID 或 RefSeq ID，这些要用scopes="ensembl.gene"或scopes="refseq"单独查。
5.4 热图聚类把亚群打乱
现象：plot_pathway_heatmap画出来的热图，亚群顺序和预期不一致，聚类树把不同处理组的亚群混在一起。原因：clustermap默认对行和列都做聚类，列聚类会按代谢打分相似度重排亚群，不按你指定的顺序。解决：如果想让亚群按指定顺序排列，把clustermap的col_cluster设为False，row_cluster保持True让通路聚类。或者用col_linkage传入自定义的聚类树，强制亚群按分组排列。
5.5 内存溢出
现象：跑打分或差异检验时进程被 kill，报MemoryError。原因：单细胞数据细胞数超过 5 万时，adata.obs里存几十个通路的打分列，每个列是 float64，内存占用会到几个 G。解决：打分时把adata.obs的打分列转成 float32，adata.obs[score_name] = adata.obs[score_name].astype(np.float32)，内存直接减半。另外差异检验不要一次性把所有通路的中间结果存 list，用生成器逐通路处理，处理完一个写一个到磁盘。
6. 进阶技巧：把代谢打分嵌进细胞通讯分析
代谢分析做到后面，经常会遇到一个问题：某个亚群的代谢通路活性高，但这个亚群和别的亚群之间有没有代谢物交换？单细胞转录组本身测不到代谢物，但可以通过代谢酶基因的表达来推断。源码里预留了一个cellchat_metab.py，把代谢通路打分和 CellChat 的细胞通讯结果做联合分析。具体做法是：先跑 CellChat 拿到配体-受体对，再筛选出配体或受体是代谢酶的通讯对，看这些通讯对在代谢高活性亚群里是不是更活跃。
import pandas as pd def link_metab_to_cellchat(adata, cellchat_df, metab_pathway="Glycolysis"): """把代谢通路打分和细胞通讯结果关联""" # 取代谢通路打分 score_col = f"metab_{metab_pathway}" if score_col not in adata.obs.columns: raise ValueError(f"{score_col} not found") # 按亚群取平均打分 group_scores = adata.obs.groupby("leiden")[score_col].mean() # 筛选配体或受体是代谢酶的通讯对 metab_genes = set(adata.var_names[adata.var_names.str.contains("Aldo|Eno|Pkm|Ldha")]) cellchat_df["is_metab"] = cellchat_df["ligand"].isin(metab_genes) | \ cellchat_df["receptor"].isin(metab_genes) # 关联：代谢高活性亚群发出的代谢相关通讯对 metab_comms = cellchat_df[cellchat_df["is_metab"]].copy() metab_comms["source_score"] = metab_comms["source"].map(group_scores) metab_comms["target_score"] = metab_comms["target"].map(group_scores) return metab_comms.sort_values("source_score", ascending=False)逻辑说明：metab_genes是手动挑的糖酵解关键酶基因，实际用的时候可以按通路基因集动态取。is_metab标记通讯对里配体或受体是否属于代谢酶。source_score和target_score把代谢打分映射到通讯对的发送方和接收方，这样就能看“代谢活性高的亚群是不是更倾向于发出代谢相关信号”。这个分析在肿瘤微环境研究里很有用，比如看糖酵解高的肿瘤细胞是不是通过乳酸相关信号影响免疫细胞。
参数说明：metab_pathway可以换成任意通路名，但建议选基因数在 20 到 100 之间的通路，太短的通路打分噪声大，太长的通路特异性差。cellchat_df的列名要按实际 CellChat 输出调整，不同版本的 CellChat 列名有差异，跑之前先print(cellchat_df.columns)确认一下。
从那以后我每次跑代谢分析，都会先把基因符号转换那一步单独跑一遍，确认匹配率在 80% 以上再往下走，不然后面全是白费功夫。希望帮到你。
本文还有配套的精品资源，点击获取
