# ============================================================
# External public-database analyses for the gastric-cancer LVI study
# From GEO/GDC download to statistical analysis and plotting
# ============================================================
# Evidence roles from the supplied manuscript:
#   Direct LVI: GSE62254 (RNA), GSE261635 (EPIC methylation)
#   Surrogate/sensitivity: TCGA-STAD, GSE246963, GSE237876
#   Single-cell sensitivity: GSE163558, GSE239676
#   In-silico prioritization: DepMap and LINCS (optional input files)
#
# IMPORTANT: the manuscript does not contain the original source code or the
# complete six gene lists. Replace the six PROGRAMS below with the exact
# prespecified lists used in the study before interpreting numerical results.
# GEO clinical fields also differ between releases; inspect pData() and set
# the mapping variables in section 4 if automatic matching is not correct.
# ============================================================

options(stringsAsFactors = FALSE)
set.seed(20260919)

# ---- 1. Packages -----------------------------------------------------------
cran <- c("tidyverse", "data.table", "pheatmap", "survival", "survminer",
          "broom", "sandwich", "lmtest", "effectsize", "ggrepel", "patchwork")
bioc <- c("GEOquery", "TCGAbiolinks", "SummarizedExperiment", "limma",
          "Biobase", "AnnotationDbi", "org.Hs.eg.db", "GSVA", "GSEABase")
install_if_missing <- function(pkgs, bioc = FALSE) {
  miss <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
  if (!length(miss)) return(invisible(NULL))
  if (bioc) {
    if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
    BiocManager::install(miss, ask = FALSE, update = FALSE)
  } else install.packages(miss, repos = "https://cloud.r-project.org")
}
install_if_missing(cran)
install_if_missing(bioc, bioc = TRUE)
invisible(lapply(c(cran, bioc), library, character.only = TRUE))

# ---- 2. Folders ------------------------------------------------------------
ROOT <- "public_external_analysis"
dir.create(ROOT, showWarnings = FALSE)
dir.create(file.path(ROOT, "raw"), showWarnings = FALSE)
dir.create(file.path(ROOT, "rds"), showWarnings = FALSE)
dir.create(file.path(ROOT, "tables"), showWarnings = FALSE)
dir.create(file.path(ROOT, "figures"), showWarnings = FALSE)

# ---- 3. Prespecified programs ---------------------------------------------
# Keep gene symbols uppercase. Replace the illustrative vectors with the exact
# lists used in the study; KEGG lists are downloaded reproducibly below.
PROGRAMS <- list(
  ketone_body = c("ACAT1","BDH1","BDH2","HMGCS2","HMGCL","OXCT1","SLC16A1","SLC16A7","ACAT2","HADHA","HADHB","PPARA"),
  butanoate = character(),          # filled from KEGG hsa00650
  fatty_acid_degradation = character(), # filled from KEGG hsa00071
  EMT_ECM = c("VIM","CDH2","SNAI1","SNAI2","TWIST1","TWIST2","ZEB1","ZEB2","MMP2","MMP9","ITGA5","ITGB1","FN1","COL1A1","COL3A1"),
  endothelial_angiogenic = c("KDR","FLT1","VEGFA","VEGFC","ANGPT1","ANGPT2","TEK","ESM1","PECAM1","VWF","ENG","EMCN","RAMP2","PLVAP","KLF2","EFNB2","DLL4","ROBO4"),
  epithelial_integrity = c("EPCAM","KRT8","KRT18","KRT19","CDH1","MUC1","CLDN3","CLDN4")
)
# KEGG REST is public and does not require an API key.
kegg_genes <- function(pathway) {
  u <- paste0("https://rest.kegg.jp/link/hsa/", pathway)
  z <- tryCatch(readLines(u, warn = FALSE), error = function(e) character())
  if (!length(z)) return(character())
  unique(sub("hsa:", "", sub(".*\\t", "", z)))
}
if (!length(PROGRAMS$butanoate)) PROGRAMS$butanoate <- kegg_genes("hsa00650")
# typo-safe assignment for the name used above
if (!length(PROGRAMS$fatty_acid_degradation)) PROGRAMS$fatty_acid_degradation <- kegg_genes("hsa00071")

# Convert KEGG Entrez IDs to symbols when needed.
id_to_symbol <- function(x) {
  x <- unique(as.character(x)); x <- x[nzchar(x)]
  out <- AnnotationDbi::select(org.Hs.eg.db, keys = x, keytype = "ENTREZID", columns = "SYMBOL")
  unique(na.omit(out$SYMBOL))
}
PROGRAMS$butanoate <- id_to_symbol(PROGRAMS$butanoate)
PROGRAMS$fatty_acid_degradation <- id_to_symbol(PROGRAMS$fatty_acid_degradation)

# ---- 4. GEO download and expression helpers -------------------------------
get_geo <- function(accession, destdir = file.path(ROOT, "raw")) {
  message("Downloading ", accession)
  g <- GEOquery::getGEO(accession, GSEMatrix = TRUE, getGPL = TRUE, destdir = destdir)
  if (is.list(g)) g <- g[[1L]]
  saveRDS(g, file.path(ROOT, "rds", paste0(accession, ".rds")))
  g
}
geo_feature_symbols <- function(eset) {
  fd <- as.data.frame(Biobase::fData(eset))
  nms <- tolower(colnames(fd))
  hit <- which(nms %in% c("gene symbol","gene_symbol","genesymbol","symbol","gene_assignment","gene"))
  if (!length(hit)) hit <- grep("symbol|gene_assignment", nms)
  if (!length(hit)) stop("No gene-symbol column found in featureData; inspect fData(eset).")
  s <- as.character(fd[[hit[1L]]])
  s <- sub(" ///.*$", "", s); s <- sub(" //.*$", "", s)
  s <- trimws(s); s[s %in% c("", "NA", "---")] <- NA
  s
}
eset_to_gene_matrix <- function(eset, log2_transform = TRUE) {
  x <- Biobase::exprs(eset)
  if (log2_transform && max(x, na.rm = TRUE) > 50) x <- log2(x + 1)
  sym <- geo_feature_symbols(eset)
  keep <- !is.na(sym)
  x <- x[keep, , drop = FALSE]; sym <- sym[keep]
  # collapse duplicate probes by median
  z <- rowsum(x, group = sym, reorder = FALSE)
  n <- table(sym)[rownames(z)]
  z / as.numeric(n)
}

# Public accessions in the supplied manuscript.
GEO_IDS <- c("GSE62254","GSE246963","GSE237876","GSE163558","GSE239676","GSE261635")
geo <- setNames(lapply(GEO_IDS, get_geo), GEO_IDS)

# ---- 5. Program scoring ----------------------------------------------------
score_programs <- function(expr, programs = PROGRAMS) {
  expr <- as.matrix(expr)
  expr <- expr[!duplicated(rownames(expr)), , drop = FALSE]
  out <- lapply(programs, function(gs) {
    gs <- intersect(toupper(gs), toupper(rownames(expr)))
    if (!length(gs)) return(rep(NA_real_, ncol(expr)))
    a <- t(scale(t(expr[gs, , drop = FALSE])))
    colMeans(a, na.rm = TRUE)
  })
  out <- as.data.frame(out, check.names = FALSE)
  rownames(out) <- colnames(expr)
  attr(out, "coverage") <- vapply(programs, function(gs) sum(toupper(gs) %in% toupper(rownames(expr))), integer(1))
  out
}

# robust extraction of a clinical variable from GEO pData
auto_field <- function(pd, patterns) {
  hit <- grep(paste(patterns, collapse = "|"), tolower(colnames(pd)))
  if (!length(hit)) return(NULL)
  pd[[hit[1L]]]
}
parse_binary <- function(x) {
  y <- tolower(trimws(as.character(x)))
  out <- ifelse(y %in% c("1","yes","positive","pos","lvi+","lvi positive","metastatic","ln+"), 1,
                ifelse(y %in% c("0","no","negative","neg","lvi-","lvi negative","non-metastatic","ln-"), 0, NA))
  as.numeric(out)
}

# ---- 6. GSE62254: actual LVI and survival -------------------------------
eset62254 <- geo[["GSE62254"]]
expr62254 <- eset_to_gene_matrix(eset62254)
pd62254 <- as.data.frame(Biobase::pData(eset62254))
write.csv(pd62254, file.path(ROOT, "tables", "GSE62254_pData_raw.csv"))
# Inspect these columns once; then set exact names if needed.
print(colnames(pd62254))
lvi62254 <- parse_binary(auto_field(pd62254, c("lymphovascular","lvi","vascular invasion")))
if (all(is.na(lvi62254))) warning("GSE62254 LVI column not auto-detected; set lvi62254 manually after inspecting pData.")
sc62254 <- score_programs(expr62254)
common <- complete.cases(lvi62254)
dat62254 <- cbind(pd62254[common, , drop = FALSE], lvi = lvi62254[common], sc62254[common, , drop = FALSE])
# Rename/derive survival columns using the actual pData fields.
os_time <- suppressWarnings(as.numeric(auto_field(dat62254, c("overall.*survival.*time","os.*time","survival.*months","days.*death"))))
os_event <- parse_binary(auto_field(dat62254, c("overall.*survival.*event","os.*event","vital.*status","death")))
dat62254$os_time <- os_time; dat62254$os_event <- os_event

wilcox_tbl <- function(dat, group, vars) {
  purrr::map_dfr(vars, function(v) {
    x <- dat[[v]]; g <- dat[[group]]
    ok <- complete.cases(x, g)
    wt <- wilcox.test(x[ok] ~ g[ok], exact = FALSE)
    tibble(program = v, n = sum(ok), p = wt$p.value,
           median_0 = median(x[ok & g[ok] == 0]), median_1 = median(x[ok & g[ok] == 1]))
  }) |> mutate(q = p.adjust(p, "BH"))
}
program_names <- names(PROGRAMS)
if (any(common)) {
  lvi_tests <- wilcox_tbl(dat62254, "lvi", program_names)
  write.csv(lvi_tests, file.path(ROOT, "tables", "GSE62254_LVI_Wilcoxon.csv"), row.names = FALSE)
  p1 <- dat62254 |> select(lvi, all_of(program_names)) |> pivot_longer(-lvi, names_to = "program", values_to = "score") |>
    ggplot(aes(factor(lvi), score, fill = factor(lvi))) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = .15, alpha = .35) +
    facet_wrap(~program, scales = "free_y") + labs(x = "LVI (0=negative, 1=positive)", y = "Program score") + theme_bw() + theme(legend.position = "none")
  ggsave(file.path(ROOT, "figures", "Fig_GSE62254_LVI_programs.pdf"), p1, width = 12, height = 7)
  # Adjusted logistic models; add age/sex/stage only when fields are identified.
  age <- suppressWarnings(as.numeric(auto_field(dat62254, c("age"))))
  sex <- auto_field(dat62254, c("gender","sex")); stage <- auto_field(dat62254, c("stage"))
  covars <- list(age = scale(age)[,1], sex = factor(sex), stage = factor(stage))
  logit_tbl <- purrr::map_dfr(program_names, function(v) {
    d <- data.frame(lvi = dat62254$lvi, score = as.numeric(scale(dat62254[[v]])), covars, check.names = FALSE)
    d <- d[complete.cases(d), , drop = FALSE]
    fit <- glm(lvi ~ score + age + sex + stage, data = d, family = binomial())
    broom::tidy(fit, conf.int = TRUE, exponentiate = TRUE) |> filter(term == "score") |> mutate(program = v)
  }) |> mutate(q = p.adjust(p.value, "BH"))
  write.csv(logit_tbl, file.path(ROOT, "tables", "GSE62254_LVI_logistic.csv"), row.names = FALSE)
}

# Survival analysis in GSE62254 when time/event columns were found.
if (sum(complete.cases(dat62254$os_time, dat62254$os_event)) > 20) {
  cox_tbl <- purrr::map_dfr(program_names, function(v) {
    d <- data.frame(time = dat62254$os_time, event = dat62254$os_event, score = as.numeric(scale(dat62254[[v]]))) |> na.omit()
    fit <- survival::coxph(Surv(time, event) ~ score, data = d)
    broom::tidy(fit, exponentiate = TRUE, conf.int = TRUE) |> filter(term == "score") |> mutate(program = v)
  }) |> mutate(q = p.adjust(p.value, "BH"))
  write.csv(cox_tbl, file.path(ROOT, "tables", "GSE62254_survival_cox.csv"), row.names = FALSE)
}

# ---- 7. TCGA-STAD from GDC -----------------------------------------------
q <- TCGAbiolinks::GDCquery(project = "TCGA-STAD", data.category = "Transcriptome Profiling",
                            data.type = "Gene Expression Quantification", workflow.type = "STAR - Counts",
                            sample.type = "Primary Tumor")
TCGAbiolinks::GDCdownload(q, directory = file.path(ROOT, "raw", "GDCdata"))
tcga_se <- TCGAbiolinks::GDCprepare(q, directory = file.path(ROOT, "raw", "GDCdata"), summarizedExperiment = TRUE)
assay_name <- if ("fpkm_unstrand" %in% SummarizedExperiment::assayNames(tcga_se)) "fpkm_unstrand" else SummarizedExperiment::assayNames(tcga_se)[1]
tcga_expr <- SummarizedExperiment::assay(tcga_se, assay_name)
if ("gene_name" %in% colnames(rowData(tcga_se))) rownames(tcga_expr) <- rowData(tcga_se)$gene_name
tcga_expr <- log2(t(t(tcga_expr) / pmax(colSums(tcga_expr), 1) * 1e6) + 1)
tcga_scores <- score_programs(tcga_expr)
tcga_clin <- TCGAbiolinks::GDCquery_clinic(project = "TCGA-STAD", type = "clinical")
write.csv(tcga_clin, file.path(ROOT, "tables", "TCGA_STAD_clinical.csv"), row.names = FALSE)
write.csv(tcga_scores, file.path(ROOT, "tables", "TCGA_STAD_program_scores.csv"))
# Correlation heatmap and pair plot.
p_cor <- cor(tcg_scores, use = "pairwise.complete.obs", method = "spearman")
pheatmap::pheatmap(p_cor, filename = file.path(ROOT, "figures", "Fig_TCGA_program_correlation.pdf"), width = 7, height = 6)
pair <- data.frame(fatty_acid_degradation = tcga_scores$fatty_acid_degradation, EMT_ECM = tcga_scores$EMT_ECM)
pair <- pair[complete.cases(pair), ]
p2 <- ggplot(pair, aes(fatty_acid_degradation, EMT_ECM)) + geom_point(alpha = .6) + geom_smooth(method = "lm", se = FALSE) + theme_bw()
ggsave(file.path(ROOT, "figures", "Fig_TCGA_fatty_acid_vs_EMT.pdf"), p2, width = 5, height = 4)

# ---- 8. Surrogate GEO cohorts ---------------------------------------------
score_geo <- function(accession, group_patterns) {
  e <- geo[[accession]]; x <- eset_to_gene_matrix(e); pd <- as.data.frame(Biobase::pData(e))
  g <- parse_binary(auto_field(pd, group_patterns))
  s <- score_programs(x); write.csv(s, file.path(ROOT, "tables", paste0(accession, "_program_scores.csv")))
  if (!all(is.na(g))) {
    d <- cbind(group = g, s) |> as.data.frame() |> pivot_longer(-group, names_to = "program", values_to = "score")
    p <- ggplot(d, aes(factor(group), score, fill = factor(group))) + geom_boxplot() + facet_wrap(~program, scales = "free_y") + theme_bw() + theme(legend.position = "none")
    ggsave(file.path(ROOT, "figures", paste0("Fig_", accession, "_programs.pdf")), p, width = 12, height = 7)
  }
  invisible(list(expr = x, pData = pd, group = g, scores = s))
}
g246963 <- score_geo("GSE246963", c("lymph.node", "ln.status", "metastasis"))
g237876 <- score_geo("GSE237876", c("site", "metastatic", "primary", "route"))

# ---- 9. GSE261635 EPIC methylation ----------------------------------------
# GEOquery downloads the matrix and sample annotations; probe-to-promoter
# annotation must be supplied from the EPIC manifest used by the study.
e261635 <- geo[["GSE261635"]]
meth <- Biobase::exprs(e261635)
pd261635 <- as.data.frame(Biobase::pData(e261635))
write.csv(pd261635, file.path(ROOT, "tables", "GSE261635_pData_raw.csv"))
write.csv(as.data.frame(meth), file.path(ROOT, "tables", "GSE261635_beta_matrix.csv"))
# Required user-supplied file: EPIC_manifest_promoters.csv with columns probe,gene.
manifest_file <- file.path(ROOT, "raw", "EPIC_manifest_promoters.csv")
if (file.exists(manifest_file)) {
  man <- fread(manifest_file)
  stopifnot(all(c("probe","gene") %in% names(man)))
  man <- man[probe %in% rownames(meth) & !is.na(gene)]
  gene_beta <- sapply(split(man$probe, man$gene), function(p) colMeans(meth[p, , drop = FALSE], na.rm = TRUE))
  gene_beta <- as.matrix(gene_beta); gene_beta <- t(gene_beta)
  rownames(gene_beta) <- colnames(meth)
  meth_scores <- score_programs(-t(gene_beta)) # negative z = putative activation proxy
  write.csv(meth_scores, file.path(ROOT, "tables", "GSE261635_program_scores.csv"))
  lvi_m <- parse_binary(auto_field(pd261635, c("lvi","lymphovascular","vascular invasion")))
  if (!all(is.na(lvi_m))) {
    md <- cbind(lvi = lvi_m, meth_scores) |> as.data.frame() |> pivot_longer(-lvi, names_to = "program", values_to = "score")
    pm <- ggplot(md, aes(factor(lvi), score, fill = factor(lvi))) + geom_boxplot() + facet_wrap(~program, scales = "free_y") + theme_bw() + theme(legend.position = "none")
    ggsave(file.path(ROOT, "figures", "Fig_GSE261635_methylation_programs.pdf"), pm, width = 12, height = 7)
  }
} else message("Provide ", manifest_file, " to run promoter-level methylation scoring.")

# ---- 10. Single-cell GEO datasets: download and save raw objects ---------
# These datasets are intentionally not treated as independent cells. Create
# sample x lineage pseudobulk after applying the manuscript's QC and lineage
# marker rules; the exact matrix format differs by GEO release.
for (id in c("GSE163558", "GSE239676")) {
  e <- geo[[id]]
  saveRDS(list(exprs = Biobase::exprs(e), pData = Biobase::pData(e), fData = Biobase::fData(e)),
          file.path(ROOT, "rds", paste0(id, "_raw_components.rds")))
}

# ---- 11. Optional DepMap and LINCS prioritization --------------------------
# Download the current public files from DepMap/LINCS portals manually or via
# their official URLs, then set these paths. No causal claim is made here.
depmap_file <- Sys.getenv("DEPMAP_CRISPR_FILE", unset = "")
if (nzchar(depmap_file) && file.exists(depmap_file)) {
  dm <- fread(depmap_file)
  # Adapt column names to the downloaded release; expected columns are gene,
  # lineage/model and gene-effect score.
  print(names(dm))
  fwrite(dm, file.path(ROOT, "tables", "DepMap_input_snapshot.csv"))
}
lincs_file <- Sys.getenv("LINCS_SIGNATURE_FILE", unset = "")
if (nzchar(lincs_file) && file.exists(lincs_file)) {
  lc <- fread(lincs_file); print(names(lc)); fwrite(lc, file.path(ROOT, "tables", "LINCS_input_snapshot.csv"))
}

# ---- 12. Session information ----------------------------------------------
writeLines(capture.output(sessionInfo()), file.path(ROOT, "sessionInfo.txt"))
message("Finished. Tables and figures are in: ", normalizePath(ROOT))