1.6 单细胞差异分析 + 功能富集分析流程

# ================================================================================
# 单细胞差异分析 1: 两组差异分析 (FindMarkers)
# 版本: 优先 Seurat V5 (你日常使用的版本) / 兼容 V4
# 来源: Seurat V5 教程 6_DEGs_analysis.R (整理为可复用模板)
# 场景: 按细胞类型分层, 比较两组(如 GO vs WT)的差异基因
#
# ┌─ Seurat V5 差异分析要点 ──────────────────────────────────────────────┐
# │ 1. 安装 V5: remotes::install_github("satijalab/seurat", ref="seurat5") │
# │ 2. V5 的 assay 结构: RNA 和 SCT 是两个独立 assay, 用 DefaultAssay() 切换 │
# │ 3. SCT 数据必须 PrepSCTFindMarkers() 否则 FindMarkers 直接报错         │
# │ 4. V5 差异结果列名: avg_log2FC (V4 是 avg_logFC), 注意区分             │
# │ 5. V5 推荐 test.use: "wilcox"(默认)/ "MAST" / "DESeq2" / "negbinom"    │
# └───────────────────────────────────────────────────────────────────────┘
# ================================================================================

rm(list = ls())
library(Seurat)

# --------------------------------------------------------------------------------
# 0. 加载单细胞对象
#    需包含: celltype(细胞类型) + orig.ident/group(条件分组)
# --------------------------------------------------------------------------------
# load("~/data_analysis/seuratV5_test/scRNA_analysis/sce_cca_cellltype.RData")
# sce_cca <- readRDS("path/to/seurat_object.RDS")

## 检查 Seurat 版本 (V5 会显示 Version 5)
packageVersion("Seurat")     # [1] '5.1.0' 等

## V5 关键: 确认当前使用的 assay
DefaultAssay(sce_cca)        # 查看当前 assay (RNA 或 SCT)
# RNA 数据 → 用 RNA assay 直接 FindMarkers
# SCT 数据 → 必须先 PrepSCTFindMarkers (见模式2)
# 切换 assay:
# DefaultAssay(sce_cca) <- "SCT"

## ================= 模式1: 按细胞类型分层比较两组 =================
# 所有细胞类型中两组的差异
celltypes <- unique(sce_cca$celltype)
DEGs_celltype <- list()

for (i in 1:length(celltypes)) {
  data <- subset(sce_cca, celltype == celltypes[i])
  deg <- FindMarkers(data,
                     group.by = "orig.ident",
                     ident.1 = 'GO',     # 处理组(做分子)
                     ident.2 = "WT",     # 对照组
                     logfc.threshold = 0.25,
                     min.pct = 0.25)
  DEGs_celltype[[i]] <- deg
  names(DEGs_celltype)[i] <- celltypes[i]
}
# 查看单个细胞类型的差异结果
head(DEGs_celltype[[1]])
# p_val  avg_log2FC  pct.1  pct.2  p_val_adj

## ================= 模式2: SCT 数据的差异分析 (V5 重点) =================
# SCT标准化数据直接 FindMarkers 会报错:
# Object contains multiple models with unequal library sizes.
# Run `PrepSCTFindMarkers()` before running `FindMarkers()`.
# 原因: V5 的 SCT assay 对每个细胞用独立模型, 需先合并模型参数

# 必须先运行 (V5 独有, V4 无此步骤):
sce_sct_cca <- PrepSCTFindMarkers(sce_sct_cca)

# 方法A: 按 条件+cluster 组合分组比较
# sce_sct_cca$newgroup <- paste(sce_sct_cca$orig.ident, sce_sct_cca$seurat_clusters, sep = "_")
# Idents(sce_sct_cca) <- "newgroup"
# DEGs_0 <- FindMarkers(sce_sct_cca,
#                       ident.1 = "GO_0",    # 处理组的cluster 0
#                       ident.2 = "WT_0",    # 对照组的cluster 0
#                       logfc.threshold = 0.25,
#                       min.pct = 0.25)

## V5 差异结果列名注意:
# V5: p_val, avg_log2FC, pct.1, pct.2, p_val_adj
# V4: p_val, avg_logFC,  pct.1, pct.2, p_val_adj
#    (V5 把 avg_logFC 改为 avg_log2FC, 筛选时注意列名)

## ================= 3. 结果整理与保存 =================
# 合并所有细胞类型差异结果
DEGs_all <- do.call(rbind, lapply(names(DEGs_celltype), function(ct) {
  df <- DEGs_celltype[[ct]]
  df$celltype <- ct
  df$gene <- rownames(df)
  df
}))

# 筛选显著差异基因
library(dplyr)
# V5 列名是 avg_log2FC; 若用 V4 换成 avg_logFC
DEGs_sig <- DEGs_all %>%
  filter(p_val_adj < 0.05, abs(avg_log2FC) > 0.25)

# 保存
write.csv(DEGs_all, "result/scRNA_DEGs_all_celltype.csv", row.names = F)
write.csv(DEGs_sig, "result/scRNA_DEGs_sig_celltype.csv", row.names = F)
save(DEGs_celltype, DEGs_all, DEGs_sig, file = "result/scRNA_DEGs.RData")

# ================================================================================
# 单细胞差异分析 2: 多组差异分析 (FindAllMarkers 找各群 marker)
# 版本: 优先 Seurat V5 / 兼容 V4
# 来源: Seurat V4 流程 + KS 单细胞多组差异基因函数
# 场景: 为每个细胞群/cluster 找特异性 marker 基因
#
# V5 注意:
# - FindAllMarkers 在 V5 中用法不变, 但 SCT 数据需先 PrepSCTFindMarkers()
# - 结果列名: V5 是 avg_log2FC, V4 是 avg_logFC
# ================================================================================

rm(list = ls())
library(Seurat)
library(dplyr)

# --------------------------------------------------------------------------------
# 0. 加载单细胞对象 (需已聚类/注释)
# --------------------------------------------------------------------------------
# sce <- readRDS("path/to/seurat_object.RDS")

## ================= 1. FindAllMarkers 找每群 marker =================
# 对每个 cluster/细胞类型找特异性 marker (仅上调)
all.markers <- FindAllMarkers(scRNAdata,
                              logfc.threshold = 0.25,
                              min.pct = 0.25,
                              only.pos = TRUE)   # 只保留上调基因

# 结果结构: gene, cluster, avg_log2FC, pct.1, pct.2, p_val, p_val_adj
head(all.markers)

# 筛选显著 marker
markers_sig <- all.markers %>%
  filter(p_val_adj < 0.05, avg_log2FC > 0.5)

# 每个群 top10 marker
top10 <- markers_sig %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = avg_log2FC)
top10

# 保存
write.csv(markers_sig, "result/scRNA_markers_all.csv", row.names = F)
write.csv(top10, "result/scRNA_markers_top10.csv", row.names = F)

## ================= 2. 多组差异基因可视化(KS原创函数参考) =================
# KS_scRNA_multiVlnvo_plot: 单细胞多组差异基因 + 小提琴展示
# 来源: KS代码/5-KS原创函数合集/2023-8-25 单细胞多组差异基因函数
# 功能: 输入seurat对象+差异结果, 绘制多组差异基因小提琴图
# 核心参数: min.pct, logfc.threshold, test.use, ident.1, ident.2,
#           logFC_cut, top_gene, label_gene, point_color, box_color 等
if (FALSE) {
  source("单细胞多组差异基因函数.R")
  KS_scRNA_multiVlnvo_plot(
    Seurat_object = scRNAdata,
    DEGs_outdir = "result/",
    DEGs_list = TRUE,
    dges = DEGs_celltype,
    min.pct = 0.25, logfc.threshold = 0.25, test.use = "wilcox",
    group = "orig.ident", ident.1 = "GO", ident.2 = "WT",
    logFC_cut = 0.5, top_gene = TRUE
  )
}

# ================================================================================
# 单细胞差异分析 3: 差异基因可视化 (火山图/VlnPlot/热图)
# ================================================================================

rm(list = ls())
library(Seurat)
library(ggplot2)
library(dplyr)
library(ggrepel)

# 加载差异结果 (由 02/03 脚本生成)
# load("result/scRNA_DEGs.RData")
# DEGs_all, DEGs_sig, DEGs_celltype

## ================= 1. 火山图 =================
## 以某个细胞类型的差异结果为例
df <- DEGs_celltype[["T cell"]]    # 替换为你的细胞类型
df$gene <- rownames(df)

# 标记显著上下调
df$regulated <- "NS"
df$regulated[df$p_val_adj < 0.05 & df$avg_log2FC > 0.25] <- "Up"
df$regulated[df$p_val_adj < 0.05 & df$avg_log2FC < -0.25] <- "Down"

# 挑选标注 top 基因
label_genes <- df %>%
  filter(regulated != "NS") %>%
  arrange(p_val_adj) %>%
  head(10)

p <- ggplot(df, aes(x = avg_log2FC, y = -log10(p_val_adj), color = regulated)) +
  geom_point(alpha = 0.6, size = 1.2) +
  geom_point(data = label_genes, size = 2.5, shape = 1, color = "black") +
  geom_text_repel(data = label_genes, aes(label = gene), size = 3, color = "black") +
  scale_color_manual(values = c(Up = "#e64b35", Down = "#4dbbd6", NS = "#999999")) +
  geom_vline(xintercept = c(-0.25, 0.25), linetype = "dashed", color = "grey50") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
  labs(x = "avg_log2FC", y = "-log10(p_val_adj)",
       title = "DEGs Volcano Plot") +
  theme_bw() +
  theme(panel.grid = element_blank(), legend.position = "right")
p
ggsave("result/volcano_plot.png", p, width = 7, height = 6, dpi = 300)

## ================= 2. VlnPlot 小提琴(单基因分组表达) =================
# 展示某个 marker 基因在两组中的表达差异
# 需加载 seurat 对象
# sce_cca$group <- sce_cca$orig.ident  # 确保有分组列
p2 <- VlnPlot(sce_cca,
              features = "TP53",          # 替换为目标基因
              group.by = "orig.ident",    # 分组
              split.by = "celltype",      # 按细胞类型分面(可选)
              pt.size = 0.1) +
  theme(legend.position = "none")
p2
ggsave("result/vlnplot_gene.png", p2, width = 6, height = 5, dpi = 300)

## ================= 3. DoHeatmap 差异基因热图 =================
# Top 差异基因在细胞中表达热图
top_genes <- df %>%
  filter(regulated != "NS") %>%
  arrange(p_val_adj) %>%
  head(30) %>%
  pull(gene)

p3 <- DoHeatmap(sce_cca,
                features = top_genes,
                group.by = "celltype",
                size = 3) +
  scale_fill_gradient2(low = "#4dbbd6", mid = "white", high = "#e64b35")
p3
ggsave("result/deg_heatmap.png", p3, width = 8, height = 8, dpi = 300)

# ================================================================================
# 单细胞功能富集 1: ORA 过代表分析 (GO / KEGG / Hallmark)
# 来源: Seurat V4 流程富集段 + clusterProfiler
# 输入: 差异基因列表 (由 01_差异分析 得到)
# ================================================================================

rm(list = ls())
library(clusterProfiler)
library(org.Hs.eg.db)   # 人; 小鼠用 org.Mm.eg.db
library(GSEABase)       # 读取 gmt
library(ggplot2)
library(dplyr)

# --------------------------------------------------------------------------------
# 0. 加载差异基因列表
# --------------------------------------------------------------------------------
# load("result/scRNA_DEGs.RData")
# 从某个细胞类型的显著差异基因取基因名 (SYMBOL)
DEG <- DEGs_celltype[["T cell"]]       # 替换
DEG <- DEG[DEG$p_val_adj < 0.05 & abs(DEG$avg_log2FC) > 0.25, ]
DEG_genes <- rownames(DEG)
head(DEG_genes)

## ================= 1. GO 富集分析 =================
# ont: BP(生物学过程) / MF(分子功能) / CC(细胞组分) / ALL
ego_BP <- enrichGO(gene = DEG_genes,
                   OrgDb = "org.Hs.eg.db",
                   keyType = "SYMBOL",        # 基因ID类型
                   ont = "BP",
                   pvalueCutoff = 1,          # 输出全部, 后续筛选
                   qvalueCutoff = 1,
                   readable = TRUE)
head(ego_BP@result)

# 若需要 ENTERZID 输入(另一模式):
# ego <- enrichGO(gene = DEG_entrez, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP")

## ================= 2. KEGG 富集分析 =================
# KEGG 需要 ENTREZID, 先转换
genelist <- bitr(DEG_genes,
                 fromType = "SYMBOL",
                 toType = "ENTREZID",
                 OrgDb = org.Hs.eg.db)
head(genelist)

ekegg <- enrichKEGG(gene = genelist$ENTREZID,
                    organism = 'hsa',         # 人; 小鼠 'mmu'
                    pvalueCutoff = 1,
                    qvalueCutoff = 1)
# 将 ENTREZID 转回 SYMBOL 方便阅读
ekegg <- DOSE::setReadable(ekegg, OrgDb = 'org.Hs.eg.db', keyType = 'ENTREZID')
head(ekegg@result)

## ================= 3. Hallmark 基因集富集 (自定义 gmt) =================
# 从 MSigDB 下载 h.all.v*.Hs.symbols.gmt
# geneset <- read.gmt("data/MsigDB/h.all.v2023.1.Hs.symbols.gmt")
# my_path <- enricher(gene = DEG_genes,
#                     pvalueCutoff = 1, qvalueCutoff = 1,
#                     TERM2GENE = geneset,
#                     minGSSize = 1, maxGSSize = 10000)

## ================= 4. 保存结果 =================
# 筛选显著富集
ego_BP_sig <- ego_BP@result[ego_BP@result$p.adjust < 0.05, ]
ekegg_sig <- ekegg@result[ekegg@result$p.adjust < 0.05, ]

write.csv(ego_BP@result, "result/enrich_GO_BP.csv", row.names = F)
write.csv(ekegg@result, "result/enrich_KEGG.csv", row.names = F)
save(ego_BP, ekegg, file = "result/enrich_ORA.RData")

# ================================================================================
# 单细胞功能富集 2: GSEA 基因集富集分析
# 来源: KS_GSEA 通用函数 (clusterProfiler 版) + Seurat V4 流程
# 输入: 所有基因的 log2FC 排序列表 (不限于显著基因)
# ================================================================================

rm(list = ls())
library(clusterProfiler)
library(org.Hs.eg.db)
library(GSEABase)
library(ggplot2)

# --------------------------------------------------------------------------------
# 0. 加载差异分析结果 (全基因, 含 logFC)
# --------------------------------------------------------------------------------
# load("result/scRNA_DEGs.RData")
# DEG 需含: 所有基因 + avg_log2FC 列 (行名为基因名)
DEG <- DEGs_celltype[["T cell"]]
DEG <- DEG[!is.na(DEG$avg_log2FC), ]

# 去重(避免基因名重复)
DEG <- DEG[!duplicated(rownames(DEG)), ]

# 构造排序基因列表 (按 avg_log2FC 排序)
geneList <- DEG$avg_log2FC
names(geneList) <- rownames(DEG)
geneList <- sort(geneList, decreasing = T)   # 降序
head(geneList)

## ================= 1. 读取 gmt 基因集 =================
# MSigDB Hallmark (需自行下载)
geneset <- read.gmt("data/MsigDB/h.all.v2023.1.Hs.symbols.gmt")
# 或用 KEGG legacy:
# geneset <- read.gmt("data/MsigDB/c2.cp.kegg_legacy.v2023.1.Hs.symbols.gmt")

## ================= 2. 运行 GSEA =================
egmt <- GSEA(geneList,
             TERM2GENE = geneset,
             pvalueCutoff = 1,    # 输出全部, 后续筛选
             minGSSize = 1,
             maxGSSize = 10000)
head(egmt@result)
# 列: ID, Description, setSize, enrichmentScore, NES, pvalue, p.adjust

# 筛选显著通路
gsea_sig <- egmt@result[egmt@result$p.adjust < 0.05, ]
gsea_sig <- gsea_sig[order(gsea_sig$NES, decreasing = T), ]
head(gsea_sig)

## ================= 3. 保存 =================
write.csv(egmt@result, "result/enrich_GSEA_all.csv", row.names = F)
write.csv(gsea_sig, "result/enrich_GSEA_sig.csv", row.names = F)
save(egmt, file = "result/enrich_GSEA.RData")

# ================================================================================
# 单细胞功能富集 3: 富集结果可视化
# dotplot / barplot / gseaplot2 / 分组对比
# ================================================================================

rm(list = ls())
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(dplyr)

# 加载富集结果
# load("result/enrich_ORA.RData")    # ego_BP, ekegg
# load("result/enrich_GSEA.RData")   # egmt

## ================= 1. ORA 结果可视化 =================
# GO BP 气泡图 (横轴为 geneRatio, 颜色为 p.adjust)
p1 <- dotplot(ego_BP, showCategory = 15, font.size = 8) +
  ggtitle("GO Biological Process") +
  theme(axis.text.y = element_text(size = 7))
p1
ggsave("result/dotplot_GO_BP.png", p1, width = 8, height = 6, dpi = 300)

# KEGG 条形图
p2 <- barplot(ekegg, showCategory = 15, font.size = 8) +
  ggtitle("KEGG Pathway")
p2
ggsave("result/barplot_KEGG.png", p2, width = 7, height = 5, dpi = 300)

# GO 三个本体 (BP/MF/CC) 组合图
if (FALSE) {
  ego_CC <- enrichGO(gene = DEG_genes, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL",
                     ont = "CC", pvalueCutoff = 1, qvalueCutoff = 1)
  ego_MF <- enrichGO(gene = DEG_genes, OrgDb = "org.Hs.eg.db", keyType = "SYMBOL",
                     ont = "MF", pvalueCutoff = 1, qvalueCutoff = 1)
  library(patchwork)
  p_all <- (barplot(ego_BP, showCategory = 10) + ggtitle("BP")) /
    (barplot(ego_CC, showCategory = 10) + ggtitle("CC")) /
    (barplot(ego_MF, showCategory = 10) + ggtitle("MF"))
  ggsave("result/GO_BP_CC_MF.png", p_all, width = 8, height = 12, dpi = 300)
}

## ================= 2. GSEA 结果可视化 =================
# 单个通路富集图 (running enrichment score)
p3 <- gseaplot2(egmt, geneSetID = 1,     # 第几条通路
                title = egmt@result$Description[1],
                pvalue_table = TRUE)
p3
ggsave("result/gseaplot_single.png", p3, width = 7, height = 6, dpi = 300)

# 多条通路富集图
if (FALSE) {
  p4 <- gseaplot2(egmt, geneSetID = 1:6, pvalue_table = FALSE,
                  color = RColorBrewer::brewer.pal(6, "Set1"))
  ggsave("result/gseaplot_multiple.png", p4, width = 8, height = 6, dpi = 300)
}

## ================= 3. GSEA 通路排序点图 (高颜值) =================
# 类似 Science 风格: x=通路排序, y=NES, 点大小=setSize, 颜色=显著性
gsea_df <- egmt@result
gsea_df <- gsea_df[order(gsea_df$NES, decreasing = T), ]
gsea_df$ID <- factor(gsea_df$ID, levels = gsea_df$ID)
gsea_df$xlab <- 1:nrow(gsea_df)

p5 <- ggplot(gsea_df, aes(x = xlab, y = NES)) +
  geom_point(aes(size = setSize, color = -log10(pvalue)), shape = 16) +
  scale_color_gradient(low = "grey80", high = "#e64b35") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  labs(x = "Pathways", y = "NES", color = "-log10(p)") +
  theme_classic(base_size = 12) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
p5
ggsave("result/gsea_dotplot_NES.png", p5, width = 8, height = 5, dpi = 300)