#####clusterprofile-选用kegg数据库-GSEA函数--可
rm(list = ls())
library(dplyr)
library(data.table)
library(GSEABase)
library(clusterProfiler)
library(DOSE)
library(org.Mm.eg.db)
library(ggplot2)
library(stringr)
library(enrichplot)
library(msigdbr)
options(stringsAsFactors = F)
setwd("/data/wanglei_lab/zhangyu/2.zy_fdy/2.TBILJN_vs_MC/tbi-20221017/6.diff/1.GSEA/GSEAhs")
###预定基因集
mdb_c2 <- msigdbr(species = "Mus musculus", category = "C2")
mdb_c2 = mdb_c2 [grep("^KEGG",mdb_c2 $gs_name),] ##R语言中的grep函数可以在给定的字符串向量中搜索某个子字符串
mdb_c2$gs_name <- gsub('KEGG_','',mdb_c2$gs_name)   #去除前缀KEGG_
mdb_c2$gs_name <- tolower(mdb_c2$gs_name)  #将大写换为小写
mdb_c2$gs_name <- gsub('_',' ',mdb_c2$gs_name)  #将_转化为空格
mdb_c2 <- mdb_c2 %>% as.data.frame %>% dplyr::select(gs_name,entrez_gene,gene_symbol) %>% as.data.frame
head(mdb_c2)

##数据准备
diff_Ast<-read.csv("/data/wanglei_lab/zhangyu/2.zy_fdy/2.TBILJN_vs_MC/tbi-20221017/6.diff/1.GSEA/diff_Ast.csv", header=TRUE,row.names= 1)  #读取数据时将基因列读作行名
diff_Ast$SYMBOL = rownames(diff_Ast)  
df <- bitr(unique(diff_Ast$SYMBOL), fromType = "SYMBOL",toType = c( "ENTREZID"),OrgDb = org.Mm.eg.db) 
#去重
df.id <- dplyr::distinct(df.id,SYMBOL,.keep_all=TRUE)
merged=merge(diff_Ast,df,by='SYMBOL',all=F)
head(merged)
data_all_sort <- merged %>%  arrange(desc(avg_log2FC))  #排序
GSEA_input = data_all_sort$avg_log2FC  #把foldchange按照从大到小提取出来
names(GSEA_input) <- data_all_sort$ENTREZID  #给上面提取的foldchange加上对应上ENTREZID
head(GSEA_input) 

##GSEA富集
egmt <- GSEA(GSEA_input, TERM2GENE= mdb_c2,verbose=F,pvalueCutoff = 0.5)
kkx_C2 <- setReadable(egmt, 'org.Mm.eg.db', 'ENTREZID') %>% as.data.frame()
write.csv(kkx_C2,file = 'diff_Ast.GSEA.result20221101-333.csv',quote = F,row.names = T)
save(GSEA_KEGG,file = 'GSEA_KEGG.Rdata')
#作图
gseap1 <- gseaplot2(egmt,
                    'toll like receptor signaling pathway',#富集的ID编号
                    title = "toll like receptor signaling pathway",#标题
                    color = "red", #GSEA线条颜色
                    base_size = 20,#基础字体大小
                    rel_heights = c(1.5, 0.5, 1),#副图的相对高度
                    subplots = 1:3,   #要显示哪些副图 如subplots=c(1,3) #只要第一和第三个图
                    ES_geom = "line", #enrichment score用线还是用点"dot"
                    pvalue_table = T) #显示pvalue等信息
ggsave(gseap1, filename = 'GSEA_up_1.pdf', width =10, height =8)
gseaplot2(Hgsea,'HALLMARK_HYPOXIA',color = 'red',pvalue_table = T)