定义：GSEA（Gene Set Enrichment Analysis）是一种基于基因集的富集分析方法， 用来评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势，从而判断其对表型的贡献

基本原理： 使用预定义的基因集 (可以是GO注释、MsigDB的注释或其它符合格式的基因集定义)，将基因按照在两类样本中的差异表达程度排序，然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。基因集合富集分析检测基因集合而不是单个基因的表达变化，因此可以包含这些细微的表达变化，预期得到更为理想的结果

与GO\KEGG差异基因富集分析区别：差异基因富集分析是先筛选差异基因，再判断差异基因在哪些注释的通路存在富集；这涉及到阈值的设定，存在一定主观性并且只能用于表达变化较大的基因，即我们定义的显著差异基因。而GSEA则不局限于差异基因，从基因集的富集角度出发，理论上更容易囊括细微但协调性的变化对生物通路的影响

GSEA-1

🌟 GSEA中关键概念
1）ES（Enrichment Score）：富集得分 ES反应基因集成员s在排序列表L的两端富集的程度。计算方式是，从基因集L的第一个基因开始，计算一个累计统计值。当遇到一个落在s里面的基因，则增加统计值。遇到一个不在s里面的基因，则降低统计值。 每一步统计值增加或减少的幅度与基因的表达变化程度（fold-change值）是相关的。富集得分ES最后定义为最大的峰值。正值ES表示基因集在列表的顶部富集，负值ES表示基因集在列表的底部富集。 p-value用来评估富集得分(ES)的显著性，通过排列检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。
2）NES（Normalized Enrichment Score）：标准化富集得分 每个基因子集s计算得到的ES根据基因集的大小进行标准化得到标准化富集得分Normalized Enrichment Score (NES)。随后会针对NES计算假阳性率FDR。
3）Leading-edge subset：领头基因亚集 对富集贡献最大的基因成员
一般认为|NES|>1，p-value<0.05，FDR<0.25的通路是显著富集的。 |NES|值越大，FDR值就越小，说明分析的结果可信度越高。

🌟 MSigDB数据库是GSEA官网提供的基因集数据库，其中包含了9种人的分类模式：
H：包含了由多个已知的基因集构成的超基因集，每个H 类别的基因集都对应多个基础的其他类别的基因集
C1：包含人类每条染色体上的不同cytoband区域对应的基因集合。根据不染色体编码进行二级分类
C2：包含了在线通路数据库、PubMed 出版物和领域专家知识的精选基因集
C3：基于对 miRNA 种子序列和预测的转录因子结合位点的基因靶点预测的调控靶基因组
C4：包含了计算机软件预测出来的基因集合，主要是和癌症相关的基因
C5：包含了由相同GO术语注释的基因集
C6：包含了致癌特征基因集：直接从来自癌症基因扰动的微阵列基因表达数据中定义
C7：包含了免疫特征基因集：代表免疫系统内的细胞状态和扰动
C8：包含了细胞类型特征基因组：从人体组织单细胞测序研究中确定的簇标记中收集

1. 使用包进行GSEA分析

1.1 读取差异表达的数据

# 加载所需的R包
library(fgsea)          # GSEA分析主程序
library(data.table)     # 数据处理
library(ggplot2)        # 画图处理
library(dplyr)          # 数据处理
library(msigdb)         # 包含基因集合，通常和GSEA分析共同使用
library(GSEABase)       # 可以提供GSEA基础结构和函数,也会被其他包调用


# 读取差异表达的数据，做排序和权重（一般是根据log2FC来作为排序依据）
head(DEG_DESeq2)
##              baseMean log2FoldChange     lfcSE      stat       pvalue         padj
## VCX2         4.848175       4.720165 0.6353906  7.428762 1.096188e-13 2.998016e-12
## AL590681.1   7.142233       4.616208 0.5637220  8.188803 2.638373e-16 9.854187e-15
## LINC00383    9.759471       4.402470 0.4719261  9.328727 1.071498e-20 6.688977e-19
## CARTPT     130.886246       4.301951 0.3291537 13.069734 4.903700e-39 4.891201e-36
## CLDN6      736.600962       4.240786 0.3126447 13.564236 6.525606e-42 1.443294e-38
## GAGE2A       9.285152       4.182749 0.7523736  5.559404 2.706976e-08 3.263125e-07

# 得到一个排序权重变量，一维有名称变量，后续分析要用
FCgenelist <- DEG_DESeq2$log2FoldChange                 # 把数值保存到一维向量中
names(FCgenelist) <- row.names(DEG_DESeq2)          # 对一维向量数值命名

1.2 准备基因集

# 获取msigdb数据集, 指定物种和基因标识符类型 ,如果第一次运行会下载一个大的缓存
msigdb.hs  <-  getMsigdb(org = 'hs',id = c("SYM", "EZID"))

# 追加KEGG集合（默认的msigdb没有KEGG）
msigdb.hs  <-  appendKEGG(msigdb.hs)

# 如果想查看集合,或者集合的子集可以使用以下函数
listCollections(msigdb.hs)
# [1] "c1" "c2" "c3" "c4" "c5" "c6" "c7" "c8" "h" 
listSubCollections(msigdb.hs)
##  [1] "CGP"             "CP:BIOCARTA"     "CP:PID"          "CP"              "CP:WIKIPATHWAYS" "MIR:MIRDB"      
##  [7] "MIR:MIR_Legacy"  "TFT:GTRD"        "TFT:TFT_Legacy"  "CGN"             "CM"              "HPO"            
## [13] "IMMUNESIGDB"     "VAX"             "CP:REACTOME"     "GO:BP"           "GO:CC"           "GO:MF"          
## [19] "CP:KEGG" 

# 这里使用msigdb 的"h"子集
hallmarks  <-  subsetCollection(msigdb.hs, 'h')

# 获取基因集合
msigdb_ids <- geneIds(hallmarks)

1.3 富集分析并作图

# 富集分析
fgseaRes <- fgsea(pathways = msigdb_ids, 
                  stats = FCgenelist,
                  eps = 0.0,
                  minSize = 15,
                  maxSize = 500)

head(fgseaRes[order(pval), ])
##                                       pathway        pval       padj   log2err        ES      NES size                              leadingEdge
## 1:                   HALLMARK_APICAL_JUNCTION 0.003223129 0.04512381 0.4317077 0.3835874 1.964507   36    CLDN6,ACTG2,ACTN2,FLNC,MYL9,CNTN1,...
## 2:                        HALLMARK_MYOGENESIS 0.012002671 0.07468170 0.3807304 0.2827196 1.723659   70      CASQ2,ACTN2,DES,MYH11,LDB3,PYGM,...
## 3:                      HALLMARK_ADIPOGENESIS 0.016003222 0.07468170 0.3524879 0.4362824 1.781732   18 ADIPOQ,FABP4,MYLK,HSPB8,CIDEA,SORBS1,...
## 4:                 HALLMARK_KRAS_SIGNALING_UP 0.029782360 0.10423826 0.2820134 0.3568450 1.667048   27 SNAP25,NAP1L2,CIDEA,PCP4,CPE,SPARCL1,...
## 5: HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION 0.128258603 0.33939394 0.1221079 0.2263475 1.352838   64      SFRP1,MYLK,MYL9,TAGLN,MGP,CRLF1,...
## 6:                       HALLMARK_COAGULATION 0.145454545 0.33939394 0.1238422 0.3165289 1.319653   19        APOA1,TF,FGG,LEFTY2,APOC3,FGA,...

# 绘制富集图
plotEnrichment(msigdb_ids[["HALLMARK_APICAL_JUNCTION"]],
               FCgenelist) + labs(title="HALLMARK APICAL JUNCTION")

# 筛选最top的几个通路并绘制富集图
topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n = 10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n = 10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))

plotGseaTable(msigdb_ids[topPathways], 
              FCgenelist, 
              fgseaRes, 
              gseaParam = 0.5)

# 只展示独立的通路，可以使用collapsePathways函数去冗余
collapsedPathways <- collapsePathways(fgseaRes[order(pval)][1:20], msigdb_ids, FCgenelist)
mainPathways <- fgseaRes[pathway %in% collapsedPathways$mainPathways][order(-NES), pathway]
plotGseaTable(msigdb_ids[mainPathways], 
              FCgenelist, 
              fgseaRes, 
              gseaParam = 0.5)

# 以文本格式保存结果
fwrite(fgseaRes, file = "fgseaRes.txt", sep = "\t", sep2 = c("", " ", ""))

GSEA-2

2. 使用包进行GSEA分析

2.1 富集分析

# 加载所需的R包
library(org.Hs.eg.db) # human的OrgDB
library(clusterProfiler)

# ID转化
gene_entrezid <- bitr(geneID = rownames(DEG_DESeq2), 
                      fromType = "SYMBOL", 
                      toType = "ENTREZID", # 转成ENTREZID
                      OrgDb = "org.Hs.eg.db"
)

head(gene_entrezid)
##      SYMBOL  ENTREZID
## 1      VCX2     51480
## 3 LINC00383 103689913
## 4    CARTPT      9607
## 5     CLDN6      9074
## 6    GAGE2A    729447
## 8     XAGE5    170627

gene_entrezid$logFC <- DEG_DESeq2$log2FoldChange[match(gene_entrezid$SYMBOL, rownames(DEG_DESeq2))]
genelist = gene_entrezid$logFC
names(genelist) = gene_entrezid$ENTREZID 

library(msigdbr)

m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, entrez_gene)
head(m_t2g)
## # A tibble: 6 × 2
##   gs_name               entrez_gene
##   <chr>                       <int>
## 1 HALLMARK_ADIPOGENESIS          19
## 2 HALLMARK_ADIPOGENESIS       11194
## 3 HALLMARK_ADIPOGENESIS       10449
## 4 HALLMARK_ADIPOGENESIS          33
## 5 HALLMARK_ADIPOGENESIS          34
## 6 HALLMARK_ADIPOGENESIS          35

gsea_res <- GSEA(genelist, 
                 TERM2GENE = m_t2g,
                 minGSSize = 10,
                 maxGSSize = 500,
                 pvalueCutoff = 1,
                 pAdjustMethod = "BH"
)
## preparing geneSet collections...
## GSEA analysis...
## leading edge analysis...
## done...

# 第一个条目的所有基因
gsea_res[[gsea_res$ID[[1]]]]

2.2 可视化

2.2.1 峰峦图（X轴是核心基因的表达倍数变化，正值表示上调，负值表示下调）

library(enrichplot)
library(ggplot2)

ridgeplot(gsea_res,
          showCategory = 20,
          fill = "pvalue", #填充色 "pvalue", "p.adjust", "qvalue" 
          core_enrichment = TRUE,#是否只使用 core_enriched gene
          label_format = 30,
          orderBy = "NES",
          decreasing = FALSE
) +
  theme(axis.text.y = element_text(size=8))

2.2.2 gseadist（展示基因集的logFC分布）

ids <- gsea_res@result$ID[1:5]

gseadist(gsea_res,
         IDs = ids,
         type = "density" # boxplot
) +
  theme(legend.direction = "vertical")

GSEA-3

2.2.3 gsearank（展示基因的排序以及富集分数的变化）

gsearank(gsea_res,
         geneSetID = 1 # 要展示的基因集
)

2.2.4 gseaplot

gseaplot函数可以画两个图：ES或者ranked-gene-list，通过参数by设置，默认是两个图都画出来，如果by="runningScore"，则是画出ES的图，如果是by = "preranked",则是画出ranked gene list的图

gseaplot(gsea_res, geneSetID = 1, by = "runningScore", 
              title = gsea_res$Description[1])

gseaplot(gsea_res, geneSetID = 1, by = "preranked", 
              title = gsea_res$Description[1]) + 
  theme(plot.title = element_text(size = 10, color = "blue"))

# 如果两个子图都画的话返回的是一个ggplist对象，此时如果要修改图形细节，可以使用取子集的方法提取其中的子图形，此时的子图形是ggplot对象，又可以使用ggplot2语法修改了
p <- gseaplot(gsea_res, geneSetID = 1, title = gsea_res$Description[1])
p

#取子集进行修改
p[[1]] <- p[[1]]+theme(plot.title = element_text(size = 6))
p

GSEA-4

2.2.5 gseaplot2

# 默认subplots = 1:3，把3个图放一起
gseaplot2(gsea_res,geneSetID = 1,title = "title",
          subplots = 1:3,
          base_size = 10)

gseaplot2(gsea_res, geneSetID = 1, subplots = 1)
gseaplot2(gsea_res, geneSetID = 1, subplots = 1:2)

#把entrezid变为symbol
gsea_res_symbol <- setReadable(gsea_res, "org.Hs.eg.db", "ENTREZID")

p <- gseaplot2(gsea_res_symbol,geneSetID = 1,
               title = gsea_res_symbol$Description[1])

p[[1]] <- p[[1]]+
  theme(title = element_text(color = "red"))
p

GSEA-5

2.2.6 展示多条通路

tmp <- as.data.frame(gsea_res_symbol)
colnames(tmp)
##  [1] "ID"              "Description"     "setSize"         "enrichmentScore" "NES"             "pvalue"          "p.adjust"       
##  [8] "qvalues"         "rank"            "leading_edge"    "core_enrichment"

head(tmp,3)
##                                                ID              Description setSize enrichmentScore      NES     pvalue  p.adjust   qvalues
## HALLMARK_APICAL_JUNCTION HALLMARK_APICAL_JUNCTION HALLMARK_APICAL_JUNCTION      36       0.3629219 1.835884 0.01019684 0.2330426 0.2026457
## HALLMARK_MYOGENESIS           HALLMARK_MYOGENESIS      HALLMARK_MYOGENESIS      70       0.2692732 1.620309 0.02868852 0.2330426 0.2026457
## HALLMARK_ADIPOGENESIS       HALLMARK_ADIPOGENESIS    HALLMARK_ADIPOGENESIS      18       0.4200174 1.691056 0.03039686 0.2330426 0.2026457
##                          rank                   leading_edge
## HALLMARK_APICAL_JUNCTION  885 tags=53%, list=31%, signal=37%
## HALLMARK_MYOGENESIS       416 tags=29%, list=15%, signal=25%
## HALLMARK_ADIPOGENESIS     671 tags=50%, list=24%, signal=38%
##                                                                                                                          core_enrichment
## HALLMARK_APICAL_JUNCTION CLDN6/ACTG2/ACTN2/FLNC/MYL9/CNTN1/CADM3/CADM2/SPEG/CLDN11/CLDN19/NEXN/NEGR1/NRXN2/SLIT2/KRT31/ACTA1/NFASC/ACTC1
## HALLMARK_MYOGENESIS        CASQ2/ACTN2/DES/MYH11/LDB3/PYGM/MYLK/CASQ1/TAGLN/HSPB8/FHL1/SGCA/CRYAB/PTGIS/DTNA/SPEG/SORBS1/IGF1/TPM2/GNAO1
## HALLMARK_ADIPOGENESIS                                                              ADIPOQ/FABP4/MYLK/HSPB8/CIDEA/SORBS1/SPARCL1/RETN/OMD

p <- gseaplot2(gsea_res, geneSetID = 1:6)
p

p[[1]] <- p[[1]] + scale_color_viridis_d(labels = c("aaa","bbb","ccc","ddd","eee","fff"))+
  geom_hline(yintercept = 0,color = "grey75", linewidth = 0.8, linetype = 2)+
  theme(legend.position = "top", legend.direction = "horizontal")
p[[2]] <- p[[2]] + scale_color_viridis_d()
p[[3]] <- p[[3]] + geom_hline(yintercept = 0, color = "steelblue", linewidth = 0.5, linetype = 2)
p

GSEA-6