手把手教你用R做GSEA分析

GSEA是非常常见的富集分析方式，以前我们做GSEA需要用依赖java的GSEA软件，那个时候准备分析的文件可能要花上很长时间，报错还不知道如何处理。现在我们来学习一下R语言进行GSEA分析。

加载R包

rm(list = ls())library(ReactomePA)library(tidyverse)library(data.table)library(org.Hs.eg.db)library(clusterProfiler)library(biomaRt)library(enrichplot)

导入文件

导入的文件是两组间差异分析的结果，有基因名和logFC

genelist_input <- fread(file="GENE.txt", header = T, sep='\t', data.table = F)genename <- as.character(genelist_input[,1]) #提取第一列基因名

基因名转换

将基因名转换成ENTREZID

gene_map <- select(org.Hs.eg.db, keys=genename, keytype="SYMBOL", columns=c("ENTREZID"))colnames(gene_map)[1]<-"Gene"write.csv(as.data.frame(gene_map),"基因转换.csv",row.names =F)#导出结果至默认路径下

将ENTREZID与logFC结合，并根据logFC的值降序排列

aaa<-inner_join(gene_map,genelist_input,by = "Gene")aaa<-aaa[,-1]aaa<-na.omit(aaa)aaa$logFC<-sort(aaa$logFC,decreasing = T)

GSEA文件准备

整理成GSEA分析的格式

geneList = aaa[,2]names(geneList) = as.character(aaa[,1])geneList

开始分析

接下来可以进行富集分析了，并保存结果文件

#GSEA分析——GOGo_gseresult <- gseGO(geneList, 'org.Hs.eg.db', keyType = "ENTREZID", ont="all", nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=1)#GSEA分析——KEGGKEGG_gseresult <- gseKEGG(geneList, nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=1)#GSEA分析——ReactomeGo_Reactomeresult <- gsePathway(geneList, nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutoff=1)#保存文件write.table (Go_gseresult, file ="Go_gseresult.csv", sep =",", row.names =TRUE)write.table (KEGG_gseresult, file ="KEGG_gseresult.csv", sep =",", row.names =TRUE)write.table (Go_Reactomeresult, file ="Go_Reactomeresult.csv", sep =",", row.names =TRUE)

可视化

用波浪图展示GO的前10个结果

#波浪图ridgeplot(Go_gseresult,10) #输出前十个结果

单个GSEA结果的展示

方法1

#富集曲线图类型1：gseaplot(Go_Reactomeresult,1,pvalue_table = TRUE) #输出第1个结果

方法2

#富集曲线图类型2：gseaplot2(Go_Reactomeresult,212,pvalue_table = TRUE)#输出第212个结果

方法3 同时展示多个结果

#gseaplot2还可以同时显示复数个功能组的富集曲线，并标记P值：gseaplot2(Go_Reactomeresult, 1:4, pvalue_table = TRUE)

本文的示例文件及全部代码