rm(list = ls())
library(Seurat)
library(ggplot2)
library(cowplot)
library(Matrix)
library(dplyr)
library(ggsci)
library(scater)
library(SingleR)
library(hdf5r)
library(tidyverse)
library(RColorBrewer)
library(dplyr)
library(Seurat)
library(ComplexHeatmap)
library(GSVA)
library(GSEABase)
library(limma)
library(ggplot2)
load("merge.RDATA")

score=read.table("score.txt",header = T,sep = "\t")


genesets <- getGmt("h.all.v2023.1.Hs.symbols.gmt")
str(genesets)

gsvascore <- gsva(outTab, 
                  genesets, 
                  min.sz=10, 
                  max.sz=500, 
                  verbose=TRUE,
                  parallel.sz=1)


aaa=as.data.frame(t(gsvascore))
aaa$id=rownames(aaa)
aaa=inner_join(score,aaa,by="id")

gene="score"
y <- as.numeric(aaa[,gene])
bbb=aaa[,3:52]
cor_data_df <- data.frame(colnames(bbb))
for (a in 1:length(colnames(bbb))){
  test <- cor.test(as.numeric(bbb[,a]),y,method = "pearson")
  cor_data_df[a,2] <- test$estimate
  cor_data_df[a,3] <- test$p.value
}
names(cor_data_df) <- c("symbol","correlation","pvalue")
write.csv (cor_data_df, file =paste(gene,'GSVA correlation analysis.csv',sep=' '), row.names =FALSE) 
cor_data_df<-na.omit(cor_data_df)
posneg<-cor_data_df %>%
  filter(symbol != gene)
posneg$group<-NA
posneg$group[ posneg$pvalue<0.05& posneg$correlation<0]<-"negative"
posneg$group[ posneg$pvalue>=0.05]<-"pvalue>0.05"
posneg$group[ posneg$pvalue<0.05& posneg$correlation>0]<-"positive"
posneg$symbol<-str_replace_all(posneg$symbol,'_',' ')
posneg$symbol<-str_replace_all(posneg$symbol,'HALLMARK ','')
library(ggpubr)
p<-ggbarplot(posneg, x = "symbol", y = "correlation",
             fill = "group",           # change fill color by mpg_level
             color = "white",            # Set bar border colors to white
             palette = "jco",            # jco journal color palett. see ?ggpar
             sort.val = "asc",           # Sort the value in ascending order
             sort.by.groups = FALSE,     # Don't sort inside each group
             #x.text.angle = 90,          # Rotate vertically x axis texts
             ylab = "GSVA correlation",
             xlab = "Term",
             legend.title = "group",
             lab.col = "black",
             rotate = TRUE,
             ggtheme = theme_minimal()
)
ggsave(p,filename = paste(gene,'GSVA Correlation Plot.pdf',sep=' '),width = 10,height = 10)