# Clear the workspace
rm(list = ls())

# Load required packages
library(reshape2)
library(ggpubr)
library(limma)
library(GSEABase)
library(GSVA)

# Define input file names
expFile <- "merge.txt"
gmtFile <- "immune.gmt"
clusterFile <- "Cluster.txt"

# Read the expression input file and preprocess
load("merge.RDATA")
rt <- outTab
exp <- outTab
data <- exp

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

# Perform ssGSEA analysis
ssgseaScore <- gsva(data, geneSets, method = 'ssgsea', kcdf = 'Gaussian', abs.ranking = TRUE)

# Normalize ssGSEA scores
normalize <- function(x) {
  return((x - min(x)) / (max(x) - min(x)))
}
ssgseaScore <- normalize(ssgseaScore)

# Output ssGSEA score results
ssgseaOut <- rbind(id = colnames(ssgseaScore), ssgseaScore)
write.table(ssgseaOut, file = "ssGSEA.result.txt", sep = "\t", quote = FALSE, col.names = FALSE)

# Read the clustering file
cluster <- read.table(clusterFile, header = TRUE, sep = "\t", check.names = FALSE, row.names = 1)

# Merge data
ssgseaScore <- t(ssgseaScore)
sameSample <- intersect(row.names(ssgseaScore), row.names(cluster))
ssgseaScore <- ssgseaScore[sameSample,,drop = FALSE]
cluster <- cluster[sameSample,,drop = FALSE]
scoreCluster <- cbind(ssgseaScore, cluster)

# Convert data to ggplot2 input format
data <- melt(scoreCluster, id.vars = c("cluster"))
colnames(data) <- c("cluster", "Immune", "Fraction")

# Create a boxplot
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 <- rotate_x_text(p, angle = 50)

# Output the boxplot as a PDF file
pdf(file = "boxplot.pdf", width = 10, height = 6.5)
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()