手把手教你用R做GSEA分析

手把手教你用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(clusterProf iler)library(biomaRt)library(enrichplot)
导入文件
导入的文件是两组间差异分析的结果，有基因名和logFC
genelist_input <- f read(f ile="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.f rame(gene_map),"基因转换.csv",s =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, pvalueCutof f=1)#GSEA 分析——KEGGKEGG_gseresult <- gseKEGG(geneList, nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutof f=1)#GSEA分析——
ReactomeGo_Reactomeresult <- gsePathway(geneList, nPerm = 1000, minGSSize = 10, maxGSSize = 1000, pvalueCutof f=1)#保存文件write.table (Go_gseresult, f ile ="Go_gseresult.csv", sep =",", s =TRUE)write.table (KEGG_gseresult, f ile ="KEGG_gseresult.csv", sep =",", s =TRUE)write.table
(Go_Reactomeresult, f ile ="Go_Reactomeresult.csv", sep =",", s =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)
本文的示例文件及全部代码。