单细胞代谢分析¶

一句话概述¶

利用scFEA、Compass和scFBA等计算工具，从单细胞RNA-seq数据推断单个细胞的代谢通量活性，揭示不同细胞类型/状态间的代谢异质性，是研究肿瘤代谢重编程、免疫代谢和细胞命运决定的前沿分析方法。

核心知识点表格¶

知识点	关键内容	常用工具
代谢通量推断	从基因表达推算代谢反应活性	scFEA, Compass, scFBA
代谢网络模型	基因组尺度代谢模型(GEM)	Recon3D, Human1
通量平衡分析	稳态假设下最优化求解	COBRApy, MATLAB COBRA
代谢通路评分	ssGSEA/AUCell代谢基因集	AUCell, GSVA, scMetabolism
代谢异质性	细胞间代谢差异与功能关联	-
代谢-转录耦合	代谢酶表达与代谢物水平关系	-
实验验证	单细胞代谢组、Seahorse	-

各步骤详解¶

第一步：理解单细胞代谢分析原理¶

白话解释： 细胞的代谢状态（用什么营养、产什么能量、合成什么物质）对其功能至关重要。但直接在单细胞水平测代谢物很困难。因此，人们开发了计算方法，利用单细胞RNA-seq数据中代谢酶基因的表达来推断代谢反应的活性——就像通过工厂中工人的数量推断生产线的产能。

技术细节： - 基因组尺度代谢模型(GEM)：包含所有已知代谢反应及其对应酶的数学模型 - 通量平衡分析(FBA)：假设代谢处于稳态，用线性规划求解代谢通量分布 - scRNA-seq → 代谢的挑战：转录本水平≠蛋白水平≠酶活性≠代谢通量 - 基本假设：酶基因高表达→该反应更活跃（合理但不完美的近似）

第二步：scFEA代谢通量推断¶

白话解释： scFEA (single-cell Flux Estimation Analysis) 使用图神经网络学习基因表达到代谢通量的映射关系。它把代谢网络建模为一个图，节点是代谢物，边是反应，然后用深度学习预测每个细胞中每条反应的通量。

代码示例：

# 安装scFEA
# pip install scFEA
# 或从GitHub: https://github.com/changwn/scFEA

import scFEA
import scanpy as sc
import pandas as pd
import numpy as np

# 1. 准备scRNA-seq数据
adata = sc.read_h5ad("scRNA_processed.h5ad")

# 标准化（scFEA需要log-normalized数据）
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# 2. 导出表达矩阵（cells × genes）
expr_matrix = pd.DataFrame(
    adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X,
    index=adata.obs_names,
    columns=adata.var_names)
expr_matrix.to_csv("expression_matrix.csv")

# 3. 运行scFEA (命令行)
# python scFEA.py --data_dir ./ \
#     --input_dir expression_matrix.csv \
#     --res_dir scfea_output/ \
#     --moduleGene_file module_gene_m168.csv \
#     --stoichiometry_matrix cmMat_c70_m168.csv \
#     --sc_imputation True \
#     --output_flux_file flux.csv \
#     --output_balance_file balance.csv

# 4. 或使用Python API
from scFEA.src import scFEA_main

results = scFEA_main.run_scFEA(
    expression_file="expression_matrix.csv",
    module_gene_file="module_gene_m168.csv",
    stoichiometry_file="cmMat_c70_m168.csv",
    output_dir="scfea_output/",
    epochs=100,
    batch_size=64)

# 5. 读取结果
flux_matrix = pd.read_csv("scfea_output/flux.csv", index_col=0)
# flux_matrix: cells × metabolic_modules

# 6. 将通量添加到adata
adata.obsm["X_flux"] = flux_matrix.values

# 7. 可视化代谢通量差异
import matplotlib.pyplot as plt

# 在UMAP上展示特定代谢通量
sc.pl.umap(adata, color=["cell_type"])

# 比较不同细胞类型的糖酵解通量
glycolysis_modules = ["M_1", "M_2", "M_3"]  # 糖酵解相关模块
for mod in glycolysis_modules:
    adata.obs[mod] = flux_matrix[mod].values
    sc.pl.umap(adata, color=mod, title=f"Glycolysis module: {mod}")

第三步：Compass代谢分析¶

白话解释： Compass是另一种单细胞代谢分析工具，基于FBA方法。它为每个细胞构建一个约束的代谢模型，用基因表达数据作为反应上界约束，然后用优化方法求解每条反应的最大可能通量。输出的"Compass score"反映每个细胞利用该代谢反应的能力。

代码示例：

# 安装Compass（不在PyPI上，需从GitHub安装，且依赖Gurobi求解器）
# python -m pip install git+https://github.com/wagnerlab-berkeley/Compass.git --upgrade
# 注意：Compass需要Gurobi许可证（学术用户免费申请WLS许可）

import scanpy as sc
import pandas as pd

# 1. 准备数据（需要TPM或CPM标准化，不做log变换）
adata = sc.read_h5ad("scRNA_processed.h5ad")
sc.pp.normalize_total(adata, target_sum=1e6)  # CPM
# 注意：Compass官方文档建议输入normalized counts，不做log1p

# 导出
expr_df = pd.DataFrame(
    adata.X.toarray(),
    index=adata.obs_names,
    columns=adata.var_names).T  # Compass需要 genes × cells

expr_df.to_csv("compass_input.tsv", sep="\t")

# 2. 运行Compass (命令行工具，不是Python库导入)
# compass --data compass_input.tsv \
#     --num-processes 16 \
#     --output-dir compass_output/ \
#     --species homo_sapiens \
#     --num-neighbors 30

# 3. 解析结果
reactions = pd.read_csv("compass_output/reactions.tsv", sep="\t", index_col=0)
# reactions: reactions × cells (Compass scores)
# 分数越高，该细胞利用该反应的能力越强

# 4. 代谢通路水平汇总
# 将反应映射到KEGG通路
pathway_mapping = pd.read_csv("reaction_to_pathway.csv")

# 计算通路水平得分
pathway_scores = {}
for pathway, rxns in pathway_mapping.groupby("pathway"):
    rxn_ids = rxns["reaction_id"].tolist()
    available_rxns = [r for r in rxn_ids if r in reactions.index]
    if available_rxns:
        pathway_scores[pathway] = reactions.loc[available_rxns].mean(axis=0)

pathway_df = pd.DataFrame(pathway_scores)

# 5. 差异代谢分析（比较细胞类型）
from scipy.stats import mannwhitneyu

cell_types = adata.obs["cell_type"]
results = []
for pathway in pathway_df.columns:
    for ct in cell_types.unique():
        mask = cell_types == ct
        score_in = pathway_df.loc[mask, pathway]
        score_out = pathway_df.loc[~mask, pathway]
        stat, pval = mannwhitneyu(score_in, score_out, alternative="two-sided")
        results.append({"pathway": pathway, "cell_type": ct, 
                       "median_score": score_in.median(), "pvalue": pval})

results_df = pd.DataFrame(results)
results_df["padj"] = results_df.groupby("pathway")["pvalue"].transform(
    lambda x: x * len(x) / (x.rank()))  # BH correction

第四步：scMetabolism代谢通路评分¶

白话解释： 如果不需要精确的通量计算，可以用更简单的方法——对代谢通路的基因集做评分（类似GSVA）。scMetabolism就是这样的工具，它直接计算每个细胞在各代谢通路上的活性得分。

代码示例：

# 安装scMetabolism
devtools::install_github("wu-yc/scMetabolism")
library(scMetabolism)
library(Seurat)

# 加载Seurat对象
seurat_obj <- readRDS("seurat_processed.rds")

# 运行scMetabolism
# 方法: VISION, AUCell, ssGSEA, GSVA
metabolism_result <- sc.metabolism.Seurat(
    obj = seurat_obj,
    method = "AUCell",           # 推荐AUCell
    imputation = FALSE,
    ncores = 8,
    metabolism.type = "KEGG")    # KEGG或Reactome

# 结果存储在seurat对象的assay中
# metabolism_result@assays$METABOLISM$score

# 可视化
DotPlot.metabolism(obj = metabolism_result,
    pathway.number = 20,        # Top 20通路
    group.by = "cell_type")

# 热图
library(pheatmap)
pathway_scores <- GetAssayData(metabolism_result, assay="METABOLISM")
avg_scores <- AverageExpression(metabolism_result, 
    assays="METABOLISM", group.by="cell_type")
pheatmap(avg_scores$METABOLISM[1:30, ], scale="row",
    cluster_cols=TRUE, fontsize_row=8)

# 比较特定通路在不同细胞类型中的活性
VlnPlot(metabolism_result, 
    features = c("Glycolysis / Gluconeogenesis",
                 "Oxidative phosphorylation",
                 "Fatty acid oxidation"),
    group.by = "cell_type",
    assay = "METABOLISM")

第五步：代谢异质性分析与功能关联¶

白话解释： 得到代谢评分后，可以分析：1）同一细胞类型内部的代谢异质性（如肿瘤细胞间的代谢差异）；2）代谢状态与细胞功能状态的关系（如T细胞的糖酵解活性与激活状态的关系）。

代码示例：

import scanpy as sc
import numpy as np
from scipy.stats import spearmanr

# 1. 代谢轨迹分析
# 在T细胞中分析代谢与激活状态的关系
t_cells = adata[adata.obs["cell_type"] == "CD8_T"]

# 计算激活评分
activation_genes = ["GZMB", "PRF1", "IFNG", "TNF", "FASLG"]
sc.tl.score_genes(t_cells, gene_list=activation_genes, 
    score_name="activation_score")

# 代谢-功能相关性
glycolysis_score = t_cells.obs["Glycolysis_score"]
activation_score = t_cells.obs["activation_score"]
r, p = spearmanr(glycolysis_score, activation_score)
print(f"糖酵解 vs 激活: r={r:.3f}, p={p:.2e}")

# 2. 代谢亚群发现
# 基于代谢评分进行聚类
from sklearn.cluster import KMeans

metabolic_features = pathway_df.loc[t_cells.obs_names]
kmeans = KMeans(n_clusters=3, random_state=42)
t_cells.obs["metabolic_cluster"] = kmeans.fit_predict(metabolic_features).astype(str)

# 3. 差异代谢通路（代谢亚群vs功能表型）
sc.tl.rank_genes_groups(t_cells, groupby="metabolic_cluster")
sc.pl.rank_genes_groups_dotplot(t_cells, n_genes=10)

# 4. 代谢通量与临床结局
# 如有配对的临床信息（如治疗反应）
responders = adata[adata.obs["response"] == "responder"]
non_responders = adata[adata.obs["response"] == "non_responder"]

# 比较免疫治疗应答者vs非应答者的T细胞代谢

第六步：整合分析与验证策略¶

代码示例：

# 多方法一致性验证
# 比较scFEA、Compass、scMetabolism的结果一致性

from scipy.stats import spearmanr
import seaborn as sns

# 糖酵解通量的方法间相关性
methods_comparison = pd.DataFrame({
    "scFEA": flux_matrix["glycolysis_module"].values,
    "Compass": compass_scores["glycolysis"].values,
    "AUCell": aucell_scores["Glycolysis"].values
}, index=adata.obs_names)

# 相关性热图
corr_matrix = methods_comparison.corr(method="spearman")
sns.heatmap(corr_matrix, annot=True, cmap="RdBu_r", center=0)
plt.title("Method Agreement for Glycolysis")
plt.savefig("method_comparison.pdf")

实战命令¶

#!/bin/bash
# === 单细胞代谢分析流程 ===

# 1. scFEA通量推断
python scFEA.py --data_dir ./ --input_dir expr_matrix.csv \
    --res_dir scfea_results/ --moduleGene_file module_gene_m168.csv \
    --stoichiometry_matrix cmMat_c70_m168.csv --sc_imputation True

# 2. Compass分析（需先安装Gurobi并配置许可证）
compass --data compass_input.tsv --num-processes 16 \
    --output-dir compass_results/ --species homo_sapiens \
    --num-neighbors 30

# 3. R中运行scMetabolism + 可视化
Rscript run_scMetabolism.R

面试常问点¶

Q1: 从scRNA-seq推断代谢通量的基本假设是什么？其局限性？ A: 核心假设是酶基因表达水平与催化反应的通量正相关。局限：1）转录≠翻译≠酶活性（翻译后修饰、蛋白稳定性）；2）代谢受底物浓度和别构调节影响，非纯转录控制；3）scRNA-seq数据稀疏，零值多；4）细胞外营养物质条件未知。是一种近似而非精确方法。

Q2: scFEA和Compass的区别？ A: scFEA使用图神经网络学习表达→通量的非线性映射，输出绝对通量估计。Compass基于经典FBA框架，用表达做约束条件求解最大可能通量，输出相对代谢能力评分。scFEA更灵活但需要训练，Compass理论基础更明确但计算量大。

Q3: 如何验证计算预测的代谢状态？ A: 1）与bulk代谢组学数据比较（如同一样本的LC-MS数据）；2）Seahorse实验测量OCR/ECAR验证氧化磷酸化/糖酵解比例；3）使用代谢物reporter基因（如乳酸传感器）验证；4）单细胞代谢组学（如SpaceM）直接验证；5）药物抑制实验验证预测的代谢依赖性。

Q4: 为什么要在单细胞水平研究代谢？ A: 同一组织中不同细胞的代谢状态可能完全不同。如：肿瘤中心缺氧区域的细胞依赖糖酵解，外周细胞用氧化磷酸化；激活态T细胞转向有氧糖酵解（Warburg效应），而记忆T细胞依赖脂肪酸氧化。Bulk分析只能看到平均值，无法揭示这些异质性。

易错点¶

过度信任计算预测：代谢通量推断是近似方法，不能等同于实验测量。重要结论需要实验验证。
零值问题未处理：scRNA-seq的dropout导致大量酶基因被检测为零，但不意味着该代谢反应不活跃。需要使用imputation或考虑零膨胀模型。
混淆代谢通路评分与通量：ssGSEA/AUCell给出的是基因集活性"评分"（无单位），不是真正的代谢通量（mmol/gDW/h）。不可直接比较不同通路的"绝对值"。
忽略细胞外环境：代谢不仅取决于酶表达，还取决于营养物质供应。体外培养和体内环境代谢可能完全不同。
过度简化代谢网络：使用简化的通路基因集会遗漏重要的辅路和旁路反应。

补充知识¶