#GSEA富集分析
setwd("D:/geo/GSEA/")

library(clusterProfiler)
library(org.Hs.eg.db) #人类的数据库，如果是鼠的则是 org.Mm.eg.db
library(stringr)

gene<-read.csv("TCGA_LUAD_DEG_Tumor_Contr.csv", header = T, row.names = 1) 
gene$SYMBOL = row.names(gene)

#开始ID转换
gene_entrizid=bitr(gene$SYMBOL,fromType="SYMBOL",toType="ENTREZID",OrgDb="org.Hs.eg.db") #如果是鼠的则需要改为org.Mn.eg.db


gene_df <- data.frame(logFC=gene$logFC, 
                      SYMBOL = gene$SYMBOL)
gene_df <- merge(gene_df,gene_entrizid,by="SYMBOL")

gene_df = gene_df[!duplicated(gene_df$ENTREZID) & gene_df$ENTREZID!="" & !is.na(gene_df$ENTREZID),]

geneList<-gene_df$logFC #第二列可以是folodchange，也可以是logFC
names(geneList)=gene_df$ENTREZID #使用转换好的ID
geneList=sort(geneList,decreasing = T) #从高到低排序


kegmt<-read.gmt("c2.all.v6.2.entrez.gmt") #读gmt文件
KEGG<-GSEA(geneList,TERM2GENE = kegmt) #GSEA分析

library(ggplot2)
pdf(file = "GSEA_dotplot.pdf")
dotplot(KEGG) #出点图 
dev.off()

pdf(file = "GSEA_dotplot_pvalue.pdf")
dotplot(KEGG,color="pvalue")  #按p值出点图 
dev.off()

pdf(file = "GSEA_dotplot_pvalue_sign.pdf",width=16,height=8)
dotplot(KEGG,split=".sign")+facet_grid(~.sign) #出点图，并且分面激活和抑制
dev.off()

pdf(file = "GSEA_dotplot_pvalue_sign_free.pdf",width=16,height=8)
dotplot(KEGG,split=".sign")+facet_wrap(~.sign,scales = "free") #换个显示方式
dev.off()


library(enrichplot)
#特定通路作图
pdf(file = "GSEA_plot2.pdf")
gseaplot2(KEGG,1,color="red",pvalue_table = T) # 按第一个做二维码图，并显示p值
dev.off()