简介这份源码资源面向从事单细胞转录组与代谢研究的生信分析人员及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 python3.9 -y conda activate sc_metab # scanpy 锁在 1.9.x1.10 之后 score_genes 的默认参数变了 pip install scanpy1.9.6 pip install mygene3.2.2 pip install scipy1.10.1 pip install pandas1.5.3 pip install matplotlib3.7.1 pip install seaborn0.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_namesgene_symbols) # 基础质控线粒体基因比例过滤 adata.var[mt] adata.var_names.str.startswith(mt-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], inplaceTrue) 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_sum1e4) sc.pp.log1p(adata) # 高变基因与降维代谢分析前必须做否则聚类不准 sc.pp.highly_variable_genes(adata, n_top_genes2000, flavorseurat) adata adata[:, adata.var[highly_variable]].copy() sc.pp.scale(adata, max_value10) sc.tl.pca(adata, n_comps30) sc.pp.neighbors(adata, n_neighbors15, n_pcs30) sc.tl.leiden(adata, resolution0.8)参数说明pct_counts_mt 20是小鼠数据的经验阈值人类数据一般卡 10小鼠线粒体基因少、比例天然偏高卡太严会丢掉真实细胞。n_top_genes2000是 scanpy 教程的默认值但代谢分析里我建议提到 3000因为代谢基因本身表达低2000 个高变基因里可能只覆盖到几十个代谢基因。resolution0.8是聚类分辨率做代谢异质性比较时如果亚群分得太粗代谢差异会被平均掉我一般会试 0.6 到 1.2 三个值看哪个分辨率下代谢通路打分在亚群间有显著差异。注意如果你的数据来自 Smart-seq2 而不是 10xnormalize_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, scopessymbol, speciesmouse, fieldshomologene, symbol, as_dataframeTrue ) # 取 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的scopessymbol表示输入是小鼠基因符号speciesmouse限定物种fieldshomologene取同源基因信息。返回的 DataFrame 里homologene字段是一个嵌套字典人类同源符号在human.symbol路径下。这个查询走的是 NCBI 的接口一次查几千个基因大概要 2 到 3 分钟源码里加了缓存机制第二次跑同一批基因直接读本地 pickle不用重复请求。参数说明as_dataframeTrue让返回结果直接是 DataFrame方便后续遍历。如果查询的基因里有别名或旧符号mygene 会自动做模糊匹配但会在notfound字段里标记出来源码里对notfound的基因做了二次查询用scopesalias再试一次能再捞回来 5% 左右。3.2 KEGG 代谢通路基因集的加载与过滤源码的pathways/目录下有一个kegg_metabolism.json里面是手动整理过的 KEGG 代谢通路基因集去掉了疾病通路和信号通路只保留碳水化合物代谢、脂质代谢、氨基酸代谢、能量代谢这几大类。加载和过滤的代码import json def load_metabolic_pathways(pathpathways/kegg_metabolism.json, min_genes10): 加载代谢通路基因集过滤掉基因数太少的通路 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_genes10是个经验值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_size50): 对每个代谢通路打分结果写入 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 fmetab_{name} sc.tl.score_genes( adata, gene_listvalid_genes, ctrl_sizectrl_size, score_namescore_name, random_state42 ) return adata逻辑说明valid_genes过滤掉数据里不存在的基因如果通路里有效基因少于 5 个就跳过避免噪声。ctrl_size50是背景基因集的基因数源码默认 50我试过 100打分结果更平滑但亚群间差异也被抹平了50 是个平衡点。random_state42固定随机种子保证每次跑结果一致这个在做多次分析对比时很重要。参数说明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_keyleiden, score_prefixmetab_): 比较各亚群间代谢通路打分的差异 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[i1:]: 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, alternativetwo-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], methodfdr_bh)[1] return df[df[padj] 0.05].sort_values(padj)逻辑说明Mann-Whitney U 检验不假设正态分布适合单细胞打分这种偏态数据。alternativetwo-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_keyleiden, figsize(12, 8)): 画亚群-代谢通路热图 # 按亚群取每个通路的平均打分 score_cols [fmetab_{p} for p in pathways if fmetab_{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, cmapRdBu_r, center0, figsizefigsize, dendrogram_ratio0.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_inchestight) return g逻辑说明mean_scores先按亚群取平均得到“亚群 × 通路”矩阵。Z-score 标准化是按通路做的axis0方向让每个通路的打分在亚群间可比不然高表达通路的绝对值会压过低表达通路。cmapRdBu_r是红蓝配色center0让零值对应白色正负差异一目了然。dendrogram_ratio0.15控制聚类树占图的比例太大热图区域会被压缩。参数说明figsize根据通路数调80 条通路建议至少 12 英寸高不然通路名会挤成一团。cbar_pos是 colorbar 的位置seaborn 的 clustermap 默认 colorbar 位置经常和热图重叠手动调一下。保存用 PDF 格式矢量图放大不糊投稿时直接能用。4.3 代谢通路的富集方向与生物学解释打分和差异分析跑完最后一步是解释。源码里有一个interpret.py把差异通路的基因集和方向性输出成表格方便写文章时直接引用。比如某个亚群的糖酵解通路打分显著高同时氧化磷酸化打分低这通常提示该亚群偏向糖酵解代谢模式在肿瘤细胞里是 Warburg 效应的典型表现。源码会把每个差异通路的 leading edge 基因对打分贡献最大的基因列出来这些基因就是后续做实验验证的候选。注意代谢通路打分是相对值不能直接说“这个亚群糖酵解活性是另一个的 2 倍”只能说“显著高于”。绝对活性需要代谢组学数据来验证单细胞转录组只能给方向性提示。5. 避坑与排查代谢分析里最容易翻车的五个地方5.1 打分结果全是 NaN 或零现象跑完score_pathwaysadata.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这些要用scopesensembl.gene或scopesrefseq单独查。5.4 热图聚类把亚群打乱现象plot_pathway_heatmap画出来的热图亚群顺序和预期不一致聚类树把不同处理组的亚群混在一起。原因clustermap默认对行和列都做聚类列聚类会按代谢打分相似度重排亚群不按你指定的顺序。解决如果想让亚群按指定顺序排列把clustermap的col_cluster设为Falserow_cluster保持True让通路聚类。或者用col_linkage传入自定义的聚类树强制亚群按分组排列。5.5 内存溢出现象跑打分或差异检验时进程被 kill报MemoryError。原因单细胞数据细胞数超过 5 万时adata.obs里存几十个通路的打分列每个列是 float64内存占用会到几个 G。解决打分时把adata.obs的打分列转成 float32adata.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_pathwayGlycolysis): 把代谢通路打分和细胞通讯结果关联 # 取代谢通路打分 score_col fmetab_{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, ascendingFalse)逻辑说明metab_genes是手动挑的糖酵解关键酶基因实际用的时候可以按通路基因集动态取。is_metab标记通讯对里配体或受体是否属于代谢酶。source_score和target_score把代谢打分映射到通讯对的发送方和接收方这样就能看“代谢活性高的亚群是不是更倾向于发出代谢相关信号”。这个分析在肿瘤微环境研究里很有用比如看糖酵解高的肿瘤细胞是不是通过乳酸相关信号影响免疫细胞。参数说明metab_pathway可以换成任意通路名但建议选基因数在 20 到 100 之间的通路太短的通路打分噪声大太长的通路特异性差。cellchat_df的列名要按实际 CellChat 输出调整不同版本的 CellChat 列名有差异跑之前先print(cellchat_df.columns)确认一下。从那以后我每次跑代谢分析都会先把基因符号转换那一步单独跑一遍确认匹配率在 80% 以上再往下走不然后面全是白费功夫。希望帮到你。本文还有配套的精品资源点击获取
简介这份源码资源面向从事单细胞转录组与代谢研究的生信分析人员及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 python3.9 -y conda activate sc_metab # scanpy 锁在 1.9.x1.10 之后 score_genes 的默认参数变了 pip install scanpy1.9.6 pip install mygene3.2.2 pip install scipy1.10.1 pip install pandas1.5.3 pip install matplotlib3.7.1 pip install seaborn0.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_namesgene_symbols) # 基础质控线粒体基因比例过滤 adata.var[mt] adata.var_names.str.startswith(mt-) sc.pp.calculate_qc_metrics(adata, qc_vars[mt], inplaceTrue) 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_sum1e4) sc.pp.log1p(adata) # 高变基因与降维代谢分析前必须做否则聚类不准 sc.pp.highly_variable_genes(adata, n_top_genes2000, flavorseurat) adata adata[:, adata.var[highly_variable]].copy() sc.pp.scale(adata, max_value10) sc.tl.pca(adata, n_comps30) sc.pp.neighbors(adata, n_neighbors15, n_pcs30) sc.tl.leiden(adata, resolution0.8)参数说明pct_counts_mt 20是小鼠数据的经验阈值人类数据一般卡 10小鼠线粒体基因少、比例天然偏高卡太严会丢掉真实细胞。n_top_genes2000是 scanpy 教程的默认值但代谢分析里我建议提到 3000因为代谢基因本身表达低2000 个高变基因里可能只覆盖到几十个代谢基因。resolution0.8是聚类分辨率做代谢异质性比较时如果亚群分得太粗代谢差异会被平均掉我一般会试 0.6 到 1.2 三个值看哪个分辨率下代谢通路打分在亚群间有显著差异。注意如果你的数据来自 Smart-seq2 而不是 10xnormalize_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, scopessymbol, speciesmouse, fieldshomologene, symbol, as_dataframeTrue ) # 取 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的scopessymbol表示输入是小鼠基因符号speciesmouse限定物种fieldshomologene取同源基因信息。返回的 DataFrame 里homologene字段是一个嵌套字典人类同源符号在human.symbol路径下。这个查询走的是 NCBI 的接口一次查几千个基因大概要 2 到 3 分钟源码里加了缓存机制第二次跑同一批基因直接读本地 pickle不用重复请求。参数说明as_dataframeTrue让返回结果直接是 DataFrame方便后续遍历。如果查询的基因里有别名或旧符号mygene 会自动做模糊匹配但会在notfound字段里标记出来源码里对notfound的基因做了二次查询用scopesalias再试一次能再捞回来 5% 左右。3.2 KEGG 代谢通路基因集的加载与过滤源码的pathways/目录下有一个kegg_metabolism.json里面是手动整理过的 KEGG 代谢通路基因集去掉了疾病通路和信号通路只保留碳水化合物代谢、脂质代谢、氨基酸代谢、能量代谢这几大类。加载和过滤的代码import json def load_metabolic_pathways(pathpathways/kegg_metabolism.json, min_genes10): 加载代谢通路基因集过滤掉基因数太少的通路 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_genes10是个经验值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_size50): 对每个代谢通路打分结果写入 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 fmetab_{name} sc.tl.score_genes( adata, gene_listvalid_genes, ctrl_sizectrl_size, score_namescore_name, random_state42 ) return adata逻辑说明valid_genes过滤掉数据里不存在的基因如果通路里有效基因少于 5 个就跳过避免噪声。ctrl_size50是背景基因集的基因数源码默认 50我试过 100打分结果更平滑但亚群间差异也被抹平了50 是个平衡点。random_state42固定随机种子保证每次跑结果一致这个在做多次分析对比时很重要。参数说明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_keyleiden, score_prefixmetab_): 比较各亚群间代谢通路打分的差异 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[i1:]: 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, alternativetwo-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], methodfdr_bh)[1] return df[df[padj] 0.05].sort_values(padj)逻辑说明Mann-Whitney U 检验不假设正态分布适合单细胞打分这种偏态数据。alternativetwo-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_keyleiden, figsize(12, 8)): 画亚群-代谢通路热图 # 按亚群取每个通路的平均打分 score_cols [fmetab_{p} for p in pathways if fmetab_{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, cmapRdBu_r, center0, figsizefigsize, dendrogram_ratio0.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_inchestight) return g逻辑说明mean_scores先按亚群取平均得到“亚群 × 通路”矩阵。Z-score 标准化是按通路做的axis0方向让每个通路的打分在亚群间可比不然高表达通路的绝对值会压过低表达通路。cmapRdBu_r是红蓝配色center0让零值对应白色正负差异一目了然。dendrogram_ratio0.15控制聚类树占图的比例太大热图区域会被压缩。参数说明figsize根据通路数调80 条通路建议至少 12 英寸高不然通路名会挤成一团。cbar_pos是 colorbar 的位置seaborn 的 clustermap 默认 colorbar 位置经常和热图重叠手动调一下。保存用 PDF 格式矢量图放大不糊投稿时直接能用。4.3 代谢通路的富集方向与生物学解释打分和差异分析跑完最后一步是解释。源码里有一个interpret.py把差异通路的基因集和方向性输出成表格方便写文章时直接引用。比如某个亚群的糖酵解通路打分显著高同时氧化磷酸化打分低这通常提示该亚群偏向糖酵解代谢模式在肿瘤细胞里是 Warburg 效应的典型表现。源码会把每个差异通路的 leading edge 基因对打分贡献最大的基因列出来这些基因就是后续做实验验证的候选。注意代谢通路打分是相对值不能直接说“这个亚群糖酵解活性是另一个的 2 倍”只能说“显著高于”。绝对活性需要代谢组学数据来验证单细胞转录组只能给方向性提示。5. 避坑与排查代谢分析里最容易翻车的五个地方5.1 打分结果全是 NaN 或零现象跑完score_pathwaysadata.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这些要用scopesensembl.gene或scopesrefseq单独查。5.4 热图聚类把亚群打乱现象plot_pathway_heatmap画出来的热图亚群顺序和预期不一致聚类树把不同处理组的亚群混在一起。原因clustermap默认对行和列都做聚类列聚类会按代谢打分相似度重排亚群不按你指定的顺序。解决如果想让亚群按指定顺序排列把clustermap的col_cluster设为Falserow_cluster保持True让通路聚类。或者用col_linkage传入自定义的聚类树强制亚群按分组排列。5.5 内存溢出现象跑打分或差异检验时进程被 kill报MemoryError。原因单细胞数据细胞数超过 5 万时adata.obs里存几十个通路的打分列每个列是 float64内存占用会到几个 G。解决打分时把adata.obs的打分列转成 float32adata.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_pathwayGlycolysis): 把代谢通路打分和细胞通讯结果关联 # 取代谢通路打分 score_col fmetab_{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, ascendingFalse)逻辑说明metab_genes是手动挑的糖酵解关键酶基因实际用的时候可以按通路基因集动态取。is_metab标记通讯对里配体或受体是否属于代谢酶。source_score和target_score把代谢打分映射到通讯对的发送方和接收方这样就能看“代谢活性高的亚群是不是更倾向于发出代谢相关信号”。这个分析在肿瘤微环境研究里很有用比如看糖酵解高的肿瘤细胞是不是通过乳酸相关信号影响免疫细胞。参数说明metab_pathway可以换成任意通路名但建议选基因数在 20 到 100 之间的通路太短的通路打分噪声大太长的通路特异性差。cellchat_df的列名要按实际 CellChat 输出调整不同版本的 CellChat 列名有差异跑之前先print(cellchat_df.columns)确认一下。从那以后我每次跑代谢分析都会先把基因符号转换那一步单独跑一遍确认匹配率在 80% 以上再往下走不然后面全是白费功夫。希望帮到你。本文还有配套的精品资源点击获取
