380_单细胞基因集评分方法比较¶

场景A：快速探索 → AddModuleScore（Seurat内置，最省事）
场景B：稀疏性问题严重 → AUCell 或 UCell
场景C：发表级别严格分析 → UCell + 多方法交叉验证
场景D：有完整调控网络数据 → decoupler（ULM/MLM）
场景E：空间转录组联合 → VISION

library(Seurat)      # 单细胞分析
library(UCell)       # 单细胞优化的基因集评分
library(AUCell)      # AUC评分
library(ggplot2)
library(dplyr)

# ============================================================
# 准备测试基因集（肿瘤相关，典型应用场景）
# ============================================================
gene_sets <- list(
  Hypoxia     = c("HIF1A", "VEGFA", "SLC2A1", "LDHA", "PGK1", "PGAM1", "ENO1"),
  Proliferation = c("MKI67", "TOP2A", "PCNA", "MCM2", "CDK4", "CCND1"),
  EMT         = c("VIM", "CDH2", "SNAI1", "SNAI2", "ZEB1", "FN1", "COL1A1")
)

# ============================================================
# 方法1：AddModuleScore
# ============================================================
seurat_obj <- AddModuleScore(seurat_obj, features = gene_sets["Hypoxia"],
                              name = "AMS_Hypoxia", seed = 42)
seurat_obj <- AddModuleScore(seurat_obj, features = gene_sets["Proliferation"],
                              name = "AMS_Prolif", seed = 42)
seurat_obj <- AddModuleScore(seurat_obj, features = gene_sets["EMT"],
                              name = "AMS_EMT", seed = 42)

# ============================================================
# 方法2：UCell（推荐替代）
# ============================================================
seurat_obj <- AddModuleScore_UCell(
  seurat_obj,
  features = gene_sets,     # 一次性输入所有基因集
  maxRank = 1500,           # 参与排名的最大基因数（影响敏感度）
  ncores = 4
)
# 结果存储在 meta.data 中，列名格式：基因集名_UCell

# ============================================================
# 方法3：AUCell
# ============================================================
# 提取表达矩阵
expr_matrix <- GetAssayData(seurat_obj, layer = "counts")  # 使用raw counts

# 计算基因排名（每个细胞内基因按表达量排序）
cells_rankings <- AUCell_buildRankings(
  expr_matrix,
  nCores = 4,
  plotStats = FALSE
)

# 计算 AUC 分数
cells_AUC <- AUCell_calcAUC(
  gene_sets,               # 基因集
  cells_rankings,          # 细胞内基因排名
  aucMaxRank = ceiling(0.05 * nrow(cells_rankings))  # 用前5%基因计算AUC
)

# 提取分数矩阵
aucell_scores <- as.data.frame(t(getAUC(cells_AUC)))   # 细胞×基因集

# 将AUCell分数添加到元数据
seurat_obj$AUCell_Hypoxia <- aucell_scores[rownames(seurat_obj@meta.data), "Hypoxia"]
seurat_obj$AUCell_EMT     <- aucell_scores[rownames(seurat_obj@meta.data), "EMT"]

# ============================================================
# 方法比较：相关性分析
# ============================================================
method_compare <- seurat_obj@meta.data %>%
  select(
    AMS_Hypoxia1,     # AddModuleScore结果（注意有数字后缀）
    Hypoxia_UCell,    # UCell结果
    AUCell_Hypoxia    # AUCell结果
  ) %>%
  rename(
    AddModuleScore = AMS_Hypoxia1,
    UCell = Hypoxia_UCell,
    AUCell = AUCell_Hypoxia
  )

# 计算三种方法之间的相关性
cor_matrix <- cor(method_compare, method = "spearman", use = "pairwise.complete.obs")
print(round(cor_matrix, 3))    # 通常三种方法高度相关（r>0.8）

# 散点图比较
library(GGally)
GGally::ggpairs(method_compare,
  upper = list(continuous = wrap("cor", method = "spearman")),
  title = "三种基因集评分方法比较（缺氧通路）"
)

import decoupler as dc   # pip install decoupler
import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

adata = sc.read_h5ad("preprocessed.h5ad")

# 从 MSigDB 获取通路基因集（需要 pydeseq2 或直接下载gmt）
# 简单起见，手动定义基因集
gene_sets_dict = {
    "Hypoxia":       ["HIF1A", "VEGFA", "SLC2A1", "LDHA", "PGK1"],
    "Proliferation": ["MKI67", "TOP2A", "PCNA", "MCM2", "CDK4"],
    "EMT":           ["VIM", "CDH2", "SNAI1", "ZEB1", "FN1"]
}

# 转换为 decoupler 需要的 net DataFrame
net_rows = []
for pathway, genes in gene_sets_dict.items():
    for gene in genes:
        if gene in adata.var_names:    # 只保留数据中存在的基因
            net_rows.append({"source": pathway, "target": gene, "weight": 1.0})
net_df = pd.DataFrame(net_rows)

# ============================================================
# 方法1：ULM（单变量线性模型）
# ============================================================
dc.run_ulm(mat=adata, net=net_df, source="source", target="target",
           weight="weight", use_raw=False)
ulm_scores = pd.DataFrame(adata.obsm["ulm_estimate"],
                           index=adata.obs_names,
                           columns=adata.uns["ulm_estimate_names"])

# ============================================================
# 方法2：MLM（多变量线性模型）
# ============================================================
dc.run_mlm(mat=adata, net=net_df, source="source", target="target",
           weight="weight", use_raw=False)
mlm_scores = pd.DataFrame(adata.obsm["mlm_estimate"],
                           index=adata.obs_names,
                           columns=adata.uns["mlm_estimate_names"])

# ============================================================
# 方法比较图
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# ULM 分数
axes[0].scatter(ulm_scores["Hypoxia"], ulm_scores["Proliferation"],
                alpha=0.3, s=5, c="steelblue")
axes[0].set_xlabel("缺氧评分 (ULM)")
axes[0].set_ylabel("增殖评分 (ULM)")
axes[0].set_title("ULM 基因集评分")

# MLM 分数
axes[1].scatter(mlm_scores["Hypoxia"], mlm_scores["Proliferation"],
                alpha=0.3, s=5, c="coral")
axes[1].set_xlabel("缺氧评分 (MLM)")
axes[1].set_ylabel("增殖评分 (MLM)")
axes[1].set_title("MLM 基因集评分")

plt.tight_layout()
plt.savefig("method_comparison.pdf", dpi=300)

# 方法间相关性
from scipy.stats import spearmanr
r, p = spearmanr(ulm_scores["Hypoxia"], mlm_scores["Hypoxia"])
print(f"ULM vs MLM Spearman r = {r:.3f}, p = {p:.2e}")