基本概念

Gene Set Enrichment Analysis (GSEA)是一种的计算方法，用于评估在一个预设订的基因集中，两个生物学数据集间，是否存在统计学上显著且一致的差异。

Gene Set Enrichment Analysis (GSEA) is a computational method that determines whether an a priori defined set of genes shows statistically
significant, concordant differences between two biological states.

GSEA富集分析和传统富集分析的比较

传统富集分析特点：

传统的富集分析的主要分析过程包括：采取固定阈值的筛选差异基因（p值，差异倍数），然后对其进行功能/通路富集分析（GO/KEGG），简单来说就是一种归类的方法。这种方法存在一个问题，一个通路的各个基因的表达量差异很大，生物调控也是递进关系，上游基因的微小变化可能就会导致下游基因的明显变化，如果直接一刀切用定阈值筛选，会损失很多基因信息，导致后续的富集不显著（基因数量太少做这种富集分析没有生物学意义）。

GSEA富集分析的特点

GSEA的分析流程主要包括：计算指标、排序、标注基因集、计算基因集得分、置换检验。从流程就可以看出，GSEA富集分析不需要用固定阈值去过滤基因，是基于所有基因分析的一种方法，规避了传统富集分析方法的短板。

因此，简单来说GSEA分析适用于，传统分析方法筛选后样本过少的数据集。

工具

GSEA软件下载：https://www.gsea-msigdb.org/gsea/index.jsp 直接选择你操作系统对应的版本，下载安装即可。

数据准备

基因表达量表（.gct）

基因表达量表

横坐标为样本，纵坐标为基因，陈列基因表达量信息；第一行#1.2为信息注释，不需要改动；第二行两个数字分别为：基因数量和样本数量（如实填写）；.gct文件，直接把表格用excel保存为制表符分隔文件（.txt），然后直接改后缀为gct即可。

样本分类关系表（.cls）

样本分类关系表

第一行三个数字分别为：样品总数，分组数，恒定数字1；第二行为样品分组名称；第三行为样品分组信息。需要注意：第二行和第三行需严格和.gct文件对应。

基因集合（.gmt）

基因集文件

第一列为基因集编号，第二列为基因集名称，第三列以后为集合中所包含的所有基因名称，注意基因名称需要与.gct文件中的名称一致。

基因集文件的准备相对比较复杂，我用到的是KEGG的通路基因集，也可以用GO的的基因集，甚至可以自己创造一个，只需要满足上图所示的文件格式。下面介绍一下如何从KEGG网站上下载某物种（小鼠）的基因集信息。

首先，在KEGG数据库中下载小鼠的KEGG通路数据库文件

找到小鼠的KEGG库
找到小鼠的KEGG库2.0
找到小鼠的KEGG库3.0
下载.keg文件
.keg文件

大家可以看到，这个.keg文件和我们需要的.gmt文件还有不小的差异，那么，我们需要按如下方法，先在Linux下解析keg文件，然后用R制作.gmt文件：

#解析信息
perl -alne '{if(/^C/){/PATH:mmu(\d+)/;$kegg=$1}else{print "$kegg\t$F[1]" if /^D/ and $kegg;}}' mmu00001.keg > mmu00001_kegg2gene.txt
grep '^C' mmu00001.keg | grep "mmu" | sed 's/^C    //g' | sed 's/\([0-9]*\)\s/\1\t/' > mmu00001_kegg2description.txt

#转换ID同时制作.gmt文件
#安装转换ID需要用到的各个R包
#安装data.table
install.packages("data.table")
#安装org.Xx.eg.db系列包中小鼠资源包
BiocManager::install("org.Mm.eg.db")
#安装clusterProfiler包，安装时需要“科学上网”
BiocManager::install("clusterProfiler")

#调用各个R包
library(data.table)
library(org.Mm.eg.db)
library(clusterProfiler)

#读入数据
path2gene_file<-"mmu00001_kegg2gene.txt"
tmp=read.table(path2gene_file,sep="\t",colClasses=c('character'))
dt <- tmp
dt$geneid <- rownames(dt)

# transform id  
map_dt <- bitr(dt$V2, fromType = "ENTREZID",toType = c( "SYMBOL"),OrgDb = org.Mm.eg.db)
dt_merge <- merge(map_dt,dt, by.y = "V2", by.x = "ENTREZID")

# 可以把转换ID后的表输出
#write.table(dt_merge, file = "example42.txt", sep = "\t", col.names = T, row.names = T, quote = F)

# first column is kegg ID, second column is entrez ID
GeneID2kegg_list<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
kegg2symbol_list<- tapply(dt_merge$SYMBOL,as.factor(dt_merge$V1),function(x) x)

#输出数据，但描述列为NA
write.gmt <- function(geneSet=kegg2symbol_list,gmt_file='kegg2symbol.gmt'){
  sink( gmt_file )
  for (i in 1:length(geneSet)){
    cat(names(geneSet)[i])
    cat('\tNA\t')
    cat(paste(geneSet[[i]],collapse = '\t'))
    cat('\n')
  }
  sink()
}

write.gmt()

GSEA分析

导入数据

打开GSEA软件，点击左侧菜单栏的Load data选项，将数据按如下图所示方法导入。

导入数据

设置参数

按照下图标记的批注设置合适的参数。

Required fields
Basic fields
Advanced fields

运行GSEA

设置好参数后点击Run运行即可。

点击RUN运行

GSEA结果展示

GSEA的所有结果都保存在结果文件夹中，完整报告为index.html文件，双击打开即可查看。

index.html文件总览信息

172/230表示在这一分组中，一共有230个功能基因集，其中172个上升
97个基因集的FDE小于25%
76个基因的名义P值小于1%
87个基因的名义P值小于5%
点击snapshot可以看富集结果
点击enrichment result in html 可以查看所有的富集分析结果，进去之后可以查看每个Enrichment plot的参数
点击enrichment result in excel就可以直接下载附带结果的excel

我们重点关注enrichment result in html中的内容，点开后会出现如下图所示的表格：

Enrichment result in html

图中：
SIZE：表示基因集里的基因数量
ES（enrichment score）：富集分数（重点）
NES（normalized enrichment score）:表示校正后的富集分数（重点）
NOM p-val (nominal p value ): 名义P值（重点）
FDR q-val（false discovery rate）：错误发现率（重点）
FWER p-val：用bonferonni校正后的P值
RANK AT AMX：ES值对应的通路基因排名
Leading-edge subset：对富集贡献最大的基因成员,即领头亚集，用于定义Leading-edge subset的参数有：Tags，List，Signal。

点击每个通路的Details即可看到如下图所示的典型GSEA数据展示图：

GSEA数据展示图

绿色曲线就是gene set里面对应的每个基因的enrichment score值（ES），开始时为零，从左到右每遇到一个基因就计算出一个ES值，连成一条绿线。当ES值大于0时，表示某一功能基因富集在排序序列的前端，若为小于0时，则某一功能基因富集在排序序列的后端，ES值越高说明这些基因在通路中有富集，非散在分布。
中间条形码似的黑线是gene set里面的基因在背景基因里的位置，每条竖线代表该通路下的基因。
蝴蝶图：按相关系数对所有基因排序，大于0正相关，小于0位负相关，越靠近两级相关性越弱。

对于结果的分析，通常认为|NES|>1，NOM p-val<0.05，FDR q-val<0.25的通路下的基因集合是有意义的；NES的绝对值越大，FDR值就越小，说明分析的结果可信度越高。NOM p-val是针对某一功能基因集得到的ES值的统计显著性，P值越小，说明基因的富集性越好，但P值很小时，FDR值也可能很大，这说明和其他功能基因子相比较，它的富集并不是很显著，原因可能是数据样本量较少、杂交信号微弱或者是选择的功能基因子集并未很好得反映样本的物理学意义。

参考：
如何实现GSEA-基因富集分析？