#Reference Package
rm(list = ls())
library(limma)
library(estimate)
library(tidyverse)
inputFile="merge.txt"                                          #Enter the file name

#Read the file and organize the input file
load("merge.RDATA")
rt=outTab


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

#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 the score for each sample
scores=read.table("estimateScore.gct",skip = 2,header = T)
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=F,col.names=F)

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

#Data Merge
sameSample=intersect(row.names(scores), row.names(cluster))
scores=scores[sameSample,,drop=F]
cluster=cluster[sameSample,,drop=F]
scoreCluster=cbind(scores, cluster)
library(reshape2)
#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)
library(ggpubr)
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=4, height=4.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()