379_单细胞代谢通路活性评分¶

library(Seurat)     # 单细胞分析
library(ggplot2)

# 定义代谢通路基因集（KEGG糖酵解通路）
glycolysis_genes <- list(
  Glycolysis = c("HK1", "HK2", "GPI", "PFKL", "PFKM", "ALDOA",
                 "ALDOC", "GAPDH", "PGK1", "PGAM1", "ENO1",
                 "PKM", "LDHA", "LDHB")   # 糖酵解关键酶基因
)

# 氧化磷酸化通路基因集
oxphos_genes <- list(
  OXPHOS = c("NDUFA1", "NDUFB1", "SDHB", "UQCRB",
             "COX5A", "ATP5A1", "ATP5B", "ATP5D")  # 线粒体呼吸链复合物基因
)

# 计算糖酵解评分
seurat_obj <- AddModuleScore(
  seurat_obj,
  features = glycolysis_genes,    # 目标基因集列表
  name = "Glycolysis_Score",      # 评分列名前缀
  ctrl = 100,                     # 背景基因数（与目标基因集大小相匹配）
  seed = 42                       # 随机种子（保证可重现）
)
# 注意：实际列名会是 "Glycolysis_Score1"（自动加数字后缀）

# 计算氧化磷酸化评分
seurat_obj <- AddModuleScore(
  seurat_obj,
  features = oxphos_genes,
  name = "OXPHOS_Score",
  ctrl = 100,
  seed = 42
)

# 在 UMAP 上可视化代谢评分
p1 <- FeaturePlot(
  seurat_obj,
  features = "Glycolysis_Score1",    # 注意加数字后缀
  reduction = "umap",
  order = TRUE,                       # 高分细胞显示在前景
  pt.size = 0.5
) + scale_color_gradient2(
  low = "blue", mid = "white", high = "red",
  midpoint = 0   # 0为中点（正值=高于背景，负值=低于背景）
) + ggtitle("糖酵解活性")

p2 <- FeaturePlot(
  seurat_obj,
  features = "OXPHOS_Score1",
  reduction = "umap",
  order = TRUE
) + scale_color_gradient2(low = "blue", mid = "white", high = "red",
                           midpoint = 0) + ggtitle("氧化磷酸化活性")

p1 | p2   # 并排显示

# 按细胞类型比较代谢活性
VlnPlot(
  seurat_obj,
  features = c("Glycolysis_Score1", "OXPHOS_Score1"),
  group.by = "cell_type",    # 按细胞类型分组
  pt.size = 0,               # 不显示点（细胞太多时）
  ncol = 2
)

import scanpy as sc
import pandas as pd
import numpy as np
from pyscenic.aucell import aucell
from pyscenic.utils import load_motifs

# 准备基因集（字典格式）
gene_sets = {
    "Glycolysis": ["HK1", "HK2", "GPI", "PFKL", "PFKM", "ALDOA",
                   "ALDOC", "GAPDH", "PGK1", "PGAM1", "ENO1", "PKM"],
    "TCA_Cycle":  ["CS", "ACON", "IDH1", "IDH2", "OGDH", "SUCLA2",
                   "SDH", "FH", "MDH1", "MDH2"],
    "Fatty_Acid_Oxidation": ["CPT1A", "CPT2", "ACADL", "HADH",
                              "HADHA", "HADHB", "ACAA2"]
}

# 加载已预处理的 AnnData
adata = sc.read_h5ad("preprocessed.h5ad")

# 只保留目标基因在数据中存在的
for pathway, genes in gene_sets.items():
    available = [g for g in genes if g in adata.var_names]
    print(f"{pathway}: {len(available)}/{len(genes)} 个基因在数据中")
    gene_sets[pathway] = available   # 更新为可用基因

# 转换为 AUCell 所需的 Regulon 格式
from ctypes import Union
from pyscenic.transform import df2regulons
import collections

# 简单方法：直接用 decoupler 包（更现代的替代品）
import decoupler as dc   # 安装：pip install decoupler

# 从 PROGENy 或自定义基因集获取通路网络
# 自定义通路基因集（转为DataFrame格式）
net_data = []
for pathway, genes in gene_sets.items():
    for gene in genes:
        net_data.append({"source": pathway, "target": gene, "weight": 1})

net_df = pd.DataFrame(net_data)   # 源(通路)-目标(基因)-权重

# 用 ULM（单变量线性模型）评分
dc.run_ulm(
    mat=adata,              # AnnData 对象
    net=net_df,             # 通路-基因网络
    source="source",        # 通路列名
    target="target",        # 基因列名
    weight="weight",        # 权重列名
    use_raw=False
)
# 结果存储在 adata.obsm["ulm_estimate"]

# 提取评分矩阵
scores_df = pd.DataFrame(
    adata.obsm["ulm_estimate"],
    index=adata.obs_names,
    columns=adata.uns["ulm_estimate_names"]
)
print(scores_df.head())

# 将评分添加到 adata.obs 供可视化
for col in scores_df.columns:
    adata.obs[f"pathway_{col}"] = scores_df[col].values

# 可视化
sc.pl.umap(adata, color=["pathway_Glycolysis", "pathway_TCA_Cycle"],
           cmap="RdBu_r",    # 红蓝配色（正=红/活跃，负=蓝/抑制）
           vmin=-2, vmax=2,
           save="_metabolism_scores.pdf")