# ============================================================
# 教程配套脚本 09：独立验证（GSE14520 验证 TCGA 发现）
# 核心思想：训练集（TCGA）构建的结论，在独立队列（GEO）中重复检验
# ============================================================

library(survival)
library(survminer)
library(tidyverse)

# ---- 1. 读取 GSE14520 数据（若尚未处理，先跑脚本 02/03）----
expr_geo <- read.csv("data/processed/GSE14520_expr_symbol.csv",
                     row.names = 1, check.names = FALSE)
pheno_geo <- read.csv("data/raw/GSE14520_pheno_raw.csv",
                      row.names = 1, check.names = FALSE)

# GSE14520 的生存信息在表型列中（列名以实际为准，常见为
# "survival months" 与 "survival" 或需要从 supplementary 获取）
# 示例（占位，正式分析按实际列名调整）：
# OS.time <- as.numeric(pheno_geo$`survival months:ch1`)
# OS <- as.numeric(pheno_geo$`survival:ch1`)

# ---- 2. 验证 hub 基因表达（示例：GSDMD 在肿瘤 vs 癌旁）----
# group 来自脚本 03 保存的分组
group_geo <- readRDS("data/processed/GSE14520_group.rds")
keep <- group_geo %in% c("Tumor", "Normal")

if ("GSDMD" %in% rownames(expr_geo)) {
  df_box <- data.frame(
    expression = as.numeric(expr_geo["GSDMD", keep]),
    group = factor(group_geo[keep], levels = c("Normal", "Tumor"))
  )
  p_box <- ggplot(df_box, aes(x = group, y = expression, fill = group)) +
    geom_boxplot(alpha = 0.8) +
    geom_jitter(width = 0.15, size = 0.6, alpha = 0.4) +
    stat_compare_means(method = "wilcox.test", label = "p.signif") +
    scale_fill_manual(values = c("#4C72B0", "#C44E52")) +
    labs(title = "GSDMD expression: GSE14520 validation",
         x = NULL, y = "Expression (log2)") +
    theme_minimal() +
    theme(legend.position = "none")
  ggsave("figures/R_validation_boxplot.png", p_box, width = 5.5, height = 5, dpi = 300)
}

# ---- 3. 验证风险模型（若有生存数据）----
# 思路：把 TCGA 训练好的 LASSO 系数应用到 GEO 队列表达，
#       计算 GEO 患者风险评分 → 分组 → KM 曲线 + AUC。
# 示例骨架（正式分析按实际数据结构补全）：
# X_geo <- t(expr_geo[selected_genes, ])           # 样本 × 基因
# risk_geo <- as.numeric(X_geo %*% lasso_coef)     # 应用训练系数
# clin_geo <- data.frame(OS.time = ..., OS = ..., risk = risk_geo)
# fit_km_geo <- survfit(Surv(OS.time, OS) ~ risk_high, data = clin_geo)
# ggsurvplot(fit_km_geo, pval = TRUE)

cat("独立验证脚本就绪：请按实际数据列名补全生存字段后运行。\n")
