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)本文的示例文件及全部代碼
