# Load required packages
rm(list = ls())
library(limma)
library(estimate)
library(tidyverse)

# Define input file name
inputFile <- "merge.txt"

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

# Output the processed matrix file
rt <- rbind(ID = colnames(rt), rt)
write.table(rt, file = "uniq.symbol.txt", sep = "\t", quote = FALSE, col.names = FALSE)

# Run the estimate package
filterCommonGenes(input.f = "uniq.symbol.txt",
                  output.f = "commonGenes.gct",
                  id = "GeneSymbol")

estimateScore(input.ds = "commonGenes.gct",
              output.ds = "estimateScore.gct",
              platform = "illumina")

# Output scores for each sample
scores <- read.table("estimateScore.gct", skip = 2, header = TRUE)
rownames(scores) <- scores[, 1]
scores <- t(scores[, 3:ncol(scores)])
rownames(scores) <- gsub("\\.", "\\-", rownames(scores))
out <- rbind(ID = colnames(scores), scores)
write.table(out, file = "scores.txt", sep = "\t", quote = FALSE, col.names = FALSE)

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

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

# Convert data to ggplot2 input format
library(reshape2)
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(ggpubr)
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 = 4, height = 4.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()