用clusterProfiler进行GSEA分析

****今天来关注如何使用clusterProfiler系列包来分析GSEA**

导入包，导入数据

library(clusterProfiler)
library(enrichplot)
library(ReactomePA)
library(data.table)
genelist_input <- fread(file="input.txt", header = T, sep='\t', data.table = F)
head(genelist_input)  NCBI gene ID           R1            1  0.051631112            2  0.031268333            9  0.120690144           10  0.026080255           12 -0.156296536           13 -0.06157593

今天比昨天多导入了一个ReactomePA包，这个包也是Y叔开发用来分析Reactome数据库中基因功能组的，包比较大，五百多M，装的时候要等等。注意今天使用的数据实例是以NCBI的基因数字ID作为基因名。

待分析数据格式处理

如果直接拿来分析就会报错：

Error in `[.data.frame`(x, order(x, na.last = na.last, decreasing = decreasing)) :   undefined columns selected

研究发现gsea读取的数据不包括表头，即首行列名，所以要处理一下（如此处理后数据也变成降序排列）：

geneList = genelist_input[,2]names(geneList) = as.character(genelist_input[,1])geneList = sort(geneList, decreasing = TRUE)

GSEA分析——GO

#Go_gseresult <- gseGO(geneList, 'org.Hs.eg.db', keyType = "ENTREZID", ont="BP", nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=0.05)

GSEA分析——KEGG

KEGG_gseresult <- gseKEGG(geneList, nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=0.05)

GSEA分析——Reactome

Go_Reactomeresult <- gsePathway(geneList, nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=0.05)

这三种分析的结果都可以输出文本查看。有时候会出现这种情况：没有富集到任何结果。研究后发现，pvalueCutoff界定值应该背锅，虽然这里写的是P值，但这个包里实际是以FDR作为界定。经常用FDR的同学会遇到这种情况：P值很多很显著的，但是FDR算出来都是一个值，时而还都是大于0.05的相同值。所以解决办法就是pvalueCutoff设置等于1，要是不设定也不行，只要设定成1，就会输出所有结果。

注意一下这里的参数：nPerm是置换次数；minGSSize和maxGSSize是指用来分析的基因功能组中基因数目的上下界值，这里设定是10-1000，那么小于10或者大于1000基因数的功能组将不会被本次分析所使用。

可视化

以上面Reactome为例，首先查看波浪图：

ridgeplot(Go_Reactomeresult, 10) #输出前十个结果

图片

富集曲线图类型1：

gseaplot(Go_Reactomeresult, 1) #输出第一个结果

图片

富集曲线图类型2：

gseaplot2(Go_Reactomeresult, 1)

图片

还是gseaplot2看着舒服，这里提示一下，要用最新版enrichplot包才能使用gseaplot2作图。gseaplot2还可以同时显示复数个功能组的富集曲线，并标记P值：

gseaplot2(Go_Reactomeresult, 1:3, pvalue_table = TRUE)

图片

GO和KEGG的GSEA结果可视化同理，下次再来介绍如何使用clusterProfiler系列包分析自定义基因功能组~~~

参考文献

Yu G, He Q (2016). “ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization.” Molecular BioSystems, 12(12), 477-479.
Yu G, Wang L, Han Y, He Q (2012). “clusterProfiler: an R package for comparing biological themes among gene clusters.” OMICS: A Journal of Integrative Biology, 16(5), 284-287.
Guangchuang Yu (2018). enrichplot: Visualization of Functional
Enrichment Result. R package version 1.2.0.
https://github.com/GuangchuangYu/enrichplot
原文：https://mp.weixin.qq.com/s/szmquOt77e6Cad4Glhi-aA