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)
应该和文献中的图很接近了，剩下的就是用自己的数据去尝试了。
