626_单细胞代谢通路评分¶

# === scMetabolism安装 ===
devtools::install_github("wu-yc/scMetabolism")  # 从GitHub安装

# === 依赖包安装 ===
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c(
    "AUCell",       # AUCell评分方法
    "GSVA",         # GSVA评分方法
    "GSEABase",     # 基因集基础设施
    "clusterProfiler"  # 通路分析
))
install.packages(c("Seurat", "ggplot2", "pheatmap"))  # 其他依赖

library(scMetabolism)    # 加载scMetabolism
library(Seurat)          # 加载Seurat
library(AUCell)          # 加载AUCell
library(GSVA)            # 加载GSVA
library(ggplot2)         # 加载ggplot2

# === 读取Seurat对象 ===
seurat_obj <- readRDS("seurat_processed.rds")  # 读取已处理的Seurat对象

# === 使用scMetabolism计算代谢评分 ===
# 方法1：VISION（默认推荐）
seurat_metabolism <- sc.metabolism.Seurat(
    obj = seurat_obj,                # 输入Seurat对象
    method = "VISION",               # 评分方法（推荐VISION）
    imputation = FALSE,              # 是否做插补（大数据集设FALSE）
    ncores = 8,                      # 线程数
    metabolism.type = "KEGG"         # 代谢通路数据库
)

# 方法2：AUCell
seurat_aucell <- sc.metabolism.Seurat(
    obj = seurat_obj,
    method = "AUCell",               # 用AUCell方法
    imputation = FALSE,
    ncores = 8,
    metabolism.type = "KEGG"
)

# 方法3：ssGSEA
seurat_ssgsea <- sc.metabolism.Seurat(
    obj = seurat_obj,
    method = "ssgsea",               # 用ssGSEA方法
    imputation = FALSE,
    ncores = 8,
    metabolism.type = "KEGG"
)

# metabolism.type选项：
# "KEGG"     — KEGG代谢通路（85条通路）
# "REACTOME" — REACTOME代谢通路（更细致）

# === 提取代谢评分矩阵 ===
metabolism_scores <- seurat_metabolism@assays$METABOLISM$score  # 提取评分

# === 气泡图：按细胞类型展示代谢通路 ===
DotPlot.metabolism(
    obj = seurat_metabolism,           # 带代谢评分的Seurat对象
    pathway = c(
        "Glycolysis / Gluconeogenesis",     # 糖酵解
        "Citrate cycle (TCA cycle)",         # TCA循环
        "Oxidative phosphorylation",         # 氧化磷酸化
        "Fatty acid biosynthesis",           # 脂肪酸合成
        "Glutathione metabolism",            # 谷胱甘肽代谢
        "Purine metabolism"                  # 嘌呤代谢
    ),
    phenotype = "cell_type",           # 按细胞类型分组
    norm = "y"                         # y轴归一化
)

ggsave("metabolism_dotplot.pdf", width = 12, height = 6)

# === 在UMAP上着色显示代谢通路活性 ===
# 将代谢评分存入meta.data
pathway_name <- "Glycolysis / Gluconeogenesis"  # 目标通路
seurat_metabolism$glycolysis_score <- metabolism_scores[pathway_name, ]

# 画UMAP
FeaturePlot(
    seurat_metabolism,
    features = "glycolysis_score",      # 糖酵解评分
    cols = c("lightgrey", "red"),       # 颜色范围
    pt.size = 0.5                       # 点大小
) +
    labs(title = "Glycolysis Activity") +
    theme_minimal()

ggsave("glycolysis_umap.pdf", width = 8, height = 7)

# === 多个通路并排显示 ===
pathways_to_plot <- c(
    "Glycolysis / Gluconeogenesis",
    "Oxidative phosphorylation",
    "Fatty acid biosynthesis"
)

for (pw in pathways_to_plot) {
    safe_name <- gsub("[/ ()]", "_", pw)  # 安全的列名
    seurat_metabolism[[safe_name]] <- metabolism_scores[pw, ]
}

FeaturePlot(
    seurat_metabolism,
    features = gsub("[/ ()]", "_", pathways_to_plot),
    cols = c("lightgrey", "red"),
    ncol = 3,                            # 3列排列
    pt.size = 0.3
)

ggsave("metabolism_umap_multi.pdf", width = 18, height = 6)

# === 手动使用AUCell（更灵活）===
library(AUCell)  # 加载AUCell

# 第1步：构建排名矩阵
# AUCell的核心思想：对每个细胞的基因按表达量排名
# 然后看目标基因集中的基因是否排在前面
expr_matrix <- GetAssayData(seurat_obj, slot = "counts")  # 获取表达矩阵
cells_rankings <- AUCell_buildRankings(
    expr_matrix,                     # 表达矩阵
    nCores = 8,                      # 线程数
    plotStats = TRUE                 # 画排名统计图
)

# 第2步：准备基因集
# 从KEGG获取代谢通路基因集
library(clusterProfiler)  # 加载clusterProfiler
library(org.Hs.eg.db)     # 人类基因注释

# 获取KEGG通路的基因列表
kegg_pathways <- download_KEGG("hsa")  # 下载人类KEGG通路
pathway_genes <- kegg_pathways$KEGGPATHID2EXTID  # 通路-基因映射

# 转换基因ID格式（KEGG用Entrez ID，需转为Symbol）
# 构建GeneSet列表
gene_sets <- list()
for (pw_id in unique(pathway_genes$from)) {
    entrez_ids <- pathway_genes$to[pathway_genes$from == pw_id]
    symbols <- mapIds(org.Hs.eg.db, 
                      keys = entrez_ids, 
                      column = "SYMBOL",
                      keytype = "ENTREZID",
                      multiVals = "first")
    symbols <- symbols[!is.na(symbols)]  # 移除NA
    if (length(symbols) >= 5) {  # 至少5个基因
        gene_sets[[pw_id]] <- symbols
    }
}

# 第3步：计算AUC分数
cells_AUC <- AUCell_calcAUC(
    gene_sets,                       # 基因集列表
    cells_rankings,                  # 排名矩阵
    nCores = 8                       # 线程数
)

# 第4步：提取评分矩阵
auc_matrix <- getAUC(cells_AUC)      # 提取AUC矩阵
dim(auc_matrix)                       # 通路数 x 细胞数

# === 使用GSVA方法 ===
library(GSVA)  # 加载GSVA

# 准备表达矩阵（用归一化后的数据）
expr_for_gsva <- as.matrix(
    GetAssayData(seurat_obj, slot = "data")  # 归一化数据
)

# 运行GSVA
gsva_scores <- gsva(
    expr_for_gsva,                   # 表达矩阵
    gene_sets,                       # 基因集列表
    method = "gsva",                 # GSVA方法
    kcdf = "Poisson",                # 适合计数数据的核密度
    parallel.sz = 8,                 # 8线程
    verbose = TRUE                   # 显示进度
)

# GSVA vs AUCell比较：
# GSVA：对整体分布更敏感，适合样本间比较
# AUCell：对稀疏数据更鲁棒，适合单细胞
# ssGSEA：每个细胞独立计算，中间方案

# === 比较不同细胞类型的代谢活性 ===
# 将AUCell评分添加到Seurat对象
for (pw in rownames(auc_matrix)) {
    safe_name <- gsub("[^a-zA-Z0-9]", "_", pw)  # 安全列名
    seurat_obj[[safe_name]] <- auc_matrix[pw, colnames(seurat_obj)]
}

# 小提琴图比较代谢通路
VlnPlot(
    seurat_obj,
    features = c("hsa00010", "hsa00020", "hsa00190"),  # 糖酵解/TCA/氧磷
    group.by = "cell_type",            # 按细胞类型分组
    pt.size = 0,                       # 不显示散点
    ncol = 3                           # 3列
) & 
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

# === 疾病vs对照的代谢差异 ===
# 使用Wilcoxon检验比较每个通路的评分
library(dplyr)  # 数据处理
library(tidyr)  # 数据整理

# 提取元数据和代谢评分
meta <- seurat_obj@meta.data  # 元数据
auc_df <- as.data.frame(t(auc_matrix))  # 转置AUC矩阵
auc_df$condition <- meta$condition  # 添加条件列
auc_df$cell_type <- meta$cell_type  # 添加细胞类型列

# 对每个通路做统计检验
results <- data.frame()  # 存储结果
for (pw in rownames(auc_matrix)) {
    for (ct in unique(meta$cell_type)) {
        subset_data <- auc_df[auc_df$cell_type == ct, ]  # 子集
        if (nrow(subset_data) > 10) {  # 至少10个细胞
            test <- wilcox.test(
                subset_data[[pw]][subset_data$condition == "disease"],
                subset_data[[pw]][subset_data$condition == "control"]
            )
            results <- rbind(results, data.frame(
                pathway = pw,
                cell_type = ct,
                p_value = test$p.value,
                mean_disease = mean(subset_data[[pw]][subset_data$condition == "disease"]),
                mean_control = mean(subset_data[[pw]][subset_data$condition == "control"])
            ))
        }
    }
}

# FDR校正
results$fdr <- p.adjust(results$p_value, method = "BH")  # BH校正
results$log2fc <- log2(results$mean_disease / results$mean_control)  # 效应量

# 显著的代谢差异
sig_results <- results[results$fdr < 0.05, ]  # FDR < 0.05
sig_results <- sig_results[order(sig_results$fdr), ]  # 按FDR排序
head(sig_results, 20)  # 前20个

# === 热图可视化 ===
library(pheatmap)  # 加载pheatmap

# 准备热图数据
heatmap_data <- results %>%
    dplyr::filter(fdr < 0.05) %>%               # 只保留显著的
    dplyr::select(pathway, cell_type, log2fc) %>%
    tidyr::pivot_wider(names_from = cell_type,   # 宽格式
                       values_from = log2fc,
                       values_fill = 0) %>%
    tibble::column_to_rownames("pathway")

pheatmap(
    as.matrix(heatmap_data),
    color = colorRampPalette(c("blue", "white", "red"))(100),
    cluster_rows = TRUE,             # 行聚类
    cluster_cols = TRUE,             # 列聚类
    fontsize_row = 8,                # 行标签大小
    main = "Metabolic Pathway Changes (Disease vs Control)"
)