先加载R包

> library(GSEABase)
> library(limma) 
> library(clusterProfiler)
> library(enrichplot)

第一步：差异分析

首先是用limma包做差异分析的过程

#差异分析：
dat <- exprSet
group_list <- c(rep("Control",3),
                rep("Combine",3))
group_list <- factor(group_list,ordered = F)
design <- model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
options(digits = 4)
topTable(fit,coef=2,adjust='BH')
deg=topTable(fit,coef=2,adjust='BH',number = Inf)
head(deg)

这里由于做的时候是FPKM数据，所以参考了生信技能树的教程：

https://mp.weixin.qq.com/s/BEvFqu-BNA9lWFNaE9Jkdw

第二步：获取基因集

MSigDB可以方便地获取GMT文件，读进来

> geneset <- read.gmt("./h.all.v7.5.1.symbols.gmt")

第三步：进行GSEA

> geneList <- deg$logFC
> names(geneList) <- toupper(rownames(deg))
> geneList <- sort(geneList,decreasing = T)

GSEA输入的geneList是排序后的gene名，因此这里需要对geneList进行一个排序。实际上得到差异分析的结果后，只需要logFC和gene名就可以构建geneList进行GSEA富集分析了。

> gsea_results <- GSEA(
+   geneList = geneList,
+   TERM2GENE = geneset,
+   verbose = F,
+   eps=1e-10
+ )
Warning messages:
1: In serialize(data, node$con) :
  'package:stats' may not be available when loading
2: In serialize(data, node$con) :
  'package:stats' may not be available when loading
3: In fgseaMultilevel(...) :
  For some pathways, in reality P-values are less than 1e-10. You can set the `eps` argument to zero for better estimation.

P值太小时可能会报warning，不过不影响。可以把eps这个参数设置成0。得到的结果是一个gseaResult对象

> class(gsea_results)
[1] "gseaResult"
attr(,"package")
[1] "DOSE"

可视化

接下来可以对这个gseaResult对象进行可视化。

> for (i in 1:length(list)) {
+   p <- gseaplot2(x=gsea_results,geneSetID=paste("HALLMARK_",list[i],sep="")) 
+   d <- paste("./",list[i],".pdf",sep="")
+   pdf(file=d,
+       family = "Times",width=10,height = 6)
+   print(p)
+   dev.off()
+ }

这里循环读取感兴趣的每一个基因集进行作图，并输出为PDF。

GSEA

当然也可以多个图放在一起

> p <- gseaplot2(gsea_results,
+           gsea_results@result[["ID"]][1:3],
+           pvalue_table = TRUE)

比如这个样子

gseaplot2