rm(list=ls())
#Gene ID conversion
library(data.table)
library(org.Hs.eg.db)
library(clusterProfiler)
library(biomaRt)
library(enrichplot)
library(tidyverse)

genelist_input <- fread(file="interGene.txt", header = F, sep='\t', data.table = F)
genename <- as.character(genelist_input[,1]) #Extract the first column of gene names
#x - Genome Annotation R Package
#keys - List of genes to be converted
#keytype - Gene name type
#columns - Gene name type data you want to return, such as NCBI gene ID, use ENTREZID
gene_map <- biomaRt::select(org.Hs.eg.db, keys=genename, keytype="SYMBOL", columns=c("ENTREZID"))
gene_map
write.csv(as.data.frame(gene_map),"Gene conversion.csv",row.names =F)#Export results to the default path
genelist_input<-gene_map[,2]

head(genelist_input)
genelist_input<-na.omit(genelist_input)
#GO analyze
Go_result_BP <- enrichGO(genelist_input, 'org.Hs.eg.db', ont="BP", pvalueCutoff=1) #The gene ID type is the ID format of ENSEMBL, the BP functional group is selected, and the P value is 0.05 as the limit

goplot(Go_result_BP, showCategory=5) #GO Topology

p1<-dotplot(Go_result_BP, showCategory=20) #Bubble chart showing the top twenty
p1<-p1 + scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(p1,filename ='GO_BP Bubble chart.pdf',width = 7.23,height = 8)
ggsave(p1,filename ='GO_BP Bubble chart.TIFF',width = 7.23,height = 8)

barplot(Go_result_BP, showCategory=20) #Bar chart showing the top twenty
y=as.data.frame(Go_result_BP)
y$geneID=as.character(sapply(y$geneID,function(x)paste(gene_map$SYMBOL[match(strsplit(x,"/")[[1]],as.character(gene_map$ENTREZID))],collapse="/")))
write.csv(y,"GO-BP.csv",row.names =F)#Export the results to the default path.
save(y,file = 'GO-BP.RDATA')

Go_result_CC <- enrichGO(genelist_input, 'org.Hs.eg.db', ont="CC", pvalueCutoff=10) #The gene ID type is the ID format of ENSEMBL, the CC functional group is selected, and the P value is 0.05 as the limit

goplot(Go_result_CC, showCategory=5) #GO Topology

p1<-dotplot(Go_result_CC, showCategory=20) #Bubble chart showing the top twenty
p1<-p1 + scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(p1,filename ='GO_CC Bubble chart.pdf',width = 7.23,height = 8)
ggsave(p1,filename ='GO_CC Bubble chart.TIFF',width = 7.23,height = 8)

barplot(Go_result_CC, showCategory=20) #Bar chart showing the top twenty

y=as.data.frame(Go_result_CC)
y$geneID=as.character(sapply(y$geneID,function(x)paste(gene_map$SYMBOL[match(strsplit(x,"/")[[1]],as.character(gene_map$ENTREZID))],collapse="/")))
write.csv(y,"GO-CC.csv",row.names =F)#Export the results to the default path
save(y,file = 'GO-CC.RDATA')

Go_result_MF <- enrichGO(genelist_input, 'org.Hs.eg.db',ont="MF", pvalueCutoff=1000) #The gene ID type is the ID format of ENSEMBL, the MF functional group is selected, and the P value is 0.05 as the limit

goplot(Go_result_MF, showCategory=5) #GO Topology

p1<-dotplot(Go_result_MF, showCategory=20) #Bubble chart showing the top twenty
p1<-p1 + scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(p1,filename ='GO_MF Bubble chart.pdf',width = 7.23,height = 8)
ggsave(p1,filename ='GO_MF Bubble chart.TIFF',width = 7.23,height = 8)

barplot(Go_result_MF, showCategory=20) #Bar chart showing the top twenty
y=as.data.frame(Go_result_MF)
y$geneID=as.character(sapply(y$geneID,function(x)paste(gene_map$SYMBOL[match(strsplit(x,"/")[[1]],as.character(gene_map$ENTREZID))],collapse="/")))

write.csv(y,"GO-MF.csv",row.names =F)#Export the results to the default path
save(y,file = 'GO-MF.RDATA')


## GO analysis #BP MF CC analysis together
go <- enrichGO(genelist_input, OrgDb = "org.Hs.eg.db", ont="all")
library(ggplot2)
p <- dotplot(go, split="ONTOLOGY") +facet_grid(ONTOLOGY~., scale="free")
p <-p + scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(p,filename ='GO Three-in-one bubble chart.pdf',width = 7.23,height = 8)
ggsave(p,filename ='GO Three-in-one bubble chart.TIFF',width = 7.23,height = 8)
write.csv(as.data.frame(go),"GO.csv",row.names =F)#Export the results to the default path


#KEGG analyze

KEGG_result <- enrichKEGG(genelist_input, keyType = "kegg",pvalueCutoff=1,qvalueCutoff=1,pAdjustMethod = "BH", minGSSize = 5, maxGSSize = 500,organism = "hsa", use_internal_data=T)  #KEGG enrichment analysis


barplot(KEGG_result, showCategory=20)#Draw a bar chart

p1<-dotplot(KEGG_result, showCategory=20) #Bubble chart showing the top twenty
p1<-p1 + scale_y_discrete(labels=function(x) str_wrap(x, width=50))
ggsave(p1,filename ='KEGG Bubble chart.pdf',width = 7.23,height = 8)
ggsave(p1,filename ='KEGG Bubble chart.TIFF',width = 7.23,height = 8)


#Circle Diagram
pdf(file="KEGG_circos.pdf",width = 10,height = 7)
kkx=setReadable(KEGG_result, 'org.Hs.eg.db', 'ENTREZID')
cnetplot(kkx, showCategory = 5, circular = TRUE, colorEdge = TRUE,node_label="all")
dev.off()
x=as.data.frame(KEGG_result)
x$geneID=as.character(sapply(x$geneID,function(x)paste(gene_map$SYMBOL[match(strsplit(x,"/")[[1]],as.character(gene_map$ENTREZID))],collapse="/")))


write.csv(as.data.frame(x),"KEGG.csv",row.names =F)#Export the results to the default path
save(x,file = 'KEGG.RDATA')


#Sort by P value, not padjust
#First, I retrieved the data frame from the enrichment object x, which is an S4 object, using the @ symbol to retrieve it. It has 293 rows, which are unfiltered data.
y =KEGG_result@result

#There is a problem with our horizontal axis because GeneRatio here is a string, and we now want to convert it into a numeric value.
## Remove the numbers before and after the semicolon and convert them into numerical values.
forward <- as.numeric(sub("/\\d+$", "", y$GeneRatio))
backward <- as.numeric(sub("^\\d+/", "", y$GeneRatio))
## Add a column of GeneRatio representing numerical values
y$GeneRatio = forward/backward

showCategory =20
#Set a font size, which can be adjusted later.
font.size =12

library(ggplot2)
library(forcats)
library(dplyr)
#Reproduce bubble chart
y %>% 
  ## Install p-value sorting, selecting a given number of rows
  arrange(pvalue) %>% 
  slice(1:showCategory) %>% 
  ## Start ggplot2 plotting, where fct_reorder adjusts the order of factor levels
  ggplot(aes(GeneRatio,forcats::fct_reorder(Description,Count)))+ 
  ## Draw a point graph
  geom_point(aes(color=pvalue, size = Count)) +
  ## Adjust the color, guide_colorbar adjusts the direction of the color map
  scale_color_continuous(low="red", high="blue", guide=guide_colorbar(reverse=TRUE))+
  ## Adjust the size of the bubble
  scale_size_continuous(range=c(3, 8))+
  ## If you use ylab("") or there is a blank on the left
  labs(y=NULL) +
  ## If there is no such sentence, the top will reach the top
  ggtitle("")+
  ## Set the theme
  theme_bw() +
  theme(axis.text.x = element_text(colour = "black",
                                   size = font.size, vjust =1 ),
        axis.text.y = element_text(colour = "black",
                                   size = font.size, hjust =1 ),
        axis.title = element_text(margin=margin(10, 5, 0, 0),
                                  color = "black",size = font.size),
        axis.title.y = element_text(angle=90))


ggsave(filename ='KEGG Bubble chart2.pdf',width = 10,height = 8)
ggsave(filename ='KEGG Bubble chart2.TIFF',width = 10,height = 8)


#Reproduce the bar chart
## Set the number of displays
showCategory =20
## Set the font size
font.size =12
y %>% 
  ## Install p-value sorting, selecting a given number of rows
  arrange(pvalue) %>% 
  slice(1:showCategory) %>% 
  ## Start ggplot2 plotting, where fct_reorder adjusts the order of factor levels
  ggplot(aes(x=forcats::fct_reorder(Description,pvalue,.desc = T),y=Count,fill=pvalue))+ 
  ## Draw a bar chart
  geom_bar(stat="identity")+
  coord_flip()+
  ## Adjust the color, guide_colorbar adjusts the direction of the color map
  scale_fill_continuous(low="red", high="blue", guide=guide_colorbar(reverse=TRUE))+
  ## If you use ylab("") or there is a blank on the left
  labs(x=NULL,y=NULL) +
  ## If there is no such sentence, the top will reach the top
  ggtitle("")+
  ## Set the theme
  theme_bw() +
  theme(axis.text.x = element_text(colour = "black",
                                   size = font.size, vjust =1 ),
        axis.text.y = element_text(colour = "black",
                                   size = font.size, hjust =1 ),
        axis.title = element_text(margin=margin(10, 5, 0, 0),
                                  color = "black",size = font.size),
        axis.title.y = element_text(angle=90))


#GO The results are drawn as ordered bar graphs
library(ggpubr)
library(tidyverse)
GODATA<-go@result
GODATA$"-log10(Pvalue)"<-  -log10(GODATA$pvalue)
GODATA$yyy<-  -log10(GODATA$pvalue)
colnames(GODATA)


write.csv(GODATA,"GODATA.csv",row.names =F)#Export the results to the default path

#Select the first ten BP MF CCs respectively
aaa<-filter(GODATA,ONTOLOGY=='BP')
aaaa<-aaa[1:10,]

bbb<-filter(GODATA,ONTOLOGY=='CC')
bbbb<-bbb[1:10,]

ccc<-filter(GODATA,ONTOLOGY=='MF')
cccc<-ccc[1:10,]


drawdata<-rbind(aaaa,bbbb,cccc)


ggbarplot(drawdata, x = "Description", y = 'yyy',
          fill = "ONTOLOGY",               # change fill color by cyl
          color = "white",            # Set bar border colors to white
          palette = c('#5FB404','#01DFD7','#C238E5'),            # jco journal color palett. see ?ggpar
          sort.val = "desc",          # Sort the value in dscending order
          sort.by.groups = FALSE,     # Don't sort inside each group
          x.text.angle = 90           # Rotate vertically x axis texts
)
ggbarplot(drawdata, x = "Description", y = "yyy",
          fill = "ONTOLOGY",               # change fill color by cyl
          color = "white",            # Set bar border colors to white
          palette = c('#5FB404','#12B5EC','#C238E5'),            # jco journal color palett. see ?ggpar
          sort.val = "asc",           # Sort the value in dscending order
          sort.by.groups = TRUE,      # Sort inside each group
          x.text.angle = 75, # Rotate vertically x axis texts
          ylab = '-log10(P-Value)',
          xlab = 'Pathway'
)
ggsave(filename ='GO triple 2.pdf',width = 7.23,height = 10)
ggsave(filename ='GO triple 2.TIFF',width = 7.23,height = 10)