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

# Load required packages
library(limma)
library(VennDiagram)

# Define input file names
expFile <- "merge.txt"
cluFile <- "Cluster.txt"

# Read the 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)

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

# Extract intersecting samples
sameSample <- intersect(colnames(data), row.names(cluster))
data <- data[, sameSample]
cluster <- cluster[sameSample,]

# Differential analysis
logFCfilter <- 1
geneList <- list()
Type <- as.vector(cluster)
design <- model.matrix(~0 + factor(Type))
colnames(design) <- levels(factor(Type))
comp <- combn(levels(factor(Type)), 2)
allDiffGenes <- c()

for (i in 1:ncol(comp)) {
  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 all genes' differential information
  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 differential results
  diffSig <- allDiff[with(allDiff, (abs(logFC) > logFCfilter & P.Value < 0.05 )), ]
  diffSigOut <- rbind(id = colnames(diffSig), diffSig)
  write.table(diffSigOut, file = paste0(contrast, ".diff.txt"), sep = "\t", quote = FALSE, col.names = FALSE)
  geneList[[contrast]] <- row.names(diffSig)
}

# Save union genes
B.A <- geneList$`B-A`
interGenes <- unique(c(B.A))
write.table(file = "interGene.txt", interGenes, sep = "\t", quote = FALSE, col.names = FALSE, row.names = FALSE)

# Save expression levels of intersecting genes
interGeneExp <- data[interGenes,]
interGeneExp <- rbind(id = colnames(interGeneExp), interGeneExp)
write.table(interGeneExp, file = "interGeneExp.txt", sep = "\t", quote = FALSE, col.names = FALSE)

# Create volcano plot
volcano <- read.table("B-A.all.txt", header = TRUE, sep = "\t", row.names = 1)
p.cut <- 0.05
logFC.cut <- 1
volcano$type[(volcano$P.Value > p.cut | volcano$P.Value == "NA") | (volcano$logFC < logFC.cut & volcano$logFC > -logFC.cut)] <- "none significant"
volcano$type[volcano$P.Value <= p.cut & volcano$logFC >= logFC.cut] <- "up-regulated"
volcano$type[volcano$P.Value <= p.cut & volcano$logFC <= -logFC.cut] <- "down-regulated"

# Create a volcano plot
library(tidyverse)
p <- ggplot(volcano, aes(logFC, -1 * log10(P.Value), color = type))
p <- p + geom_point()

x_lim <- max(volcano$logFC, -allDiff$logFC) 
gg <- p + geom_point(aes(size = abs(logFC)), alpha = 0.4) + xlim(-x_lim, x_lim) + labs(x = "log2(FC)", y = "-log10(P.Value)") +
  scale_color_manual(values = c("green", "grey", "red")) +
  geom_hline(aes(yintercept = -1 * log10(p.cut)), colour = "black", linetype = "dashed") +
  geom_vline(xintercept = c(-logFC.cut, logFC.cut), colour = "black", linetype = "dashed")

# Print the volcano plot
pdf(file = "Volcano_Plot.pdf", width = 6, height = 6)
gg
dev.off()