网页上关于GSEA富集分析的教程很多，对自己这个生信小白来说比较简单通俗易懂的当数生信宝典和生信技能树的教程了。感谢大佬们的付出。下面是自己整理的个人学习笔记，可能有理解不对或者断章取义的地方，有问题还是参考大佬们的教程。
什么是GSEA?
GSEA:gene set enrichment analysis即基因集富集分析，是常见富集分析之一，类似GO/KEGG.
先给定一个排序的基因表L和一个预先定义的基因集S (比如编码某个代谢通路的产物的基因, 基因组上物理位置相近的基因，或同一GO注释下的基因)，GSEA的目的是判断S里面的成员s在L里面是随机分布还是主要聚集在L的顶部或底部。这些基因排序的依据是其在不同表型状态下的表达差异，若研究的基因集S的成员显著聚集在L的顶部或底部，则说明此基因集成员对表型的差异有贡献，也是我们关注的基因集。
来自:生信宝典富集分析 - 界面操作
原理：
测序或者芯片分析一般是根据研究目的设置两个组，然后通过统计学检验方法（t-test、或者方差分析等）寻找差异表达的基因，然后根据差异表达基因推断相关的生物学机制。
差异表达基因需满足以下条件：
•统计学上有意义：, p ≤0.01
•控制假阳性率：FDR，一般采用BH法进行校正，FDR ≤ 0.1 or 0.2
•FoldChange足够大： |FC| ≥ 1.5或者2
•探针平均表达水平较高组大部分表达值都高于背景：＞80%但如果组内样本仅为3个时必须3个都满足；如不满足这些条件，不宜进行统计学检验分析
•另外，根据研究目的和测序平台等，筛选差异表达基因的条件可以调整。
但是进行差异表达分析，可能遇到的问题是：找到了一堆DEGs，却无法得到统一的生物学推论，因为DEG的标准是人为设置较为随意，又或者很少DEG满足人为设置的标准。
关于在DEG的列表中进行GSEA需注意以下几点:
1.检验DEGs的成员是否在给定pathway和其他基因集中显著高表达？
基因满足组间比较中显著上调的标准；
基因与某一临床变量高度正（负）相关
基因在组间比较中有较大的差异倍数
2.进行分析的工具：DAVID/GoMinor/GenMAPP/Onto-Express/Gostat/ GATHER等
3.分开提交上调和下调的基因列表
4.基因集富集分析只能得到验证过的在数据库中的基因集相关的结果。
举个例子：如果基因芯片共检测10，000个基因，而其中200个待测基因在已知基因集S中，测到了50个差异表达基因。
那么一般来说随机至少有1个基因落在基因集S上（200/10000*50);但最后结果得到了6个差异表达基因落在基因集S中，根据超几何分布推测，该事件发生的p值为0.00045，即随机抽取50个基因，有相同或者更多基因落在基因集S的概率。
基因集富集方法：
•这些方法为利用选定的度量标准计算基因集中每个基因提供了统计学公式，增加了统计学效力如
1.组A和组B比较的T-评分
2.组Avs组B的FC值
3.基因表达与某一临床变量的pearson相关系数
•数据集中所有的基因都被运用到富集分析中
•可被运用到多个基因集中
•值得注意的是，结果取决于选取的geneset；另外存在多次检验导致的错误；建议分开进行不同种类的基因集分析。
GSEA步骤分为三步：
1.计算富集分数（ES）： ES反映基因集成员s在排序列表L的两端富集的程度。计算方式是，从基因集L的第一个基因开始，计算一个累计统计值。当遇到一个落在s里面的基因，则增加统计值。遇到一个不在s里面的基因，则降低统计值。每一步统计值增加或减少的幅度与基因的表达变化程度（更严格的是与基因和表型的关联度）是相关的。富集得分ES最后定义为最大的峰值。正值ES表示基因集在列表的顶部富集，负值ES表示基因集在列表的底部富集。
2.估计富集分数的显著性水平：通过基于表型而不改变基因之间关系的排列检验 (permutation test)计算观察到的富集得分(ES)出现的可能性。若样品量少，也可基于基因集做排列检验 (permutation test)，计算p-value。
3.矫正多重假设检验：首先对每个基因子集s计算得到的ES根据基因集的大小进行标准化得到Normalized Enrichment Score (NES)。随后针对NES计算假阳性率。（计算NES也有另外一种方法，是计算出的ES除以排列检验得到的所有ES的平均值）
经过GSEA富集分析后可以得到对genesets贡献最大的部分基因，这部分基因被称为：Leading-edge subset：对富集得分贡献最大的基因成员。根据这些基因，后续进行Leading-edge analysis
GSEA和GO/KEGG富集分析的差异：
1.GO/KEGG分析强调差异基因，是针对统计学存在显著差异的一部分基因的分析 ，容易遗漏一些统计学上无显著差异但具有生物学意义的基因；而GSEA是从全体基因的表达矩阵中找出具有协同差异 (concordant differences)的基因集，故能兼顾差异较小的基因。
2.GO富集是定性的分析，GSEA考虑到了表达或其它度量水平的值的影响。
3.对于时间序列数据或样品有定量属性时，GSEA的优势会更明显，不需要每个分组分别进行富集，直接对整体进行处理。GO/KEGG更适合两组间的分析。
GSEA软件
软件安装：教程：http://software.broadinstitute.org/gsea/doc/desktop_tutorial.jsp ，需要java环境
输入数据一般格式：http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats
差异表达数据： *.gct
表型、分组 ：*.cls
基因集： *.gmt
芯片注释：*.chip
输出格式：点击success就能进入报告主页，里面的链接可以进入任意一个分报告。一般是基于MSigDB数据库的富集结果。解读参照网页：：http://www.broadinstitute.org/cancer/software/gsea/wiki/index.php/Main_Page
具体操作：GSEA富集分析 - 界面操作
R语言进行GSEA分析
用到的R包依旧是是做富集分析的神器-Y叔的clusterProfiler
同样需要排序基因表L和预先定义的基因集S这两个东西。
先看排序的基因表L:这里需要构建一个genelist，可以是来自自己测序数据差异分析的结果或者是GEO数据集，anyway，这个genelist由两列构成，第一列表示差异表达的基因ID（基因ID不能重复的,形式同样为entrezID），第二列为基因对应的表达量或者是FC等数值型向量，注意按照数值从高到低排列（为什么要由高到低因为要制作排序 的基因表L）
gene <- names(geneList) ##提取排序后基因名
接下来需要基因集S，可以直接读取自己提前准备好的gmt文件或者clusterProfiler的自带基因集文件或者从R包msigdf中提取
##gmtfile <- system.file("extdata", "c5.cc.v5.0.entrez.gmt", package="clusterProfiler")
S<-read.gmt(gmtfile） ##读取gmt文件得到基因集S
好了，准备好了排序基因集L----gene，和基因集S ------S，就可以进行超几何分布检验和GSEA富集分析
egmt <- enricher(gene, TERM2GENE=S)##超几何分布检验
head(egmt) ##查看GSEA分析结果
write.csv(egmt,file ="egmt.csv") ##保存超几何分布检验
egmt2 <-GSEA(geneList,TERM2GENE=S,pvalueCutoff=0.5,#verbose=FALSE)
head(egmt2)
gseaplot2(egmt2,geneSetID=1,title=egmt2$Description[1],pvalue_table=T)
barplot(egmt）##可视化，展示富集的基因集
另外，也对于GO和KEGG的GSEA富集分析，可以直接用gseGO()、gseKEGG()等。
还可以对leading-edge subset进行可视化分析： dotplot() cnetplot()
具体操作参考Y叔的dotplot for GSEA以及enrichplot: 让你们对clusterProfiler系列包无法自拔
同时前面提到，可以通过安装R包msigdf可以免去从MSigDB下载gmt文件的步骤。但需注意msigdf的更新情况并不和数据库一致。参考Y叔的msigdf + clusterProfiler全方位支持MSigDb一文。
devtools::install_github("ToledoEM/msigdf")
library(msigdf)
library(dplyr)
c2 <- msigdf.human %>% filter(category_code == "c2") %>% select(geneset, symbol) %>% as.data.frame
这里得到基因集c2，注意这个包中的基因ID都是SYMBOLID，需要转换为ID，再用clusterProfiler进行GSEA分析。
再次感到clusterProfiler的强大！！！
如何自己制作geneset做GSEA分析
如果用GSEA软件进行GSEA分析，需要GCT,CLS和GMT文件。其中GMT文件可以在数据库MSigDB中找到相应的geneset。但是不是所有想要的geneset都能在MSigDB找到，这时就需要我们自己制作geneset的GMT格式文件。
首先来看看一般GSEA用到的GMT格式文件长什么样子？
可以看到，表格每一行代表一个geneset，第一列是geneset的名称，不能重复命名；第二列是geneset的描述，没有是可以用“na”填充；之后每一列是该geneset包含的基因名，在Excel编辑时应当注意基因名是否正确。每个geneset包含的基因数可以不同。
制作前要准备好感兴趣的geneset的genelist，可以通过查阅文献和相关数据库如KEGG等等得到，基因名必须是以hugo规定的标准symbolID。
生信技能树给出的示范：
##R语言中将kegg的文件kegg2gene.txt读入path2gene_file文件中
tmp=read.table(path2gene_file,sep="\t",colClasses=c('character'))
#tmp=toTable(org.Hs.egPATH)
## 读入的文件，第一列是kegg的编号，第2列是基因的entrezID
GeneID2kegg_list<<- tapply(tmp[,1],as.factor(tmp[,2]),function(x)x)
kegg2GeneID_list<<-tapply(tmp[,2],as.factor(tmp[,1]),function(x)x)
##最后需要将基因的entrezID转换为symbolID，得到kegg2symbol_list
具体操作可以参照gene的各种ID转换终结者-bioconductor系列包（来自生信技能树）
最后将genelist保存为GMT格式文件
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()}
这样就可以把想要的geneset转换为GSEA分析的GMT文件了。制作重点还是在于弄清输入的GMT文件的格式。
参考：GSEA的分析汇总
