# Clear the R environment
rm(list = ls())

# Load required packages
library(limma)
library(ggpubr)

# Expression data file
expFile = "merge.txt"
# Score grouping file
scoreFile = "score.group.txt"

# Read the expression data file
load("merge.RDATA")
rt = outTab
data = outTab
colnames(data) = gsub("(.*?)\\_(.*?)", "\\2", colnames(data))
data2 = data

# Target gene's standard name
gene = "CTLA4"

# Extract expression levels for the target gene
data = rbind(data2, gene = data2[gene,])
exp = t(data[c("gene", gene),])
exp = avereps(exp)

# Read the score grouping file
score = read.table(scoreFile, header = TRUE, sep = "\t", check.names = FALSE, row.names = 1)

# Merge data
sameSample = intersect(row.names(exp), row.names(score))
exp = exp[sameSample,]
score = score[sameSample,]
data = cbind(as.data.frame(exp), as.data.frame(score))

# Set up comparison groups
data$group = factor(data$group, levels = c("Low", "High"))
group = levels(factor(data$group))
comp = combn(group, 2)
my_comparisons = list()
for (i in 1:ncol(comp)) {
  my_comparisons[[i]] <- comp[, i]
}

# Create a boxplot
boxplot = ggboxplot(data, x = "group", y = "gene", fill = "group",
                    xlab = "score",
                    ylab = paste(gene, "expression"),
                    legend.title = "score",
                    palette = c("#00AFBB", "#E7B800")) + 
  stat_compare_means(comparisons = my_comparisons)

# Output the plot as a PDF file
pdf(file = paste0(gene, ".pdf"), width = 5, height = 4.5)
print(boxplot)
dev.off()