# ============================================================
# BLCA Cuproptosis Prognostic Model — Fully Reproducible Code
# Author: Yuhan Zhang
# Affiliation: Guangling College, Yangzhou University
# ============================================================

# 0. Environment setup =======================================
setwd("C:/Users/云丘/Desktop")
library(data.table)
library(survival)
library(glmnet)

# 1. Load expression matrix ===================================
expr <- fread("TCGA.BLCA.sampleMap_HiSeqV2.gz", data.table = FALSE)
rownames(expr) <- expr[, 1]
expr <- expr[, -1]

# 2. Load clinical data =======================================
clinical <- fread("TCGA.BLCA.sampleMap_BLCA_clinicalMatrix", data.table = FALSE)

# 3. Define cuproptosis-related genes =========================
cuproptosis_genes <- c("FDX1", "LIAS", "LIPT1", "DLD", "DLAT",
                       "PDHA1", "PDHB", "MTF1", "GLS", "CDKN2A")

# 4. Truncate TCGA sample IDs to 12 characters ================
colnames(expr) <- substr(colnames(expr), 1, 12)

# 5. Build survival data frame ================================
clin <- data.frame(
  sampleID = clinical$bcr_patient_barcode,
  vital_status = clinical$vital_status,
  days_to_death = as.numeric(as.character(clinical$days_to_death)),
  days_to_last_followup = as.numeric(as.character(clinical$days_to_last_followup)),
  stringsAsFactors = FALSE
)
clin$OS.time <- ifelse(clin$vital_status == "DECEASED",
                       clin$days_to_death,
                       clin$days_to_last_followup)
clin$OS.event <- ifelse(clin$vital_status == "DECEASED", 1, 0)
clin <- clin[!is.na(clin$OS.time) & clin$OS.time > 0, ]

# 6. Match samples between expression and clinical data =========
common_samples <- intersect(colnames(expr), clin$sampleID)
clin <- clin[match(common_samples, clin$sampleID), ]
expr <- expr[, common_samples]

# 7. Extract cuproptosis gene expression =======================
available_genes <- intersect(cuproptosis_genes, rownames(expr))
cu_expr <- as.data.frame(t(expr[available_genes, ]))
cu_expr$sampleID <- rownames(cu_expr)
model_df <- merge(clin, cu_expr, by = "sampleID")
rownames(model_df) <- model_df$sampleID
cat("Training cohort size:", nrow(model_df), "\n")

# 8. Univariate Cox regression =================================
uni_cox <- data.frame()
for (gene in available_genes) {
  fit <- coxph(Surv(OS.time, OS.event) ~ model_df[[gene]], data = model_df)
  s <- summary(fit)
  uni_cox <- rbind(uni_cox, data.frame(
    Gene = gene,
    HR   = s$conf.int[1],
    LCI  = s$conf.int[3],
    UCI  = s$conf.int[4],
    P    = s$coefficients[5]
  ))
}
cat("\nUnivariate Cox regression (P < 0.1):\n")
print(uni_cox[uni_cox$P < 0.1, ])

# 9. LASSO-penalized Cox regression ============================
sig_genes <- uni_cox$Gene[uni_cox$P < 0.1]
if (length(sig_genes) < 2) {
  sig_genes <- uni_cox$Gene[order(uni_cox$P)[1:3]]
}
x <- as.matrix(model_df[, sig_genes])
y <- Surv(model_df$OS.time, model_df$OS.event)
set.seed(2025)
cv_lasso <- cv.glmnet(x, y, family = "cox", alpha = 1, nfolds = 10)
cat("Optimal lambda:", cv_lasso$lambda.min, "\n")
lasso_coef <- coef(cv_lasso, s = "lambda.min")
lasso_genes <- rownames(lasso_coef)[lasso_coef[, 1] != 0]
cat("LASSO-selected genes:", paste(lasso_genes, collapse = ", "), "\n")

# 10. Multivariate Cox model & risk score =======================
multi_cox <- coxph(as.formula(paste("Surv(OS.time, OS.event) ~",
                                     paste(lasso_genes, collapse = " + "))),
                   data = model_df)
cat("\nMultivariate Cox coefficients:\n")
print(summary(multi_cox)$coefficients)

model_df$riskScore <- predict(multi_cox, type = "risk")
model_df$riskGroup <- ifelse(model_df$riskScore > median(model_df$riskScore),
                             "High", "Low")

# 11. Figure 1 — KM curve in TCGA training cohort ===============
library(survminer)
fit_tcga <- survfit(Surv(OS.time, OS.event) ~ riskGroup, data = model_df)
pdf("Figure1_TCGA_KM.pdf", width = 8, height = 6)
ggsurvplot(fit_tcga, data = model_df, pval = TRUE, risk.table = TRUE,
           palette = c("#E64B35", "#4DBBD5"),
           title = "TCGA-BLCA Training Cohort")
dev.off()

# 12. Independent prognostic analysis (Table 2) =================
clinical_sub <- clinical[, c("bcr_patient_barcode",
                             "age_at_initial_pathologic_diagnosis")]
colnames(clinical_sub)[2] <- "Age"
temp_df <- merge(model_df, clinical_sub,
                 by.x = "sampleID", by.y = "bcr_patient_barcode")
temp_df$Age <- as.numeric(as.character(temp_df$Age))
fit_indep <- coxph(Surv(OS.time, OS.event) ~ riskScore + Age, data = temp_df)
cat("\nIndependent prognostic analysis (Table 2):\n")
print(summary(fit_indep)$coefficients)

# ==============================================================
# 13. External validation — GSE32894 ============================
# ==============================================================
library(GEOquery)
gse_val <- getGEO("GSE32894", GSEMatrix = TRUE, getGPL = FALSE)
expr_val <- exprs(gse_val[[1]])
clinical_val <- pData(gse_val[[1]])

# Extract survival information
os_time <- as.numeric(as.character(clinical_val$`time_to_dod_(months):ch1`))
os_event_raw <- as.character(clinical_val$`dod_event_(yes/no):ch1`)
os_event <- ifelse(grepl("yes", os_event_raw, ignore.case = TRUE), 1, 0)

# Probe IDs for LIPT1 and DLAT in GPL6947
probes_lipt1 <- c("ILMN_1717524", "ILMN_2343105")
probes_dlat  <- "ILMN_1706583"
expr_LIPT1 <- colMeans(expr_val[probes_lipt1, , drop = FALSE], na.rm = TRUE)
expr_DLAT  <- as.numeric(expr_val[probes_dlat, ])

# Build validation data frame
val_df <- data.frame(
  sampleID = colnames(expr_val),
  OS.time  = os_time,
  OS.event = os_event,
  LIPT1    = expr_LIPT1,
  DLAT     = expr_DLAT,
  stringsAsFactors = FALSE
)
val_df <- val_df[!is.na(val_df$OS.time) & val_df$OS.time > 0 &
                 !is.na(val_df$LIPT1) & !is.na(val_df$DLAT), ]

# Apply the same risk score formula from TCGA
val_df$riskScore <- (-0.607 * val_df$LIPT1) + (0.283 * val_df$DLAT)
val_df$riskGroup <- ifelse(val_df$riskScore > median(val_df$riskScore),
                           "High", "Low")

# Figure 2 — External validation KM curve
pdf("Figure2_GSE32894_KM.pdf", width = 8, height = 6)
fit_val <- survfit(Surv(OS.time, OS.event) ~ riskGroup, data = val_df)
ggsurvplot(fit_val, data = val_df, pval = TRUE, risk.table = TRUE,
           palette = c("#E64B35", "#4DBBD5"),
           title = "GSE32894 External Validation")
dev.off()

# ==============================================================
# 14. Enrichment analysis & ROC =================================
# ==============================================================
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(timeROC)

# Filter genes with zero variance
expr_var <- apply(expr, 1, var, na.rm = TRUE)
expr_filt <- expr[expr_var > 0, ]

# Identify genes co-expressed with LIPT1 or DLAT
coexp_genes <- c()
for (i in 1:nrow(expr_filt)) {
  ct1 <- cor.test(as.numeric(expr_filt[i, ]),
                  as.numeric(expr_filt["LIPT1", ]))
  ct2 <- cor.test(as.numeric(expr_filt[i, ]),
                  as.numeric(expr_filt["DLAT", ]))
  if (!is.na(ct1$estimate) & !is.na(ct2$estimate)) {
    if ((abs(ct1$estimate) > 0.4 & ct1$p.value < 0.05) |
        (abs(ct2$estimate) > 0.4 & ct2$p.value < 0.05)) {
      coexp_genes <- c(coexp_genes, rownames(expr_filt)[i])
    }
  }
}
cat("Number of co-expressed genes:", length(coexp_genes), "\n")

# GO enrichment
ego <- enrichGO(gene = coexp_genes, OrgDb = org.Hs.eg.db,
                keyType = "SYMBOL", ont = "BP",
                pAdjustMethod = "BH", qvalueCutoff = 0.05)
gene_entrez <- bitr(coexp_genes, fromType = "SYMBOL",
                    toType = "ENTREZID", OrgDb = org.Hs.eg.db)
ekegg <- enrichKEGG(gene = gene_entrez$ENTREZID, organism = "hsa",
                    pAdjustMethod = "BH", qvalueCutoff = 0.05)

# Figure 3A — GO enrichment dotplot
pdf("Figure3A_GO.pdf", width = 10, height = 8)
print(dotplot(ego, showCategory = 15,
              title = "GO Biological Process Enrichment"))
dev.off()

# Figure 3B — KEGG enrichment dotplot
pdf("Figure3B_KEGG.pdf", width = 10, height = 6)
print(dotplot(ekegg, showCategory = 15,
              title = "KEGG Pathway Enrichment"))
dev.off()

# Figure 4 — Time-dependent ROC curves
roc_out <- timeROC(T = model_df$OS.time, delta = model_df$OS.event,
                   marker = model_df$riskScore, cause = 1,
                   weighting = "marginal", times = c(365, 1095, 1825))
pdf("Figure4_ROC.pdf", width = 6, height = 6)
plot(roc_out, time = 365, col = "red", lwd = 2)
plot(roc_out, time = 1095, col = "blue", lwd = 2, add = TRUE)
plot(roc_out, time = 1825, col = "green", lwd = 2, add = TRUE)
legend("bottomright",
       paste(c("1-year", "3-year", "5-year"),
             "AUC =", round(roc_out$AUC, 3)),
       col = c("red", "blue", "green"), lwd = 2)
dev.off()

# ==============================================================
# 15. Save tables as CSV =======================================
# ==============================================================
write.csv(uni_cox, "Table1_Univariate_Cox.csv", row.names = FALSE)
write.csv(as.data.frame(summary(fit_indep)$coefficients),
          "Table2_Multivariate_Cox.csv")

cat("\n========== All analyses completed ==========\n")