下面的代码还未测试 ，后续测试

library(Seurat)
library(tidyverse)
## 1 常规富集分析----
# 常规富集基于 top DEG
library(clusterProfiler)
library(enrichplot)

### 1.1、获得DEG----
# 同前，只需要要把默认 idents 从 seurat_clusters 改成 cell_type
Idents(sc_pbmc_int3)
Idents(sc_pbmc_int3) = "cell_type" # 把celltype设置为主要的分类标签

#
cell.markers <- FindAllMarkers(object = sc_pbmc_int3,
                               only.pos = FALSE, # 是否只保留表达相对上调的基因，设置FALSE则会保留下调的
                               test.use = "wilcox", # 默认使用 wilcox 非参数检验，其它选项可以查看说明
                               slot = "data", # 需要注意的是，默认使用 data，而不是 counts
                               min.pct = 0.25, # 设置表达比例的阈值，没有统一标准，差异基因很多的情况下可以把阈值调高，差异基因有1000个够用
                               logfc.threshold = 0.5 # 设置 log2FC 即差异倍率的阈值，没有统一标准，差异基因很多的情况下可以把阈值调高
)

cell.markers$gene2 = rownames(cell.markers)


### 1.2 准备基因列表----
colnames(cell.markers)
table(cell.markers$cluster)

# 筛选DEG
lsxx = cell.markers %>%
  filter(cluster == "CD4+ T cells") %>%
  filter(p_val_adj < 0.05) %>%
  filter(avg_log2FC > 0) %>% # 多个条件的筛选也可以使用 & 和 | 
  # filter( abs(avg_log2FC) > 1 )
  arrange(desc(avg_log2FC))

# symbol 转 ENTREZID
annot = bitr(lsxx$gene, fromType = "SYMBOL", toType = c("ENTREZID"),
             OrgDb = "org.Hs.eg.db", drop = F) %>%
  drop_na() %>%
  distinct(SYMBOL, .keep_all = T)

lsxx = lsxx %>%
  mutate(ENTREZID = annot$ENTREZID[match(gene, annot$SYMBOL)]) %>%
  na.omit(ENTREZID)


### 1.3 富集分析----
go_data <- enrichGO(OrgDb = "org.Hs.eg.db", # 种属基因数据
                    gene = lsxx$ENTREZID,
                    ont = "ALL", # CC BP MF
                    # pvalueCutoff = 0.1, # 设置校准P值筛选阈值
                    # qvalueCutoff = 0.1, # 设置Q值筛选阈值
                    readable = TRUE # 自动将 ID 转换为 SYMBOL
) %>%
  clusterProfiler::simplify(.) %>% # 删减重复条目，可选，相似结果比较多时好用
  enrichplot::pairwise_termsim() # emap 网络图需要做此转换

#
kegg_data <- enrichKEGG(gene = lsxx$ENTREZID,
                        keyType = "kegg",
                        organism = 'hsa', # 种属设置
                        # pvalueCutoff = 1, # 设置1意味着不筛选，保留所有结果，可以在后续处理中筛选
                        # qvalueCutoff = 1,
                        use_internal_data = F) %>% # 是否使用本地数据，默认 F 
  setReadable(., OrgDb = 'org.Hs.eg.db', keyType = 'ENTREZID') %>% # 将 ENTRZID 转换成 SYMBOL
  pairwise_termsim() # emap 网络图需要做此转换

#
save(go_data, kegg_data, kegg_category, file = '1.3-常规富集分析.Rdata')

rm(annot, go_data, kegg_data, kegg_category)


### 1.4 可视化


### 1.5 批量富集
unique(cell.markers$cluster)

enrich_dis = lapply(unique(cell.markers$cluster), function(xx) {
  # xx = unique( cell.markers$cluster )[2]

  lsxx = cell.markers %>%
    filter(cluster == xx) %>%
    filter(p_val_adj < 0.05) %>%
    filter(avg_log2FC > 0)  # 多个条件的筛选也可以使用 & 和 | 
  #
  annot = bitr(lsxx$gene, fromType = "SYMBOL", toType = c("ENTREZID"),
               OrgDb = "org.Hs.eg.db", drop = F) %>%
    distinct(SYMBOL, .keep_all = T)
  #
  lsxx = lsxx %>%
    mutate(ENTREZID = annot$ENTREZID[match(gene, annot$SYMBOL)]) %>%
    na.omit(ENTREZID)

  #
  go_data <- enrichGO(OrgDb = "org.Hs.eg.db", # 种属基因数据
                      gene = lsxx$ENTREZID,
                      ont = "ALL",
                      # pvalueCutoff = 0.1, # 设置校准P值筛选阈值
                      # qvalueCutoff = 0.1, # 设置Q值筛选阈值
                      readable = TRUE # 自动将 ID 转换为 SYMBOL
  ) %>%
    clusterProfiler::simplify(.) %>% # 删减重复条目，可选，相似结果比较多时好用
    enrichplot::pairwise_termsim() # emap 网络图需要做此转换

  #
  kegg_data <- enrichKEGG(gene = lsxx$ENTREZID,
                          keyType = "kegg",
                          organism = 'hsa', # 种属设置
                          # pvalueCutoff = 1, # 设置1意味着不筛选，保留所有结果，可以在后续处理中筛选
                          # qvalueCutoff = 1,
                          use_internal_data = F) %>% # 是否使用本地数据，默认 F 
    setReadable(., OrgDb = 'org.Hs.eg.db', keyType = 'ENTREZID') %>% # 将 ENTRZID 转换成 SYMBOL
    pairwise_termsim() # emap 网络图需要做此转换

  #
  list(GO = go_data, KEGG = kegg_data)

})

#
names(enrich_dis) = unique(cell.markers$cluster)


rm(cell.markers, annot, lsxx, go_data, kegg_data, kegg_category, enrich_dis)


## 2 GSEA 富集分析----
Idents(sc_pbmc_int3) = "cell_type" # 把celltype设置为主要的分类标签
# 这种方法不大推荐，直接使用所有基因的 logFC 进行 GSEA 富集有较大的风险，如果 logFC 分布不均会导致明显的结果偏倚

### 2.1 获得所有基因的logFC----
all.genes <- FindAllMarkers(object = sc_pbmc_int3,
                            only.pos = F, # 是否只保留表达相对上调的基因，设置FALSE则会保留下调的
                            test.use = "wilcox", # 默认使用 wilcox 非参数检验，其它选项可以查看说明
                            slot = "data", # 需要注意的是，默认使用 data，而不是 counts
                            min.pct = 0.01, # 设置更低可以保留更多基因，但是容易得到异常的差异基因，比如表达比例非常低但倍率非常高的情况
                            logfc.threshold = 0 # 设置为 0，保证结果里面有所有基因
)

#
save(all.genes, file = "2.1-不同细胞所有基因logFC.rda")

table(all.genes$cluster)
table(all.genes$pct.1 >= 0.01)

hist(all.genes$avg_log2FC)

colnames(all.genes)

ggplot(all.genes, aes(x = avg_log2FC, y = -log10(p_val_adj))) +
  geom_point(size = 0.1, alpha = 0.5)


### 2.2 准备基因列表----
lsxx = all.genes %>%
  filter(cluster == "CD4+ T cells") %>%
  arrange(desc(avg_log2FC)) # logFC 从大大小，逆序排列

hist(lsxx$avg_log2FC)

# 转换成 ENTREZ
annot = bitr(lsxx$gene, fromType = "SYMBOL", toType = c("ENTREZID"),
             OrgDb = "org.Hs.eg.db", drop = F) %>%
  distinct(SYMBOL, .keep_all = T)

lsxx = lsxx %>%
  mutate(ID = annot$ENTREZID[match(gene, annot$SYMBOL)]) %>%
  na.omit(ID)

# 得到命名向量，数值为从大到小排列的 logFC，向量名为 ENTREZID
lsxx = setNames(lsxx$avg_log2FC, lsxx$ID)

head(lsxx); tail(lsxx)


### 2.3 富集分析----
gs_go <- gseGO(lsxx,
               ont = "ALL",
               OrgDb = "org.Hs.eg.db",
               eps = 0) %>%
  clusterProfiler::simplify() %>%
  setReadable(., OrgDb = 'org.Hs.eg.db', keyType = 'ENTREZID')  # 将 ENTRZID 转换成 SYMBOL，GSEA 做这个转换可能导致后续 ridgeplot 出错 

#
gs_kg <- gseKEGG(lsxx,
                 keyType = "kegg",
                 organism = 'hsa',
                 # pvalueCutoff = 1,
                 # use_internal_data = TRUE, # 如果提示网络连接失败，使用这条，用本地的数据库做富集，需要安装本地数据库，具体请自行搜索
                 eps = 0) %>%
  setReadable(., OrgDb = 'org.Hs.eg.db', keyType = 'ENTREZID')  # 将 ENTRZID 转换成 SYMBOL，GSEA 做这个转换可能导致后续 ridgeplot 出错
#
save(gs_go, gs_kg, file = '2.3-GSEA富集分析.Rdata')

rm(lsxx, annot, gs_go, gs_kg)

### 2.4 可视化
# 富集分析自定义可视化

### 2.5 批量富集（方法同前，使用for循环或lapply等，需要调整）----

## 3、基于 UCell 等方法的基因集评分----
# 这种方法的优势在于同时整合了ssGSEA外的其他评分方法
# 并且生成的得分数据可以直接使用 Seurat 自带函数进行分析和可视化

library(irGSEA)
library(ComplexHeatmap)
Idents(sc_pbmc_int3) = "cell_type" # 把celltype设置为主要的分类标签

### 3.1 分析基因集评分----
sc_vec_fun <- irGSEA.score(object = sc_pbmc_int3,
                           assay = "RNA",
                           slot = "data",
                           seeds = 124,
                           ncores = 10, # 线程数，内存不够的把线程设低点
                           min.cells = 3, min.feature = 0,
                           custom = F, geneset = NULL, # 是否使用自定义的基因集
                           msigdb = T,
                           species = "Homo sapiens", # 设置种属
                           category = "H", # 设置要使用的 MsigDB 数据库子库，默认 H，即 Hallmark
                           subcategory = NULL,
                           geneid = "symbol", # 输入的基因名类型
                           method = c("AUCell", "UCell", "singscore", "ssgsea"), # 有多种方法可选
                           aucell.MaxRank = NULL, ucell.MaxRank = NULL,
                           kcdf = 'Gaussian')

# 新增打分数据：4 other assays present: AUCell, UCell, singscore, ssgsea
sc_vec_fun


### 3.2 结果分析（自带函数）----
result.def <- irGSEA.integrate(object = sc_vec_fun,
                               group.by = "cell_type", # 基于什么分组计算差异基因集/功能，可以改成别的属性，比如 vaccine
                               method = c("AUCell", "UCell", "singscore", "ssgsea")) # 

# 其中 RRA 表格中是对多种方法使用Robust Rank Aggregation整合得到的差异基因集/功能
names(result.def)

### 3.3 结果分析（Seurat函数）----
# 方法同单细胞分析
table(sc_vec_fun$cell_type)
table(sc_vec_fun$group)

DefaultAssay(sc_vec_fun) = "UCell"
#
sc_def = FindMarkers(subset(sc_vec_fun, cell_type == "CD8+ T cells"),
                     ident.1 = "elevated group (E)", ident.2 = "control group (C)",
                     group.by = "group",
                     logfc.threshold = 0 # 因为得到的倍率都比较小，建议设置 0
)


### 3.4 可视化----

# a、自带可视化函数

# 条形图
irGSEA.barplot(object = result.def, method = c("UCell"))

# 自带山脊图
irGSEA.ridgeplot(object = sc_vec_fun, method = "UCell", show.geneset = "HALLMARK-INFLAMMATORY-RESPONSE")

# 自带散点图
irGSEA.density.scatterplot(object = sc_vec_fun,
                           method = "UCell",
                           show.geneset = c("HALLMARK-INFLAMMATORY-RESPONSE"),
                           reduction = "umap")


# b、Seurat 可视化

# 因为是 Seurat 对象，还可以用 FeaturePlot 和 DotPlot 等可视化
# 注意：如果使用多种方法进行评分，会导致 Assay 下多个矩阵中有同样的基因集/功能，导致可视化报错
names(sc_vec_fun@assays)

# 可以指定用哪个Assay可视化
DefaultAssay(sc_vec_fun) = "UCell"

# FeaturePlot
FeaturePlot(sc_vec_fun, slot = "data",
            features = "HALLMARK-INFLAMMATORY-RESPONSE", reduction = "umap",
            order = T, pt.size = 0.8) +
  scale_color_gradient(low = "grey", high = "orangered3")

#
FeaturePlot(sc_vec_fun, slot = "data",
            features = "HALLMARK-INFLAMMATORY-RESPONSE", reduction = "umap",
            order = T, pt.size = 0.5) +
  scale_color_gradient2(low = "olivedrab", mid = "grey", high = "orangered3", midpoint = 0.14)


# DotPlot
head(result.def$UCell$Name)

DotPlot(sc_vec_fun, features = head(result.def$UCell$Name),
        assay = "UCell", group.by = "cell_type") +
  RotatedAxis()