#!/usr/bin/env Rscript

# ============================================================
# CRC diabetes-related metabolic phenotype scoring
# Manual weighted ssGSEA-like implementation without GSVA
# Author: Dr Henry Okuchukwu Ebili (MBBS,PhD,FMCPath); ebili.henry@oouagoiwoye.edu.ng 
# Primary phenotype score:
#   Z(ssGSEA_up) - Z(ssGSEA_down)
#
# Optional:
#   residualize phenotype scores against covariates from metadata
#
# Usage:
#   ./crc_dm_fixed_ssgsea.R expression_matrix.tsv output_prefix
#
#   ./crc_dm_fixed_ssgsea.R expression_matrix.tsv output_prefix metadata.tsv purity,stromal,batch
#
# Input expression matrix:
#   - tab-delimited
#   - first column = gene symbols
#   - remaining columns = samples
#   - values should be normalized continuous expression
#     (e.g. log2(TPM+1), log2(FPKM+1), vst/rlog)
#
# Metadata file (optional):
#   - tab-delimited
#   - must contain a column named SampleID
#   - additional covariate columns can be numeric or factor-like
# ============================================================

args <- commandArgs(trailingOnly = TRUE)

if (length(args) < 2) {
  stop(
    paste(
      "Usage:",
      "  ./crc_dm_fixed_ssgsea.R <expression_matrix.tsv> <output_prefix> [metadata.tsv] [covariates_comma_sep]",
      sep = "\n"
    )
  )
}

expr_file <- args[1]
out_prefix <- args[2]
metadata_file <- if (length(args) >= 3) args[3] else NA_character_
covariates_arg <- if (length(args) >= 4) args[4] else NA_character_

# -------------------------
# 1. Utility functions
# -------------------------
message_line <- function(...) cat(..., "\n", sep = "")

safe_dir_create <- function(path) {
  d <- dirname(path)
  if (!dir.exists(d) && d != ".") dir.create(d, recursive = TRUE)
}

zscore_vec <- function(x) {
  s <- sd(x, na.rm = TRUE)
  m <- mean(x, na.rm = TRUE)
  if (is.na(s) || s == 0) {
    return(rep(NA_real_, length(x)))
  }
  (x - m) / s
}

write_table_safe <- function(x, file, row.names = FALSE, col.names = TRUE) {
  write.table(
    x, file = file, sep = "\t", quote = FALSE,
    row.names = row.names, col.names = col.names
  )
}

clean_gene_symbols <- function(x) {
  x <- trimws(x)
  x <- sub("\\.\\d+$", "", x)     # remove Ensembl-style version suffixes
  x <- toupper(x)
  x
}

looks_like_raw_counts <- function(mat) {
  if (length(mat) == 0) return(FALSE)
  finite_vals <- mat[is.finite(mat)]
  if (length(finite_vals) == 0) return(FALSE)
  frac_integer <- mean(abs(finite_vals - round(finite_vals)) < 1e-8)
  max_val <- max(finite_vals, na.rm = TRUE)
  med_val <- median(finite_vals, na.rm = TRUE)
  frac_zero <- mean(finite_vals == 0)
  frac_integer > 0.95 && max_val > 50 && med_val >= 0 && frac_zero > 0.05
}

collapse_duplicate_genes <- function(expr_matrix) {
  if (!anyDuplicated(rownames(expr_matrix))) return(expr_matrix)
  message_line("Duplicate gene symbols detected. Collapsing duplicates by mean expression.")
  idx_split <- split(seq_len(nrow(expr_matrix)), rownames(expr_matrix))
  collapsed <- t(vapply(
    idx_split,
    function(idx) {
      if (length(idx) == 1) {
        expr_matrix[idx, ]
      } else {
        colMeans(expr_matrix[idx, , drop = FALSE], na.rm = TRUE)
      }
    },
    FUN.VALUE = numeric(ncol(expr_matrix))
  ))
  colnames(collapsed) <- colnames(expr_matrix)
  rownames(collapsed) <- names(idx_split)
  mode(collapsed) <- "numeric"
  collapsed
}

# -------------------------
# 2. Read expression matrix
# -------------------------
message_line("Reading expression matrix: ", expr_file)

expr_df <- read.table(
  expr_file,
  header = TRUE,
  sep = "\t",
  check.names = FALSE,
  stringsAsFactors = FALSE,
  quote = "",
  comment.char = ""
)

if (ncol(expr_df) < 2) {
  stop("Expression file must have at least 2 columns: gene symbol + >=1 sample columns.")
}

gene_symbols <- clean_gene_symbols(expr_df[[1]])
if (any(is.na(gene_symbols)) || any(gene_symbols == "")) {
  stop("Gene symbol column contains missing or empty values.")
}

expr <- as.matrix(expr_df[, -1, drop = FALSE])
mode(expr) <- "numeric"
rownames(expr) <- gene_symbols

if (anyNA(expr)) {
  message_line("Warning: expression matrix contains NA values. They will be tolerated where possible.")
}

# remove rows with all NA
all_na_rows <- rowSums(is.na(expr)) == ncol(expr)
if (any(all_na_rows)) {
  message_line("Removing ", sum(all_na_rows), " genes with all-NA values.")
  expr <- expr[!all_na_rows, , drop = FALSE]
}

# collapse duplicate genes
expr <- collapse_duplicate_genes(expr)

message_line(sprintf("Loaded matrix: %d genes x %d samples", nrow(expr), ncol(expr)))

if (looks_like_raw_counts(expr)) {
  warning(
    "Input looks like raw counts. This script expects normalized continuous expression ",
    "(e.g. log2(TPM+1), log2(FPKM+1), vst/rlog). Results may be misleading if raw counts are used."
  )
}

if (all(grepl("^ENSG", rownames(expr)))) {
  warning(
    "Row names look like Ensembl IDs, but this script expects gene symbols. ",
    "Most module genes will not match unless you convert IDs first."
  )
}

# -------------------------
# 3. Define modules
# -------------------------
modules <- list(
  Hyperglycaemia_up = c(
    "HK2","PFKP","PKM","PCK1","G6PC","SLC2A1","LDHA","ALDOA","ENO1"
  ),
  Hyperglycaemia_down = c(
    "SLC2A4","INSR","IRS1","IRS2","AKT2"
  ),

  InsulinResistance_up = c(
    "FOXO1","PCK1","G6PC","PTPN1","PTEN","SOCS1","SOCS3",
    "TRIB3","GSK3B","TNF","IL6","STAT3","NFKB1","JUN"
  ),
  InsulinResistance_down = c(
    "INSR","IRS1","IRS2","PIK3CA","PIK3R1","PDPK1",
    "AKT1","AKT2","MTOR","SLC2A4"
  ),

  LipidMetabolism_up = c(
    "FASN","ACACA","ACLY","SCD","ELOVL6","SREBF1",
    "DGAT1","DGAT2","FABP4","CD36","CPT1A","ACADM","ACOX1"
  ),
  LipidMetabolism_down = c(
    "PPARA","PPARGC1A","LIPE","PNPLA2","ADIPOQ"
  ),

  MitoDysfunction_up = c(
    "HIF1A","SOD2","GPX1","HMOX1","ATF4","DDIT3",
    "XBP1","BNIP3","PDK1","NQO1"
  ),
  MitoDysfunction_down = c(
    "NDUFA1","NDUFA9","NDUFS1","SDHA","UQCRC1",
    "COX4I1","COX5A","ATP5F1A","ATP5F1B","TFAM","PPARGC1A"
  ),

  MetabolicSyndrome_up = c(
    "IL6","TNF","CCL2","NFKB1","RELA","STAT3","SOCS3","CRP",
    "LEP","SREBF1","FASN","PCK1","G6PC","SERPINE1","ICAM1","VCAM1"
  ),
  MetabolicSyndrome_down = c(
    "ADIPOQ","INSR","IRS1","IRS2","SLC2A4","PPARA","PPARGC1A"
  )
)

modules <- lapply(modules, function(x) unique(clean_gene_symbols(x)))
all_signature_genes <- unique(unlist(modules))

# -------------------------
# 4. Filter low-expression genes
#    BUT ALWAYS KEEP MODULE GENES
# -------------------------
mean_expr <- rowMeans(expr, na.rm = TRUE)
threshold <- quantile(mean_expr, 0.10, na.rm = TRUE)

keep_low_filter <- mean_expr > threshold
keep_signature <- rownames(expr) %in% all_signature_genes
keep <- keep_low_filter | keep_signature

expr_filtered <- expr[keep, , drop = FALSE]

message_line(sprintf(
  "After filtering: %d genes retained (%d removed). All signature genes were protected from filtering.",
  nrow(expr_filtered), nrow(expr) - nrow(expr_filtered)
))

# -------------------------
# 5. Module coverage report
# -------------------------
module_coverage <- data.frame(
  Module = names(modules),
  OriginalSize = vapply(modules, length, integer(1)),
  DetectedSize = integer(length(modules)),
  CoveragePct = numeric(length(modules)),
  MissingGenes = character(length(modules)),
  stringsAsFactors = FALSE
)

for (i in seq_along(modules)) {
  mod <- modules[[i]]
  present <- intersect(mod, rownames(expr_filtered))
  missing <- setdiff(mod, rownames(expr_filtered))
  module_coverage$DetectedSize[i] <- length(present)
  module_coverage$CoveragePct[i] <- 100 * length(present) / length(mod)
  module_coverage$MissingGenes[i] <- paste(missing, collapse = ",")
}

message_line("\nModule coverage:")
print(module_coverage[, c("Module", "OriginalSize", "DetectedSize", "CoveragePct")], row.names = FALSE)

# stricter minimum usable size
min_genes_abs <- 3
min_genes_prop <- 0.50

module_pass <- vapply(seq_along(modules), function(i) {
  det <- module_coverage$DetectedSize[i]
  orig <- module_coverage$OriginalSize[i]
  det >= min_genes_abs && det / orig >= min_genes_prop
}, logical(1))

usable_modules <- modules[module_pass]
failed_modules <- setdiff(names(modules), names(usable_modules))

if (length(failed_modules) > 0) {
  warning(
    "The following modules failed minimum coverage criteria and will be returned as NA: ",
    paste(failed_modules, collapse = ", ")
  )
}

# -------------------------
# 6. Corrected manual weighted ssGSEA-like scoring
#    FIX: weights use expression magnitude, not rank number
# -------------------------
calculate_weighted_ssgsea <- function(expr_matrix, gene_sets, alpha = 0.25) {
  samples <- colnames(expr_matrix)
  set_names <- names(gene_sets)

  score_mat <- matrix(
    NA_real_,
    nrow = length(samples),
    ncol = length(set_names),
    dimnames = list(samples, set_names)
  )

  for (s in seq_along(samples)) {
    sample_name <- samples[s]
    x <- expr_matrix[, s]
    names(x) <- rownames(expr_matrix)

    # drop NA genes for this sample
    x <- x[!is.na(x)]
    if (length(x) == 0) next

    # sort genes by decreasing expression
    ord <- order(x, decreasing = TRUE)
    x_sorted <- x[ord]
    genes_sorted <- names(x_sorted)
    N <- length(x_sorted)

    # transform expression to non-negative weights
    # preserves within-sample ranking, prevents negative weight issues
    x_shifted <- x_sorted - min(x_sorted, na.rm = TRUE) + 1e-12

    for (g in seq_along(gene_sets)) {
      gs_name <- set_names[g]
      gs <- intersect(gene_sets[[g]], genes_sorted)
      Nh <- length(gs)

      if (Nh == 0) {
        score_mat[s, g] <- NA_real_
        next
      }

      hit_idx <- genes_sorted %in% gs
      miss_idx <- !hit_idx
      Nm <- sum(miss_idx)

      if (Nm == 0) {
        score_mat[s, g] <- NA_real_
        next
      }

      # weighted hit contribution
      hit_weights <- abs(x_shifted[hit_idx])^alpha
      hit_weights_sum <- sum(hit_weights)

      if (!is.finite(hit_weights_sum) || hit_weights_sum <= 0) {
        score_mat[s, g] <- NA_real_
        next
      }

      P_hit <- rep(0, N)
      P_hit[hit_idx] <- hit_weights / hit_weights_sum

      P_miss <- rep(1 / Nm, N)
      P_miss[hit_idx] <- 0

      running_es <- cumsum(P_hit) - cumsum(P_miss)

      # signed maximum deviation
      max_i <- which.max(abs(running_es))
      es <- running_es[max_i]

      # mild normalization by set size
      es <- es / sqrt(Nh)
      score_mat[s, g] <- es
    }
  }

  score_mat
}

message_line("\nCalculating raw module scores...")
module_scores_raw <- calculate_weighted_ssgsea(expr_filtered, usable_modules, alpha = 0.25)

# Rebuild full module matrix including failed modules as NA
full_module_scores_raw <- matrix(
  NA_real_,
  nrow = nrow(module_scores_raw),
  ncol = length(modules),
  dimnames = list(rownames(module_scores_raw), names(modules))
)
full_module_scores_raw[, colnames(module_scores_raw)] <- module_scores_raw
module_scores_raw <- full_module_scores_raw

# -------------------------
# 7. Z-standardize module scores across samples
# -------------------------
module_scores_z <- apply(module_scores_raw, 2, zscore_vec)
if (is.vector(module_scores_z)) {
  module_scores_z <- matrix(module_scores_z, ncol = 1)
  colnames(module_scores_z) <- colnames(module_scores_raw)[1]
}
rownames(module_scores_z) <- rownames(module_scores_raw)

zero_var_modules <- colnames(module_scores_raw)[apply(module_scores_raw, 2, function(x) {
  s <- sd(x, na.rm = TRUE)
  is.na(s) || s == 0
})]

if (length(zero_var_modules) > 0) {
  warning(
    "These modules had zero/undefined variance across samples after raw scoring and became NA after Z-standardization: ",
    paste(zero_var_modules, collapse = ", ")
  )
}

# -------------------------
# 8. Combine up/down modules
# -------------------------
safe_col <- function(mat, nm) {
  if (nm %in% colnames(mat)) mat[, nm] else rep(NA_real_, nrow(mat))
}

combine_phenotypes <- function(module_mat) {
  out <- data.frame(
    SampleID = rownames(module_mat),
    Hyperglycaemia =
      safe_col(module_mat, "Hyperglycaemia_up") -
      safe_col(module_mat, "Hyperglycaemia_down"),
    InsulinResistance =
      safe_col(module_mat, "InsulinResistance_up") -
      safe_col(module_mat, "InsulinResistance_down"),
    LipidMetabolism =
      safe_col(module_mat, "LipidMetabolism_up") -
      safe_col(module_mat, "LipidMetabolism_down"),
    MitoDysfunction =
      safe_col(module_mat, "MitoDysfunction_up") -
      safe_col(module_mat, "MitoDysfunction_down"),
    MetabolicSyndrome =
      safe_col(module_mat, "MetabolicSyndrome_up") -
      safe_col(module_mat, "MetabolicSyndrome_down"),
    stringsAsFactors = FALSE,
    check.names = FALSE
  )
  rownames(out) <- out$SampleID
  out
}

phenotype_scores_raw <- combine_phenotypes(module_scores_raw)
phenotype_scores_z <- combine_phenotypes(module_scores_z)

# -------------------------
# 9. Correlation with p-values and FDR
# -------------------------
pairwise_cor_test <- function(score_df, method = "spearman") {
  mat <- as.matrix(score_df[, -1, drop = FALSE])
  phenos <- colnames(mat)

  cor_mat <- matrix(NA_real_, nrow = ncol(mat), ncol = ncol(mat),
                    dimnames = list(phenos, phenos))
  p_mat <- cor_mat
  n_mat <- cor_mat

  for (i in seq_along(phenos)) {
    for (j in seq_along(phenos)) {
      x <- mat[, i]
      y <- mat[, j]
      ok <- complete.cases(x, y)
      n_use <- sum(ok)

      n_mat[i, j] <- n_use

      if (i == j) {
        cor_mat[i, j] <- 1
        p_mat[i, j] <- 0
      } else if (n_use >= 3) {
        test <- suppressWarnings(cor.test(x[ok], y[ok], method = method, exact = FALSE))
        cor_mat[i, j] <- unname(test$estimate)
        p_mat[i, j] <- test$p.value
      }
    }
  }

  # upper triangle FDR
  ut <- upper.tri(p_mat)
  fdr_vals <- p.adjust(p_mat[ut], method = "BH")
  fdr_mat <- matrix(NA_real_, nrow = nrow(p_mat), ncol = ncol(p_mat),
                    dimnames = dimnames(p_mat))
  fdr_mat[ut] <- fdr_vals
  fdr_mat[lower.tri(fdr_mat)] <- t(fdr_mat)[lower.tri(fdr_mat)]
  diag(fdr_mat) <- 0

  list(cor = cor_mat, p = p_mat, fdr = fdr_mat, n = n_mat)
}

flatten_pairwise_results <- function(cor_obj) {
  phenos <- rownames(cor_obj$cor)
  out <- data.frame(
    Phenotype1 = character(),
    Phenotype2 = character(),
    Rho = numeric(),
    Pvalue = numeric(),
    FDR = numeric(),
    N = integer(),
    stringsAsFactors = FALSE
  )

  for (i in seq_len(length(phenos))) {
    for (j in seq_len(length(phenos))) {
      if (i < j) {
        out <- rbind(out, data.frame(
          Phenotype1 = phenos[i],
          Phenotype2 = phenos[j],
          Rho = cor_obj$cor[i, j],
          Pvalue = cor_obj$p[i, j],
          FDR = cor_obj$fdr[i, j],
          N = cor_obj$n[i, j],
          stringsAsFactors = FALSE
        ))
      }
    }
  }

  out[order(out$Pvalue), ]
}

# -------------------------
# 10. Optional confounder adjustment
# -------------------------
residualize_scores <- function(score_df, metadata, covariates) {
  stopifnot("SampleID" %in% colnames(metadata))

  merged <- merge(
    data.frame(SampleID = score_df$SampleID, score_df[, -1, drop = FALSE], check.names = FALSE),
    metadata,
    by = "SampleID",
    all.x = FALSE,
    all.y = FALSE
  )

  if (nrow(merged) == 0) {
    stop("No overlapping samples between phenotype score table and metadata file.")
  }

  for (cv in covariates) {
    if (!(cv %in% colnames(merged))) {
      stop("Covariate not found in metadata: ", cv)
    }
    # try numeric conversion, otherwise factor
    suppressWarnings(num_try <- as.numeric(merged[[cv]]))
    if (!all(is.na(num_try))) {
      # keep numeric only if at least some conversion succeeded meaningfully
      if (sum(!is.na(num_try)) >= ceiling(0.8 * length(num_try))) {
        merged[[cv]] <- num_try
      } else {
        merged[[cv]] <- as.factor(merged[[cv]])
      }
    } else {
      merged[[cv]] <- as.factor(merged[[cv]])
    }
  }

  residual_df <- data.frame(SampleID = merged$SampleID, check.names = FALSE)

  pheno_cols <- setdiff(colnames(score_df), "SampleID")
  formula_rhs <- paste(covariates, collapse = " + ")

  for (ph in pheno_cols) {
    dat <- merged[, c("SampleID", ph, covariates), drop = FALSE]
    dat <- dat[complete.cases(dat), , drop = FALSE]

    if (nrow(dat) < 5) {
      residual_df[[ph]] <- NA_real_
      next
    }

    form <- as.formula(paste(ph, "~", formula_rhs))
    fit <- try(lm(form, data = dat), silent = TRUE)

    res <- rep(NA_real_, nrow(merged))
    if (!inherits(fit, "try-error")) {
      res_match <- residuals(fit)
      idx <- match(dat$SampleID, merged$SampleID)
      res[idx] <- res_match
    }

    residual_df[[ph]] <- res
  }

  rownames(residual_df) <- residual_df$SampleID
  residual_df
}

cor_raw <- pairwise_cor_test(phenotype_scores_raw, method = "spearman")
cor_z <- pairwise_cor_test(phenotype_scores_z, method = "spearman")

residualized_scores <- NULL
cor_resid <- NULL

if (!is.na(metadata_file) && file.exists(metadata_file) && !is.na(covariates_arg)) {
  message_line("\nReading metadata: ", metadata_file)
  metadata <- read.table(
    metadata_file,
    header = TRUE,
    sep = "\t",
    check.names = FALSE,
    stringsAsFactors = FALSE,
    quote = "",
    comment.char = ""
  )

  if (!("SampleID" %in% colnames(metadata))) {
    stop("Metadata file must contain a column named 'SampleID'.")
  }

  covariates <- trimws(unlist(strsplit(covariates_arg, ",")))
  covariates <- covariates[covariates != ""]

  if (length(covariates) > 0) {
    message_line("Residualizing Z-phenotype scores against covariates: ", paste(covariates, collapse = ", "))
    residualized_scores <- residualize_scores(phenotype_scores_z, metadata, covariates)
    cor_resid <- pairwise_cor_test(residualized_scores, method = "spearman")
  }
}

# -------------------------
# 11. Output summaries
# -------------------------
summarize_score_matrix <- function(mat, name) {
  message_line("\n", name)
  print(data.frame(
    Score = colnames(mat),
    Mean = apply(mat, 2, mean, na.rm = TRUE),
    SD = apply(mat, 2, sd, na.rm = TRUE),
    Min = apply(mat, 2, min, na.rm = TRUE),
    Median = apply(mat, 2, median, na.rm = TRUE),
    Max = apply(mat, 2, max, na.rm = TRUE),
    Missing = apply(mat, 2, function(x) sum(is.na(x))),
    row.names = NULL
  ))
}

summarize_score_matrix(module_scores_raw, "Raw module scores")
summarize_score_matrix(module_scores_z, "Z-standardized module scores")
summarize_score_matrix(as.matrix(phenotype_scores_z[, -1, drop = FALSE]), "Primary phenotype scores (Z-up minus Z-down)")

# -------------------------
# 12. Export files
# -------------------------
safe_dir_create(out_prefix)

write_table_safe(
  cbind(Gene = rownames(expr_filtered), expr_filtered),
  paste0(out_prefix, "_filtered_expression.tsv"),
  row.names = FALSE
)

write_table_safe(module_coverage, paste0(out_prefix, "_module_coverage.tsv"), row.names = FALSE)
write_table_safe(
  cbind(SampleID = rownames(module_scores_raw), module_scores_raw),
  paste0(out_prefix, "_module_scores_raw.tsv"),
  row.names = FALSE
)
write_table_safe(
  cbind(SampleID = rownames(module_scores_z), module_scores_z),
  paste0(out_prefix, "_module_scores_z.tsv"),
  row.names = FALSE
)
write_table_safe(phenotype_scores_raw, paste0(out_prefix, "_phenotype_scores_raw.tsv"), row.names = FALSE)
write_table_safe(phenotype_scores_z, paste0(out_prefix, "_phenotype_scores_z.tsv"), row.names = FALSE)

write_table_safe(cor_raw$cor, paste0(out_prefix, "_correlation_matrix_raw.tsv"), row.names = TRUE, col.names = NA)
write_table_safe(cor_raw$p, paste0(out_prefix, "_correlation_pvalues_raw.tsv"), row.names = TRUE, col.names = NA)
write_table_safe(cor_raw$fdr, paste0(out_prefix, "_correlation_fdr_raw.tsv"), row.names = TRUE, col.names = NA)
write_table_safe(flatten_pairwise_results(cor_raw), paste0(out_prefix, "_pairwise_correlations_raw.tsv"), row.names = FALSE)

write_table_safe(cor_z$cor, paste0(out_prefix, "_correlation_matrix_z.tsv"), row.names = TRUE, col.names = NA)
write_table_safe(cor_z$p, paste0(out_prefix, "_correlation_pvalues_z.tsv"), row.names = TRUE, col.names = NA)
write_table_safe(cor_z$fdr, paste0(out_prefix, "_correlation_fdr_z.tsv"), row.names = TRUE, col.names = NA)
write_table_safe(flatten_pairwise_results(cor_z), paste0(out_prefix, "_pairwise_correlations_z.tsv"), row.names = FALSE)

if (!is.null(residualized_scores)) {
  write_table_safe(residualized_scores, paste0(out_prefix, "_phenotype_scores_residualized.tsv"), row.names = FALSE)
}
if (!is.null(cor_resid)) {
  write_table_safe(cor_resid$cor, paste0(out_prefix, "_correlation_matrix_residualized.tsv"), row.names = TRUE, col.names = NA)
  write_table_safe(cor_resid$p, paste0(out_prefix, "_correlation_pvalues_residualized.tsv"), row.names = TRUE, col.names = NA)
  write_table_safe(cor_resid$fdr, paste0(out_prefix, "_correlation_fdr_residualized.tsv"), row.names = TRUE, col.names = NA)
  write_table_safe(flatten_pairwise_results(cor_resid), paste0(out_prefix, "_pairwise_correlations_residualized.tsv"), row.names = FALSE)
}

# module gene membership file
module_genes_long <- data.frame(
  Module = rep(names(modules), vapply(modules, length, integer(1))),
  Gene = unlist(modules),
  PresentAfterFiltering = unlist(modules) %in% rownames(expr_filtered),
  stringsAsFactors = FALSE
)
write_table_safe(module_genes_long, paste0(out_prefix, "_module_genes.tsv"), row.names = FALSE)

# analysis report
report_file <- paste0(out_prefix, "_analysis_report.txt")
sink(report_file)

cat("CRC DIABETES-RELATED METABOLIC PHENOTYPE ANALYSIS REPORT\n")
cat("========================================================\n\n")
cat("Date: ", as.character(Sys.Date()), "\n", sep = "")
cat("Time: ", as.character(Sys.time()), "\n\n", sep = "")

cat("INPUT SUMMARY\n")
cat("-------------\n")
cat(sprintf("Expression matrix: %d genes x %d samples after filtering\n", nrow(expr_filtered), ncol(expr_filtered)))
cat(sprintf("Total defined modules: %d\n", length(modules)))
cat(sprintf("Usable modules: %d\n", length(usable_modules)))
if (length(failed_modules) > 0) {
  cat("Failed modules: ", paste(failed_modules, collapse = ", "), "\n", sep = "")
}
cat("\n")

cat("PRIMARY SCORING STRATEGY\n")
cat("------------------------\n")
cat("Raw module scores were computed using a corrected weighted ssGSEA-like running-sum statistic.\n")
cat("Hit weights were based on within-sample expression magnitude after ranking by decreasing expression,\n")
cat("not on rank number itself.\n\n")
cat("Primary phenotype score = Z(module_up) - Z(module_down)\n")
cat("Sensitivity phenotype score = raw(module_up) - raw(module_down)\n\n")

cat("WHY THIS VERSION IS SAFER\n")
cat("-------------------------\n")
cat("1. Signature genes were protected from low-expression filtering.\n")
cat("2. Modules with inadequate gene coverage were flagged.\n")
cat("3. Modules were Z-standardized across samples before primary combination.\n")
cat("4. Pairwise phenotype correlations include p-values and FDR.\n")
cat("5. Optional residualized correlations can adjust for covariates.\n\n")

cat("MODULE COVERAGE\n")
cat("---------------\n")
print(module_coverage, row.names = FALSE)
cat("\n")

cat("PAIRWISE CORRELATIONS: PRIMARY Z-BASED PHENOTYPES\n")
cat("-----------------------------------------------\n")
print(flatten_pairwise_results(cor_z), row.names = FALSE)
cat("\n")

if (!is.null(cor_resid)) {
  cat("PAIRWISE CORRELATIONS: RESIDUALIZED Z-BASED PHENOTYPES\n")
  cat("------------------------------------------------------\n")
  print(flatten_pairwise_results(cor_resid), row.names = FALSE)
  cat("\n")
}

sink()

save(
  expr_filtered,
  modules,
  module_coverage,
  usable_modules,
  failed_modules,
  module_scores_raw,
  module_scores_z,
  phenotype_scores_raw,
  phenotype_scores_z,
  cor_raw,
  cor_z,
  residualized_scores,
  cor_resid,
  file = paste0(out_prefix, "_results.RData")
)

message_line("\nDone.")
message_line("Primary output for manuscript use:")
message_line("  ", paste0(out_prefix, "_phenotype_scores_z.tsv"))
message_line("  ", paste0(out_prefix, "_pairwise_correlations_z.tsv"))

if (!is.null(cor_resid)) {
  message_line("Confounder-adjusted output:")
  message_line("  ", paste0(out_prefix, "_phenotype_scores_residualized.tsv"))
  message_line("  ", paste0(out_prefix, "_pairwise_correlations_residualized.tsv"))
}