先加載R包

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

第一步：差異分析

首先是用limma包做差異分析的過(guò)程

#差異分析：
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)

這里由于做的時(shí)候是FPKM數(shù)據(jù)，所以參考了生信技能樹(shù)的教程：

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

第二步：獲取基因集

MSigDB可以方便地獲取GMT文件，讀進(jìn)來(lái)

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

第三步：進(jìn)行GSEA

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

GSEA輸入的geneList是排序后的gene名，因此這里需要對(duì)geneList進(jìn)行一個(gè)排序。實(shí)際上得到差異分析的結(jié)果后，只需要logFC和gene名就可以構(gòu)建geneList進(jìn)行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值太小時(shí)可能會(huì)報(bào)warning，不過(guò)不影響?？梢园裡ps這個(gè)參數(shù)設(shè)置成0。得到的結(jié)果是一個(gè)gseaResult對(duì)象

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

可視化

接下來(lái)可以對(duì)這個(gè)gseaResult對(duì)象進(jìn)行可視化。

> 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()
+ }

這里循環(huán)讀取感興趣的每一個(gè)基因集進(jìn)行作圖，并輸出為PDF。

GSEA

當(dāng)然也可以多個(gè)圖放在一起

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

比如這個(gè)樣子

gseaplot2