rm(list = ls())
#Reference Package
library(reshape2)
library(ggpubr)
library(limma)
library(GSEABase)
library(GSVA)
expFile="merge.txt"           #Expression Input File
gmtFile="immune.gmt"          #Immunity dataset files
clusterFile="Cluster.txt"       #Typing input file

#Read expression input files and organize the input files
load("merge.RDATA")
rt=outTab
exp=outTab
data=exp

#Read gene set file
geneSets=getGmt(gmtFile, geneIdType=SymbolIdentifier())

#ssGSEA analyze
ssgseaScore=gsva(data, geneSets, method='ssgsea', kcdf='Gaussian', abs.ranking=TRUE)
#Correction of ssGSEA scores
normalize=function(x){
  return((x-min(x))/(max(x)-min(x)))}
ssgseaScore=normalize(ssgseaScore)
#Output ssGSEA scoring results
ssgseaOut=rbind(id=colnames(ssgseaScore), ssgseaScore)
write.table(ssgseaOut,file="ssGSEA.result.txt",sep="\t",quote=F,col.names=F)

#Read typing file
cluster=read.table(clusterFile, header=T, sep="\t", check.names=F, row.names=1)

#Data Merge
ssgseaScore=t(ssgseaScore)
sameSample=intersect(row.names(ssgseaScore), row.names(cluster))
ssgseaScore=ssgseaScore[sameSample,,drop=F]
cluster=cluster[sameSample,,drop=F]
scoreCluster=cbind(ssgseaScore, cluster)

#Convert data into ggplot2 input file
data=melt(scoreCluster, id.vars=c("cluster"))
colnames(data)=c("cluster", "Immune", "Fraction")

#Draw a box plot
bioCol=c("#0066FF","#FF9900","#FF0000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length(levels(factor(data[,"cluster"])))]
library(tidyverse)
data$Immune=str_replace_all(data$Immune,"na","")
p=ggboxplot(data, x="Immune", y="Fraction", color="cluster", 
            ylab="Immune infiltration",
            xlab="",
            legend.title="cluster",
            palette=bioCol)
p=p+rotate_x_text(50)
pdf(file="boxplot.pdf", width=10, height=6.5)                          #Output image file
p+stat_compare_means(aes(group=cluster),symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),label = "p.signif")
dev.off()