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

# Load libraries
library(limma)
library(survival)
library(tidyverse)

# Define input files
expFile <- "interGeneExp.txt"   # Differential gene expression file
cliFile <- "clinicaldata.txt"    # Survival data file

# Read and preprocess the expression file
rt <- read.table(expFile, header = TRUE, sep = "\t", check.names = FALSE)
rt <- as.matrix(rt)
rownames(rt) <- rt[, 1]
exp <- rt[, 2:ncol(rt)]
dimnames <- list(rownames(exp), colnames(exp))
data <- matrix(as.numeric(as.matrix(exp)), nrow = nrow(exp), dimnames = dimnames)
data <- avereps(data)
data1 <- t(data)
rownames(data1) <- str_replace_all(rownames(data1), "TCGA_", "")
rownames(data1) <- str_replace_all(rownames(data1), "GSE17538_", "")
rownames(data1) <- str_replace_all(rownames(data1), "GSE39582_", "")

# Read survival data
cli <- read.table(cliFile, header = TRUE, sep = "\t", check.names = FALSE, row.names = 1)
cli <- cli[, 1:2]
cli$futime <- cli$futime / 365

# Merge data
sameSample <- intersect(row.names(data1), row.names(cli))
data1 <- data1[sameSample,, drop = FALSE]
cli <- cli[sameSample,, drop = FALSE]
rt <- cbind(cli, data1)

# Perform Cox analysis for prognosis-related genes
outTab <- data.frame()
sigGenes <- c()
for (i in colnames(rt[, 3:ncol(rt)])) {
  # Cox analysis
  cox <- coxph(Surv(futime, fustat) ~ rt[, i], data = rt)
  coxSummary <- summary(cox)
  coxP <- coxSummary$coefficients[,"Pr(>|z|)"]
  if (coxP < 100) {
    outTab <- rbind(outTab,
                    cbind(id = i,
                          HR = coxSummary$conf.int[,"exp(coef)"],
                          HR.95L = coxSummary$conf.int[,"lower .95"],
                          HR.95H = coxSummary$conf.int[,"upper .95"],
                          pvalue = coxSummary$coefficients[,"Pr(>|z|)"])
    )
  }
}

# Output results for single-factor analysis
write.table(outTab, file = "uniCox.txt", sep = "\t", row.names = FALSE, quote = FALSE)
outTab$pvalue <- as.numeric(outTab$pvalue)
aaa <- dplyr::filter(outTab, pvalue < 0.00001)
sigGenes <- aaa$id

# Save expression data for significant genes
sigGeneExp <- data[sigGenes,]
sigGeneExp <- rbind(id = colnames(sigGeneExp), sigGeneExp)
write.table(sigGeneExp, file = "uniSigGeneExp.txt", sep = "\t", quote = FALSE, col.names = FALSE)

# Load additional libraries
library(survminer)
library(survival)
library(survivalROC)
library(forestplot)

rt <- aaa
row.names(rt) <- rt$id
rt <- rt[, -1]
rt <- na.omit(rt)
rt <- arrange(rt, pvalue)

options(forestplot_new_page = FALSE)
clrs <- fpColors(box = "red", line = "darkblue", summary = "royalblue")

data <- as.matrix(rt)
HR <- as.data.frame(data[, 1:3])
HR[, 1] <- as.numeric(HR[, 1])
HR[, 2] <- as.numeric(HR[, 2])
HR[, 3] <- as.numeric(HR[, 3])
hr <- sprintf("%.3f", HR[,"HR"])
hrLow <- sprintf("%.3f", HR[,"HR.95L"])
hrHigh <- sprintf("%.3f", HR[,"HR.95H"])
pVal <- data[,"pvalue"]
pVal <- as.numeric(pVal)
pVal <- ifelse(pVal < 0.00001, "<0.00001", sprintf("%.3f", pVal))
tabletext <- 
  list(c(NA, rownames(HR)),
       append("pvalue", pVal),
       append("Hazard ratio", paste0(hr, "(", hrLow, "-", hrHigh, ")"))
  )

pdf(file = "OS_forest.pdf",
    width = 6,
    height = 8
)

forestplot(tabletext, 
           zero = 1,
           lwd.zero = 2,
           rbind(rep(NA, 3), HR),
           col = clrs,
           graphwidth = unit(50, "mm"),
           xlog = TRUE,
           lwd.ci = 2,
           boxsize = 0.2,
           xlab = "Hazard ratio"
)

dev.off()