---
title: "Celigo_Boxplot-streamlined"
output: html_document
---

```{r boxplot-prep}
# === ------------------------ ===
# === Run Once Before Plotting ===
# === ------------------------ ===

# === LIBRARIES ===
library(ggplot2)
library(tidyr)
library(dplyr)
library(data.table)
library(rstudioapi)

# === DATA IMPORT ===
runs_included = 1

read_multiple_files <- function(n_files = runs_included) {
  data_list <- list()

  for(i in 1:n_files) {
    cat("Selecting file", i, "of", n_files, "\n")

    user_input <- showPrompt(title = "Organizing Data:", message = "Enter dataset name. Next select the corresponding tabular .csv and condition .csv")
    
    file_path <- file.choose()
    data <- read.csv(file_path, skip = 15, header = TRUE, check.names = TRUE) %>%
      dplyr::select(1:39) %>%
      filter((Row %in% LETTERS[1:8]) & (Column %in% 1:12))
    
    condition_path <- file.choose()
    Condition <- read.csv(condition_path, check.names = FALSE)
    data <- data %>%
      left_join(Condition %>% select(Well, Condition), by = "Well")
    
    data$Source <- user_input
    data_list[[i]] <- data
    
    cat("Combining", file_path, "and", condition_path, "\n")
    cat("Added dataset:", user_input, "\n\n")
  }

  combined_data <- do.call(rbind, data_list)
  return(combined_data)
}

NET_Data_Combined <- read_multiple_files(runs_included)

NET_Data_Combined <- mutate(NET_Data_Combined, across(starts_with("X.."), ~as.numeric(gsub("%", "", .x)) / 100))
NET_Data_Combined$Condition <- gsub("\\\\n", "\n", NET_Data_Combined$Condition)
names(NET_Data_Combined) <- gsub("^X\\.\\.", "%", names(NET_Data_Combined))

NET_Data_Combined$MPO.Index <- as.numeric(NET_Data_Combined$'%MPO') * as.numeric(NET_Data_Combined$'AVG.MPO.Mean.Intensity')
NET_Data_Combined$CitH3.Index <- as.numeric(NET_Data_Combined$'%CitH3') * as.numeric(NET_Data_Combined$'AVG.CitH3.Mean.Intensity')
NET_Data_Combined$DAPI.Index <- as.numeric(NET_Data_Combined$'%DAPI') * as.numeric(NET_Data_Combined$'AVG.DAPI.Mean.Intensity')

setDT(NET_Data_Combined)

# === HELPER FUNCTIONS ===
draw_bracket <- function(p, x1, x2, y_pos, p_val, bracket_height, text_offset) {
  p +
    annotate("segment", x = x1, xend = x2, y = y_pos, yend = y_pos, linewidth = 0.5) +
    annotate("segment", x = x1, xend = x1, y = y_pos, yend = y_pos - bracket_height, linewidth = 0.5) +
    annotate("segment", x = x2, xend = x2, y = y_pos, yend = y_pos - bracket_height, linewidth = 0.5) +
    annotate("text", x = (x1 + x2) / 2, y = y_pos + text_offset,
             label = paste("p =", round(p_val, 3)), size = 4)
}

# === DATA PREPARATION ===
prepare_plot_data <- function(raw_data, runs, conditions, measures) {
  plot_data <- raw_data %>%
    filter(Condition %in% conditions, Source %in% runs) %>%
    mutate(Condition = factor(Condition, levels = conditions)) %>%
    pivot_longer(cols = all_of(measures), names_to = "Measure", values_to = "value")
  
  return(plot_data)
}

# === STATISTICAL TESTS ===
calculate_wilcoxon <- function(plot_data, condition_pairs) {
  results <- data.frame()
  
  for (pair in condition_pairs) {
    pair_data <- plot_data %>% filter(Condition %in% pair)
    if (nrow(pair_data) > 0) {
      test_result <- wilcox.test(value ~ Condition, data = pair_data)
      results <- rbind(results, 
                      data.frame(cond1 = as.character(pair[1]),
                               cond2 = as.character(pair[2]),
                               p_value = test_result$p.value))
    }
  }
  
  cat("\n=== Wilcoxon P-values ===\n")
  print(results)
  return(results)
}
```

``` {r plot-chunk}
# === MAIN SCRIPT ===

# User settings
runs_to_plot <- c("Run 9")
conditions_to_plot <- c("- PMA", "+ PMA")
measure_to_plot <- c("CitH3.Index")

# Prepare data
plot_data <- prepare_plot_data(NET_Data_Combined, runs_to_plot, conditions_to_plot, measure_to_plot)

# Generate all pairwise comparisons
condition_pairs <- combn(conditions_to_plot, 2, simplify = FALSE)

# Calculate statistics
pairwise_results <- calculate_wilcoxon(plot_data, condition_pairs)

# Format correction
measure_to_plot <- gsub("\\.", " ", measure_to_plot)

# Create plot
p <- ggplot(plot_data, aes(x = Condition, y = value)) +
  geom_boxplot(width = 0.6, alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.2, size = 2.5, alpha = 0.8) +
  labs(x = "", y = measure_to_plot, 
       title = paste0(runs_to_plot, " ", measure_to_plot)) +
  theme_classic() +
  theme(
    plot.title = element_text(size = 16, hjust = 0.5),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text.y = element_text(size = 14),
    axis.title.y = element_text(size = 16)
  )

# Add brackets
max_value <- max(plot_data$value, na.rm = TRUE)
bracket_height <- max_value * 0.05
text_offset <- max_value * 0.03

for (i in seq_along(condition_pairs)) {
  pair <- condition_pairs[[i]]
  p_val <- pairwise_results$p_value[i]
  y_pos <- max_value * 1.15 + (i - 1) * max_value * 0.1
  
  cond1_x <- which(conditions_to_plot == pair[1])
  cond2_x <- which(conditions_to_plot == pair[2])
  
  p <- draw_bracket(p, cond1_x, cond2_x, y_pos, p_val, bracket_height, text_offset)
}

# Display and save
print(p)
ggsave("boxplot.png", plot = p, width = 8, height = 6)
cat("\nPlot saved to: boxplot.png\n")
```