# ============================================================================
# MULTIPLICITY IN HIGH-DIMENSIONAL DATA - TCGA-BRCA
# Public analysis and figure generation
# ============================================================================
# Developer: Planimeter Kft.
# Version:   1.0.0
# Date:      2026-09-23
#
# Copyright (C) 2026 Planimeter Kft.
#
# License:
#   This program is free software: you can redistribute it and/or modify it
#   under the terms of the GNU General Public License as published by the
#   Free Software Foundation, either version 3 of the License, or (at your
#   option) any later version.
#
#   This program is distributed in the hope that it will be useful, but
#   WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
#   or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
#   for more details.
#
#   SPDX-License-Identifier: GPL-3.0-or-later
#
# AI use statement:
#   Generative AI was used as a support tool during code development, review
#   and validation. The analytical design, final code, interpretation and
#   responsibility for the released materials remain with Planimeter Kft.
#
# Purpose:
#   Reproduce the public Planimeter multiplicity example from compact,
#   analysis-ready TCGA-BRCA-derived CSV files. The main analysis demonstrates
#   how multiple-testing adjustment, empirical Bayes moderation and ranking by
#   effect size versus statistical evidence can change the evidence set.
#
# Required input files (relative paths):
#   data/sample_metadata.csv
#   data/independent_testing_results.csv
#   data/limma_results.csv
#   data/multiplicity_summary.csv
#
# Main outputs:
#   figures/01_multiplicity_narrowing_corridor.png
#   figures/02_empirical_bayes_signal_field.png
#   figures/03_significance_vs_effect_size_scatter.png
#   figures/04_top100_overlap_panel.png
#   figures/05_from_discovery_to_evidence_journey.png
#   output/multiplicity_summary_recalculated.csv
#   output/validation_checks.csv
#   output/sample_structure.csv
#   output/top100_overlap_features.csv
#
# Important methodological note:
#   The distributed CSV files reproduce the original exploratory demonstration.
#   The independent tests were generated by feature-wise t-tests on raw STAR
#   counts, while the limma analysis used log2(count + 1). This is deliberately
#   preserved for exact reproducibility; it is not presented as a recommended
#   RNA-seq differential-expression workflow.
# ============================================================================

# ---- 0. Packages -----------------------------------------------------------

required_packages <- c("dplyr", "ggplot2", "readr", "scales")
missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]
if (length(missing_packages) > 0) {
  stop(
    paste0(
      "Missing required package(s): ",
      paste(missing_packages, collapse = ", "),
      ". Install them before running this script."
    )
  )
}

library(dplyr)
library(ggplot2)
library(readr)
library(scales)

# ---- 1. Relative paths -----------------------------------------------------

DATA_DIR   <- "data"
FIGURE_DIR <- "figures"
OUTPUT_DIR <- "output"

dir.create(FIGURE_DIR, recursive = TRUE, showWarnings = FALSE)
dir.create(OUTPUT_DIR, recursive = TRUE, showWarnings = FALSE)

files <- c(
  metadata = file.path(DATA_DIR, "sample_metadata.csv"),
  independent = file.path(DATA_DIR, "independent_testing_results.csv"),
  limma = file.path(DATA_DIR, "limma_results.csv"),
  supplied_summary = file.path(DATA_DIR, "multiplicity_summary.csv")
)

if (!all(file.exists(files))) {
  stop("One or more required files are missing from data/.")
}

# ---- 2. Read released datasets --------------------------------------------

sample_metadata <- read_csv(files[["metadata"]], show_col_types = FALSE)
independent_results <- read_csv(files[["independent"]], show_col_types = FALSE)
limma_results <- read_csv(files[["limma"]], show_col_types = FALSE)
supplied_summary <- read_csv(files[["supplied_summary"]], show_col_types = FALSE)

required_meta <- c("sample_id", "ajcc_pathologic_stage", "group")
required_ind <- c(
  "feature_id", "effect_log2_mean_difference", "p_value_raw",
  "p_value_bh", "p_value_bonferroni"
)
required_limma <- c(
  "feature_id", "logFC", "AveExpr", "t", "P.Value", "adj.P.Val", "B"
)

stopifnot(all(required_meta %in% names(sample_metadata)))
stopifnot(all(required_ind %in% names(independent_results)))
stopifnot(all(required_limma %in% names(limma_results)))

# ---- 3. Recalculate key results -------------------------------------------

n_all   <- nrow(independent_results)
n_raw   <- sum(independent_results$p_value_raw < 0.05, na.rm = TRUE)
n_bh    <- sum(independent_results$p_value_bh < 0.05, na.rm = TRUE)
n_bonf  <- sum(independent_results$p_value_bonferroni < 0.05, na.rm = TRUE)
n_limma <- sum(limma_results$adj.P.Val < 0.05, na.rm = TRUE)

limma_plot <- limma_results %>%
  mutate(
    absFC = abs(logFC),
    neglogFDR = -log10(adj.P.Val)
  ) %>%
  filter(is.finite(absFC), is.finite(neglogFDR))

top100_p <- limma_plot %>%
  arrange(adj.P.Val) %>%
  slice_head(n = 100)

top100_fc <- limma_plot %>%
  arrange(desc(absFC)) %>%
  slice_head(n = 100)

common <- intersect(top100_p$feature_id, top100_fc$feature_id)

cor_effect_sig <- cor(
  limma_plot$absFC,
  limma_plot$neglogFDR,
  use = "complete.obs"
)

summary_recalculated <- tibble(
  metric = c(
    "Samples retained",
    "Early-stage samples",
    "Late-stage samples",
    "All tested features",
    "Raw p < 0.05",
    "FDR discoveries",
    "Bonferroni discoveries",
    "Empirical Bayes discoveries",
    "Top-100 overlap",
    "Effect-size/significance correlation"
  ),
  value = c(
    nrow(sample_metadata),
    sum(sample_metadata$group == "Early"),
    sum(sample_metadata$group == "Late"),
    n_all,
    n_raw,
    n_bh,
    n_bonf,
    n_limma,
    length(common),
    cor_effect_sig
  )
)

write_csv(
  summary_recalculated,
  file.path(OUTPUT_DIR, "multiplicity_summary_recalculated.csv")
)

validation_checks <- supplied_summary %>%
  rename(supplied_value = value) %>%
  left_join(
    summary_recalculated %>% rename(recalculated_value = value),
    by = "metric"
  ) %>%
  mutate(
    status = if_else(
      abs(as.numeric(supplied_value) - as.numeric(recalculated_value)) < 1e-8,
      "PASS",
      "FAIL"
    )
  )

write_csv(
  validation_checks,
  file.path(OUTPUT_DIR, "validation_checks.csv")
)

if (any(validation_checks$status == "FAIL", na.rm = TRUE)) {
  warning("One or more validation checks failed. Review output/validation_checks.csv.")
}

# ---- 4. Sample-structure diagnostics --------------------------------------

split_ids <- strsplit(sample_metadata$sample_id, "-", fixed = TRUE)
patient_id <- vapply(
  split_ids,
  function(x) paste(x[1:3], collapse = "-"),
  character(1)
)
sample_type_code <- vapply(
  split_ids,
  function(x) substr(x[4], 1, 2),
  character(1)
)

sample_structure <- tibble(
  sample_structure = c(
    "Retained expression samples",
    "Unique TCGA participant IDs",
    "Participants with >1 retained sample",
    "Primary tumor sample code 01",
    "Metastatic sample code 06",
    "Solid tissue normal sample code 11"
  ),
  count = c(
    nrow(sample_metadata),
    length(unique(patient_id)),
    sum(table(patient_id) > 1),
    sum(sample_type_code == "01"),
    sum(sample_type_code == "06"),
    sum(sample_type_code == "11")
  )
)

write_csv(
  sample_structure,
  file.path(OUTPUT_DIR, "sample_structure.csv")
)

write_csv(
  tibble(feature_id = sort(common)),
  file.path(OUTPUT_DIR, "top100_overlap_features.csv")
)

# ---- 5. Visual style -------------------------------------------------------

bg      <- "#07111f"
grid    <- "#24384F"
txt     <- "white"
subtxt  <- "#A9BED6"
cyan    <- "#55DDE0"
purple  <- "#7B61FF"
gold    <- "#FFD166"
muted   <- "#334A63"
red     <- "#FF6B6B"

theme_story <- function(base_size = 18) {
  theme_void(base_size = base_size) +
    theme(
      plot.background = element_rect(fill = bg, color = NA),
      panel.background = element_rect(fill = bg, color = NA),
      plot.title = element_text(color = txt, face = "bold", size = base_size + 10),
      plot.subtitle = element_text(color = subtxt, size = base_size + 2),
      plot.margin = margin(35, 45, 35, 45)
    )
}

save_story <- function(p, filename, width = 14, height = 8) {
  ggsave(
    file.path(FIGURE_DIR, filename),
    p,
    width = width,
    height = height,
    dpi = 300,
    bg = bg
  )
}

# ---- 6. Figure 1: multiplicity narrowing ----------------------------------

narrow_df <- data.frame(
  Stage = factor(
    c("Raw findings", "FDR discoveries", "Strict evidence"),
    levels = c("Raw findings", "FDR discoveries", "Strict evidence")
  ),
  Count = c(n_raw, n_bh, n_bonf),
  y = c(3, 2, 1)
) %>%
  mutate(
    width = sqrt(Count / max(Count)) * 10,
    xmin = -width / 2,
    xmax = width / 2,
    ymin = y - 0.28,
    ymax = y + 0.28,
    label = paste0(comma(Count), "\n", Stage)
  )

p_narrow <- ggplot(narrow_df) +
  geom_rect(
    aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax, fill = Stage),
    alpha = 0.96
  ) +
  geom_text(
    aes(x = 0, y = y, label = label),
    color = txt,
    fontface = "bold",
    size = 6,
    lineheight = 0.9
  ) +
  geom_segment(
    aes(x = 0, xend = 0, y = 2.7, yend = 2.35),
    arrow = arrow(length = grid::unit(0.18, "inches")),
    color = subtxt,
    linewidth = 1
  ) +
  geom_segment(
    aes(x = 0, xend = 0, y = 1.7, yend = 1.35),
    arrow = arrow(length = grid::unit(0.18, "inches")),
    color = subtxt,
    linewidth = 1
  ) +
  scale_fill_manual(values = c(cyan, purple, gold)) +
  coord_cartesian(xlim = c(-5.8, 5.8), ylim = c(0.4, 3.6)) +
  labs(
    title = "Multiplicity narrows discovery",
    subtitle = "Thousands of apparent findings collapse into a small evidence set"
  ) +
  theme_story() +
  theme(legend.position = "none")

save_story(p_narrow, "01_multiplicity_narrowing_corridor.png")

# ---- 7. Figure 2: empirical Bayes signal field ----------------------------

set.seed(123)

dots_ind <- data.frame(
  Method = "Independent testing",
  x = runif(n_bh, 0.05, 0.42),
  y = runif(n_bh, 0.25, 0.75)
)

dots_eb <- data.frame(
  Method = "Empirical Bayes",
  x = runif(n_limma, 0.58, 0.95),
  y = runif(n_limma, 0.08, 0.92)
)

dots <- bind_rows(dots_ind, dots_eb)

p_eb_dots <- ggplot(dots, aes(x, y, color = Method)) +
  geom_point(alpha = 0.55, size = 2.1) +
  annotate(
    "text", x = 0.24, y = 1.02,
    label = paste0("Independent testing\n", comma(n_bh), " discoveries"),
    color = txt, fontface = "bold", size = 6
  ) +
  annotate(
    "text", x = 0.76, y = 1.02,
    label = paste0("Empirical Bayes\n", comma(n_limma), " discoveries"),
    color = txt, fontface = "bold", size = 6
  ) +
  annotate(
    "text", x = 0.5, y = -0.06,
    label = paste0(round(n_limma / n_bh, 1), "x more discoveries from the same data"),
    color = gold, fontface = "bold", size = 7
  ) +
  scale_color_manual(values = c("Independent testing" = cyan, "Empirical Bayes" = gold)) +
  coord_cartesian(xlim = c(0, 1), ylim = c(-0.1, 1.1), clip = "off") +
  labs(
    title = "Borrowing strength reveals additional signals",
    subtitle = "Empirical Bayes moderation stabilizes variance estimation across features"
  ) +
  theme_story() +
  theme(legend.position = "none")

save_story(p_eb_dots, "02_empirical_bayes_signal_field.png")

# ---- 8. Figure 3: significance versus effect size -------------------------

limma_plot <- limma_plot %>%
  mutate(
    Highlight = case_when(
      feature_id %in% common ~ "Shared",
      feature_id %in% top100_p$feature_id ~ "Significance only",
      feature_id %in% top100_fc$feature_id ~ "Effect size only",
      TRUE ~ "Other"
    )
  )

p_scatter <- ggplot(limma_plot, aes(absFC, neglogFDR, color = Highlight)) +
  geom_point(
    data = filter(limma_plot, Highlight == "Other"),
    alpha = 0.12,
    size = 0.8
  ) +
  geom_point(
    data = filter(limma_plot, Highlight != "Other"),
    alpha = 0.85,
    size = 2.2
  ) +
  scale_color_manual(
    values = c(
      "Other" = muted,
      "Significance only" = cyan,
      "Effect size only" = purple,
      "Shared" = gold
    )
  ) +
  labs(
    title = "Significance and importance are related, not identical",
    subtitle = paste0(
      "Top 100 by significance and top 100 by effect size shared only ",
      length(common),
      " features; r = ",
      round(cor_effect_sig, 2)
    ),
    x = "Absolute effect size |logFC|",
    y = "-log10(FDR-adjusted p-value)",
    color = NULL
  ) +
  theme_minimal(base_size = 18) +
  theme(
    plot.background = element_rect(fill = bg, color = NA),
    panel.background = element_rect(fill = bg, color = NA),
    panel.grid.major = element_line(color = grid, linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.title = element_text(color = txt, face = "bold", size = 28),
    plot.subtitle = element_text(color = subtxt, size = 19),
    axis.title = element_text(color = subtxt),
    axis.text = element_text(color = subtxt),
    legend.position = "bottom",
    legend.background = element_rect(fill = bg),
    legend.text = element_text(color = subtxt),
    plot.margin = margin(30, 40, 30, 40)
  )

save_story(p_scatter, "03_significance_vs_effect_size_scatter.png")

# ---- 9. Figure 4: top-100 overlap -----------------------------------------

overlap_df <- data.frame(
  Category = factor(
    c("Significance only", "Shared", "Effect size only"),
    levels = c("Significance only", "Shared", "Effect size only")
  ),
  Count = c(
    length(setdiff(top100_p$feature_id, top100_fc$feature_id)),
    length(common),
    length(setdiff(top100_fc$feature_id, top100_p$feature_id))
  )
)

p_overlap <- ggplot(overlap_df, aes(Category, Count, fill = Category)) +
  geom_col(width = 0.58, alpha = 0.96) +
  geom_text(
    aes(label = Count),
    vjust = -0.35,
    color = txt,
    fontface = "bold",
    size = 9
  ) +
  scale_fill_manual(values = c(cyan, gold, purple)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.18))) +
  labs(
    title = "Two rankings tell different stories",
    subtitle = "Only a small core is shared between statistical certainty and magnitude",
    x = NULL,
    y = "Features"
  ) +
  theme_minimal(base_size = 18) +
  theme(
    plot.background = element_rect(fill = bg, color = NA),
    panel.background = element_rect(fill = bg, color = NA),
    panel.grid.major = element_line(color = grid, linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.title = element_text(color = txt, face = "bold", size = 28),
    plot.subtitle = element_text(color = subtxt, size = 19),
    axis.title = element_text(color = subtxt),
    axis.text = element_text(color = subtxt),
    legend.position = "none"
  )

save_story(p_overlap, "04_top100_overlap_panel.png")

# ---- 10. Figure 5: discovery to evidence ----------------------------------

journey <- data.frame(
  Step = factor(
    c("Tested", "Raw", "FDR", "Bonferroni", "Empirical Bayes"),
    levels = c("Tested", "Raw", "FDR", "Bonferroni", "Empirical Bayes")
  ),
  Count = c(n_all, n_raw, n_bh, n_bonf, n_limma)
) %>%
  mutate(
    x = as.numeric(Step),
    y = log10(Count),
    label = comma(Count)
  )

p_journey <- ggplot(journey, aes(x, y)) +
  geom_line(color = cyan, linewidth = 1.5, alpha = 0.55) +
  geom_point(aes(color = Step), size = 9) +
  geom_text(
    aes(label = label),
    vjust = -1.2,
    color = txt,
    fontface = "bold",
    size = 6
  ) +
  scale_color_manual(values = c(muted, cyan, purple, gold, red)) +
  scale_x_continuous(
    breaks = journey$x,
    labels = journey$Step
  ) +
  scale_y_continuous(
    labels = function(x) comma(round(10^x)),
    expand = expansion(mult = c(0.08, 0.22))
  ) +
  labs(
    title = "From discovery to evidence",
    subtitle = "Statistical decisions reshape the evidence landscape",
    x = NULL,
    y = "Features, log scale"
  ) +
  theme_minimal(base_size = 18) +
  theme(
    plot.background = element_rect(fill = bg, color = NA),
    panel.background = element_rect(fill = bg, color = NA),
    panel.grid.major = element_line(color = grid, linewidth = 0.4),
    panel.grid.minor = element_blank(),
    plot.title = element_text(color = txt, face = "bold", size = 28),
    plot.subtitle = element_text(color = subtxt, size = 19),
    axis.title = element_text(color = subtxt),
    axis.text = element_text(color = subtxt),
    legend.position = "none"
  )

save_story(p_journey, "05_from_discovery_to_evidence_journey.png")

# ---- 11. Console summary ---------------------------------------------------

print(summary_recalculated)
print(sample_structure)
cat("\nAnalysis complete. Figures are in figures/ and tables are in output/.\n")