---

title: “GSEA plots script” output: html_notebook

title: "GSEA plots script"
library(tidyverse)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(gridExtra)
library(dplyr)
library(cowplot)
#import edgeR DEG results
 DEG.result <-read.csv("edgeR.csv",row.names = 1,header = TRUE)
DEG.result <- DEG.result %>% rownames_to_column("genename") 
dim(DEG.result)
genelist <- DEG.result$genename
## convert gene symbol to ENSEMBL and ENTREZID
genelist <- genelist %>%  bitr(fromType = "SYMBOL", toType = c("ENSEMBL", "ENTREZID"), OrgDb = org.Hs.eg.db)
head(genelist)
colnames(genelist)[1]<-"genename"
DEG.result <-genelist %>% inner_join(DEG.result, by='genename') %>%  dplyr::select(ENTREZID, logFC, everything())
head(DEG.result)
#convert format of genelist for GSEA
geneList<- DEG.result$logFC
names(geneList)<- as.character(DEG.result$ENTREZID)
geneList <-geneList[!duplicated(names(geneList))]
geneList<- sort(geneList,decreasing=TRUE)
head(geneList)
##LOAD Msigdbr
require(enrichplot)
library(msigdbr)
msigdbr_show_species()
#Show hallmark genes from C1-C7 category
m_t2g <- msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, entrez_gene)
head(m_t2g)
em <- GSEA(geneList,nPerm = 10000 ,TERM2GENE = m_t2g,pvalueCutoff = 0.05,verbose = FALSE)
#plot dot plot , where you can increase or decrease the category number
dotplot(em, showCategory=12, split=".sign")+ facet_grid(~.sign)
#plot GSEA
gseaplot2(em, geneSetID = 1, title = em$ID[1],color="red",pvalue_table=T)
gseaplot2(em, geneSetID = 1:3, title ="Activated Pathway",pvalue_table=T)
###
p<-list()
for (i in c(1:11)){
  p[[i]]<-gseaplot2(em, geneSetID = i, title=em$ID[i],color="red",pvalue_table=T)
  ggsave(sprintf("GSEAp%s.pdf", i), width = 8, height = 6, onefile = T)##批量保存
  p[[i]]
  dev.off()
}
p3 <-marrangeGrob(p,nrow = 2,ncol=2)
p3
#save plots
ggsave("GSEA.plots.pdf",width = 14,height = 8, plot = p3)
#use lapply funtion to get many plots together
GSEAplot<-function(x){
    gseaplot2(em, geneSetID =x,title=em$ID[x],color="red",subplots = ,
              pvalue_table=T)
}
GSEAplot(1)
##
p4<-lapply(c(1:nrow(em@result)), GSEAplot)
print(p4)
dev.off()