rm(list = ls())
#Reference Package
library(limma)
library(reshape2)
library(ggpubr)
expFile="uniSigGeneExp.txt"          #Expression Input File
geneCluFile="geneCluster.txt"     #Genotyping files

#Reading expression input files
rt=read.table(expFile, header=T, sep="\t", check.names=F)
exp=rt
rownames(exp)=exp$id
exp=exp[,-1]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=t(data)

#Read genotyping files
geneClu=read.table(geneCluFile, header=T, sep="\t", check.names=F, row.names=1)

#Merge data
sameSample=intersect(row.names(data), row.names(geneClu))
expClu=cbind(data[sameSample,,drop=F], geneClu[sameSample,,drop=F])

#Convert data into ggplot2 input file
data=melt(expClu, id.vars=c("geneCluster"))
colnames(data)=c("geneCluster", "Gene", "Expression")

#Setting Color
bioCol=c("#0066FF","#FF9900","#FF0000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length(levels(factor(data[,"geneCluster"])))]

#Draw a box plot
p=ggboxplot(data, x="Gene", y="Expression", color = "geneCluster", 
            ylab="Gene expression",
            xlab="",
            legend.title="geneCluster",
            palette = bioCol,
            width=1)
p=p+rotate_x_text(60)
p1=p+stat_compare_means(aes(group=geneCluster),
                        symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),
                        label = "p.signif")

#Output box plot
pdf(file="boxplot.pdf", width=10, height=5)
print(p1)
dev.off()