rm(list = ls())

#Reference Package
library(limma)
library(GSEABase)
library(GSVA)
library(pheatmap)
expFile="merge.txt"               #Expression Input File
clusterFile="Cluster.txt"      #Typing input file
gmtFile="h.all.v2023.1.Hs.symbols.gmt"                    #Gene set files


#Read expression input files and organize the input files
load("merge.RDATA")
rt=outTab
exp=outTab
dimnames=list(rownames(exp), colnames(exp))
data=matrix(as.numeric(as.matrix(exp)), nrow=nrow(exp), dimnames=dimnames)
data=avereps(data)

#GSVA analyze
geneSets=getGmt(gmtFile, geneIdType=SymbolIdentifier())
gsvaResult=gsva(data, 
                geneSets, 
                min.sz=10, 
                max.sz=500, 
                verbose=TRUE,
                parallel.sz=1)
gsvaOut=rbind(id=colnames(gsvaResult), gsvaResult)
write.table(gsvaOut, file="gsvaOut.txt", sep="\t", quote=F, col.names=F)

#Reading cluster files
cluster=read.table(clusterFile, header=T, sep="\t", check.names=F, row.names=1)

#Data Merge
gsvaResult=t(gsvaResult)
sameSample=intersect(row.names(gsvaResult), row.names(cluster))
gsvaResult=gsvaResult[sameSample,,drop=F]
cluster=cluster[sameSample,,drop=F]
gsvaCluster=cbind(gsvaResult, cluster)
Project=gsub("(.*?)\\_.*", "\\1", rownames(gsvaCluster))
gsvaCluster=cbind(gsvaCluster, Project)
colnames(gsvaCluster)[which(colnames(gsvaCluster)=="cluster")]="cluster"
#Variance Analysis
adj.P.Val.Filter=0.05
allType=as.vector(gsvaCluster$cluster)
comp=combn(levels(factor(allType)), 2)
for(i in 1:ncol(comp)){
  #Sample grouping
  treat=gsvaCluster[gsvaCluster$cluster==comp[2,i],]
  con=gsvaCluster[gsvaCluster$cluster==comp[1,i],]
  data=rbind(con, treat)
  #Variance Analysis
  Type=as.vector(data$cluster)
  ann=data[,c(ncol(data), (ncol(data)-1))]
  data=t(data[,-c((ncol(data)-1), ncol(data))])
  design=model.matrix(~0+factor(Type))
  colnames(design)=levels(factor(Type))
  fit=lmFit(data, design)
  contrast=paste0(comp[2,i], "-", comp[1,i])
  cont.matrix=makeContrasts(contrast, levels=design)
  fit2=contrasts.fit(fit, cont.matrix)
  fit2=eBayes(fit2)
  
  #Output the differences of all paths
  allDiff=topTable(fit2,adjust='fdr',number=200000)
  allDiffOut=rbind(id=colnames(allDiff),allDiff)
  write.table(allDiffOut, file=paste0(contrast, ".all.txt"), sep="\t", quote=F, col.names=F)
  
  #Significant difference in output
  diffSig=allDiff[with(allDiff, (abs(logFC)>0.1 & adj.P.Val < adj.P.Val.Filter )), ]
  diffSigOut=rbind(id=colnames(diffSig),diffSig)
  write.table(diffSigOut, file=paste0(contrast, ".diff.txt"), sep="\t", quote=F, col.names=F)
  
  #Cluster Color
  bioCol=c("#0066FF","#FF9900","#FF0000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
  ann_colors=list()
  CluCol=bioCol[1:length(levels(factor(allType)))]
  names(CluCol)=levels(factor(allType))
  ann_colors[["cluster"]]=CluCol[c(comp[1,i], comp[2,i])]
  
  #Draw a heat map of differential pathways
  termNum=20
  diffTermName=as.vector(rownames(diffSig))
  diffLength=length(diffTermName)
  if(diffLength<termNum){termNum=diffLength}
  hmGene=diffTermName[1:termNum]
  hmExp=data[hmGene,]
  pdf(file=paste0(contrast,".heatmap.pdf"),height=6,width=10)
  pheatmap(hmExp, 
           annotation=ann,
           annotation_colors = ann_colors,
           color = colorRampPalette(c(rep("blue",2), "white", rep("red",2)))(50),
           cluster_cols =F,
           show_colnames = F,
           gaps_col=as.vector(cumsum(table(Type))),
           scale="row",
           fontsize = 10,
           fontsize_row=7,
           fontsize_col=10)
  dev.off()
}