clusterProfiler|GSEA富集分析及可视化

GSEA（Gene Set EnrichmentAnalysis），即基因集富集分析，无需设定阈值来区分上调下调基因，使用所有的基因进行分析。

GO 和 KEGG 可参考：R|clusterProfiler-富集分析

一 准备数据

1.1，加载R包，数据

library(org.Hs.eg.db)
library(clusterProfiler)
library(pathview)
library(enrichplot)

data <- read.csv("limmaOut.csv",header=TRUE)
head(data)

GSEA分析只需要两列信息，基因列和logFC（不同软件的差异分析这一列的名字会有差别）。

1.2 ID转化

基因名是symbol，将之转换为entrezid格式

gene <- data$SYMBOL
#开始ID转换，会有丢失
gene=bitr(gene,fromType="SYMBOL",toType="ENTREZID",OrgDb="org.Hs.eg.db")
#去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)

#合并data 和 entrezid

data_all <- data %>%
 inner_join(gene,by="SYMBOL")
dim(data_all)
head(data_all)

1.3 基因列表排序

将基因按照logFC进行从高到低排序，只需要基因列和logFC即可

data_all_sort <- data_all %>%
 arrange(desc(logFC))
head(data_all_sort)

geneList = data_all_sort$logFC #把foldchange按照从大到小提取出来
names(geneList) <- data_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(geneList)

#结果为排序的logGC，names为ENTREZID

345557     1308    55220     8788     6770     8745 2.425345 1.938050 1.831625 1.825617 1.786812 1.767371

二 GSEA分析

2.1 GSEA分析

需要gmt文件，http://www.gsea-msigdb.org/gsea/downloads.jsp路径下载，选择合适的

kegg_gmt <- read.gmt("c2.cp.biocarta.v7.4.entrez.gmt") #读gmt文件
gsea <- GSEA(geneList,
            TERM2GENE = kegg_gmt) #GSEA分析
head(gsea)

其中：

• ID ：信号通路

• Description ：信号通路的描述

• setSize ：富集到该信号通路的基因个数

• enrichmentScore ：富集分数，也就是ES

• NES ：标准化以后的ES，全称normalized enrichment score

• pvalue：富集的P值

• p.adjust ：校正后的P值

• qvalues ：FDR （false discovery rate）错误发现率

• rank ：排名

• core_enrichment：富集到该通路的基因列表。

可视化-点图

dotplot(gsea)

2.2 GSEA-GO分析

gse.GO <- gseGO(
 geneList, #geneList
 ont = "BP",  # 可选"BP"、"MF"和"CC"或"ALL"
 OrgDb = org.Hs.eg.db, #人 注释基因
 keyType = "ENTREZID",
 pvalueCutoff = 0.05,
 pAdjustMethod = "BH",#p值校正方法
)
head(gse.GO,2)

2.3 GSEA-KEGG分析

gse.KEGG <- gseKEGG(geneList,
                   organism = "hsa", # 人 hsa
                   pvalueCutoff = 1,
                   pAdjustMethod = "BH",) #具体参数在下面
head(gse.KEGG)
#为方便展示，最后一列暂不展示
head(gse.KEGG)[1:10]

三 GSEA可视化

使用gseaplot2函数进行可视化

3.1 简单可视化

gseaplot2(gse.KEGG,
         1) #展示第一个通路

可以进行一些调整以接近文献

1）修改GSEA线条颜色

2）添加P值的table

3）展示指定的通路

4）展示多个通路

5）只展示上两部分

3.2 展示指定通路

gseaplot2(gse.KEGG,
         title = "Olfactory transduction",  #设置title
         "hsa04740", #绘制hsa04740通路的结果
         color="red", #线条颜色
         base_size = 20, #基础字体的大小
         subplots = 1:2, #展示上2部分
         pvalue_table = T) # 显示p值

3.3 展示多个GSEA结果

A：使用数字的方式

gseaplot2(gse.KEGG,
         1:3, #绘制前3个
         pvalue_table = T) # 显示p值

B：使用向量指定通路

gseaplot2(gse.KEGG,
         c("hsa04740","hsa05310"), #指定通路向量
         pvalue_table = T) # 显示p值

C：点 形式

gseaplot2(gse.KEGG,

          1:5, #按照第一个作图

          ES_geom = "dot",

          base_size = 20,

          pvalue_table = T) 

应该和文献中的图很接近了，剩下的就是用自己的数据去尝试了。

四 参考资料： 

http://yulab-smu.top/clusterProfiler-book/