rm(list = ls())

library(org.Bt.eg.db)
library(stringr)
library(BiocGenerics)
library(clusterProfiler)
library(ggplot2)
require(DOSE)
library(Hmisc)
library(future)
library(future.apply)
library(clusterProfiler)
library(limma)
library(enrichplot)
library(psych)
library(RColorBrewer)
library(GSVA)


matrix_uniq.filt.tumor <- read.table("GSEA-GSE137943.txt", header = T, row.names = 1)
matrix_uniq.filt.tumor[1:5,1:5]
###
data01 <- c("RYR1")
data01

tar.exp <- matrix_uniq.filt.tumor[data01,]
head(tar.exp)

class(tar.exp)
# "data.frame"
y <- as.numeric(tar.exp)
class(y)
#"numeric"

data1 <- data.frame()
for (i in rownames(matrix_uniq.filt.tumor)){
  dd <- corr.test(as.numeric(matrix_uniq.filt.tumor[i,]), y, method = "pearson", adjust = "fdr")
  data1 <- rbind(data1, data.frame(gene =i, cor = dd$r, p.value=dd$p))  
}
data1 <- data1[order(data1$cor,decreasing = T),]
gene1 <- data1$cor
names(gene1) <- mapIds(org.Bt.eg.db,keys = data1$gene,column = 'ENTREZID',
                       keytype = 'SYMBOL',multiVals='filter')
gene1 <- na.omit(gene1)
save(gene1,file = paste0('RYR1.cor.Rdata'))
#分析
options(digits = 4)
#install.packages("R.utils")

library(R.utils)

R.utils::setOption("clusterProfiler.download.method","auto")

kegggsea <-  gseKEGG(
  gene1,
  organism = "bta",
  keyType = "kegg",
  pvalueCutoff = 0.05,
  pAdjustMethod = "BH",
  verbose = TRUE,
  seed = TRUE)
write.table(kegggsea,'RYR1.KEGGgsea.txt',sep = '\t',quote = F,row.names = F)

pdf('RYR1.GSEA.KEGG.pdf',width = 15,height = 10)
gseaplot2(kegggsea, geneSetID = c("bta04910","bta04020","bta04010","bta04152","bta00010"),pvalue_table = FALSE,
          base_size = 13,  
          title = "RYR1_GSEA",
          rel_heights = c(1.5, 0.5, 0.5),color= brewer.pal(10,'Paired'))

dev.off()