# Clear the current workspace
rm(list = ls())

# Load required packages
library(limma)
library(GSEABase)
library(GSVA)
library(pheatmap)

# Define input files
expFile <- "merge.txt"               # Expression input file
clusterFile <- "Cluster.txt"         # Clustering input file
gmtFile <- "c2.cp.reactome.v7.5.1.symbols.gmt"   # Gene set file

# Read the expression input file and preprocess
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)

# Perform GSVA analysis
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 = FALSE, col.names = FALSE)

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

# Merge data
gsvaResult <- t(gsvaResult)
sameSample <- intersect(row.names(gsvaResult), row.names(cluster))
gsvaResult <- gsvaResult[sameSample,,drop = FALSE]
cluster <- cluster[sameSample,,drop = FALSE]
gsvaCluster <- cbind(gsvaResult, cluster)
Project <- gsub("(.*?)\\_.*", "\\1", rownames(gsvaCluster))
gsvaCluster <- cbind(gsvaCluster, Project)
colnames(gsvaCluster)[which(colnames(gsvaCluster) == "cluster")] <- "cluster"

# Differential 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)
  
  # Differential 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 differential results for all pathways
  allDiff <- topTable(fit2, adjust = 'fdr', number = 200000)
  allDiffOut <- rbind(id = colnames(allDiff), allDiff)
  write.table(allDiffOut, file = paste0(contrast, ".all.txt"), sep = "\t", quote = FALSE, col.names = FALSE)
  
  # Output significant differential results
  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 = FALSE, col.names = FALSE)
  
  # Cluster colors
  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])]
  
  # Generate heatmap for 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("green", 2), "white", rep("red", 2)))(50),
           cluster_cols = FALSE,
           show_colnames = FALSE,
           gaps_col = as.vector(cumsum(table(Type))),
           scale = "row",
           fontsize = 10,
           fontsize_row = 7,
           fontsize_col = 10)
  dev.off()
}