# Heatmap visualization: Effect size and significance
library(ggplot2)
library(dplyr)
library(tidyr)
library(reshape2)

# Set working directory and aesthetic parameters for divergent color scaling
base_dir <- "/Desktop/"
colores <- c("white", "#FFC9B6", "#FF9576", "#E47155", "#FF7431", "#FF5119", "#FF5119")
breaks <- c(0, 1e-10, 0.2, 0.5, 0.8, 1.2, 2, Inf)
labels_cat <- c("0", "0 - 0.2", "0.2 - 0.5", "0.5 - 0.8", "0.8 - 1.2", "1.2 - 2", ">2")

# Load pre-processed effect size and significance data
hpc1 <- read.csv(paste0(base_dir, "hedgespc1spp-positivos.csv"))
hpc2 <- read.csv(paste0(base_dir, "hedgespc2spp-positivos.csv"))
hpc3 <- read.csv(paste0(base_dir, "hedgespc3spp-positivos.csv"))

# Standardize column names for merging
colnames(hpc1)[1:2] <- c("spp_ref", "spp_foc")
colnames(hpc2)[1:2] <- c("spp_ref", "spp_foc")
colnames(hpc3)[1:2] <- c("spp_ref", "spp_foc")

pval1 <- read.csv(paste0(base_dir, "pvalues_pc1.csv"))
pval2 <- read.csv(paste0(base_dir, "pvalues_pc2.csv"))
pval3 <- read.csv(paste0(base_dir, "pvalues_pc3.csv"))

# Define species order and scientific names mapping for plot labeling
orden_deseado <- c("carolinus", "aurifrons", "santacruzi", "radiolatus", "superciliaris", 
                   "uropygialis", "hoffmannii", "rubricapillus", "chrysogenys", "pygmaeus", 
                   "flavifrons", "hypopolius", "cruentatus", "chrysauchen", "pulcher", 
                   "pucherani", "herminieri", "portoricensis", "erythrocephalus", "lewis", 
                   "formicivorus", "cactorum", "candidus", "striatus", "percussus", 
                   "nuchalis", "ruber", "varius", "thyroideus")

nombre_completo <- c(
  "carolinus" = "M. carolinus", "aurifrons" = "M. aurifrons", "santacruzi" = "M. santacruzi",
  "radiolatus" = "M. radiolatus", "superciliaris" = "M. superciliaris", "uropygialis" = "M. uropygialis",
  "hoffmannii" = "M. hoffmannii", "rubricapillus" = "M. rubricapillus", "chrysogenys" = "M. chrysogenys",
  "pygmaeus" = "M. pygmaeus", "flavifrons" = "M. flavifrons", "hypopolius" = "M. hypopolius",
  "cruentatus" = "M. cruentatus", "chrysauchen" = "M. chrysauchen", "pulcher" = "M. pulcher",
  "pucherani" = "M. pucherani", "herminieri" = "M. herminieri", "portoricensis" = "M. portoricensis",
  "erythrocephalus" = "M. erythrocephalus", "lewis" = "M. lewis", "formicivorus" = "M. formicivorus",
  "cactorum" = "M. cactorum", "candidus" = "M. candidus", "striatus" = "M. striatus",
  "percussus" = "X. percussus", "nuchalis" = "S. nuchalis", "ruber" = "S. ruber",
  "varius" = "S. varius", "thyroideus" = "S. thyroideus"
)
orden_deseado_completo <- unname(nombre_completo)

# Define processing function: merges Hedges' g values with p-values, formats triangular matrices, and applies taxonomic labels
procesar_datos_heatmap <- function(hpc_data, pval_data) {
  pval_clean <- pval_data %>%
    mutate(spp_ref = factor(spp_ref, levels = orden_deseado),
           spp_foc = factor(spp_foc, levels = orden_deseado),
           row_num = as.numeric(spp_ref), col_num = as.numeric(spp_foc)) %>%
    filter(row_num <= col_num) %>%
    mutate(significant = p.value < 0.05) %>%
    mutate(spp_ref = dplyr::recode(as.character(spp_ref), !!!nombre_completo),
           spp_foc = dplyr::recode(as.character(spp_foc), !!!nombre_completo)) %>%
    mutate(spp_ref = factor(spp_ref, levels = orden_deseado_completo),
           spp_foc = factor(spp_foc, levels = orden_deseado_completo))
  
  hpc_matrix <- hpc_data %>%
    dplyr::select(spp_ref, spp_foc, hedg_g) %>%
    bind_rows(hpc_data %>% dplyr::select(spp_ref, spp_foc, hedg_g) %>% rename(spp_ref = spp_foc, spp_foc = spp_ref)) %>%
    distinct(spp_ref, spp_foc, .keep_all = TRUE) %>%
    filter(spp_ref != spp_foc) %>%
    mutate(spp_ref = factor(spp_ref, levels = orden_deseado),
           spp_foc = factor(spp_foc, levels = orden_deseado)) %>%
    complete(spp_ref, spp_foc, fill = list(hedg_g = 0)) %>%
    pivot_wider(names_from = spp_foc, values_from = hedg_g, values_fill = 0) %>%
    as.data.frame() %>% {rownames(.) <- .$spp_ref; .} %>%
    dplyr::select(-spp_ref) %>% as.matrix()
  
  hpc_clean <- reshape2::melt(hpc_matrix, varnames = c("spp_ref", "spp_foc"), value.name = "hedg_g") %>%
    mutate(spp_ref = factor(spp_ref, levels = orden_deseado),
           spp_foc = factor(spp_foc, levels = orden_deseado),
           row_num = as.numeric(spp_ref), col_num = as.numeric(spp_foc)) %>%
    filter(row_num <= col_num) %>%
    mutate(hedg_g_cat = cut(hedg_g, breaks = breaks, include.lowest = TRUE)) %>%
    mutate(spp_ref = dplyr::recode(as.character(spp_ref), !!!nombre_completo),
           spp_foc = dplyr::recode(as.character(spp_foc), !!!nombre_completo)) %>%
    mutate(spp_ref = factor(spp_ref, levels = orden_deseado_completo),
           spp_foc = factor(spp_foc, levels = orden_deseado_completo))
  
  return(list(hpc = hpc_clean, pval = pval_clean))
}

# Define visualization function: creates heatmaps with significance markers and italicized scientific names
crear_heatmap <- function(datos_list) {
  ggplot(datos_list$hpc, aes(x = spp_foc, y = spp_ref)) +
    geom_tile(aes(fill = hedg_g_cat), color = "#4a4a4a", linewidth = 0.5) +
    geom_point(data = filter(datos_list$pval, significant),
               aes(x = spp_foc, y = spp_ref), shape = 8, size = 2, color = "black", stroke = 0.8) +
    scale_fill_manual(name = "Hedges' g", values = colores, labels = labels_cat, na.value = "gray90") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 17, face = "italic"),
          axis.text.y = element_text(hjust = 1, size = 13, face = "italic"),
          panel.grid = element_blank(),
          legend.position = "right",
          text = element_text(family = "Times New Roman"),
          axis.ticks.x = element_line(linewidth = 0.9, color = "black")) +
    labs(x = "", y = "") +
    coord_fixed()
}

# Process and export high-resolution JPEG files for PC1, PC2, and PC3
# PC1
data_pc1 <- procesar_datos_heatmap(hpc1, pval1)
plot_pc1 <- crear_heatmap(data_pc1)
ggsave(paste0(base_dir, "plotpc1.jpeg"), plot = plot_pc1, dpi = 800, width = 30, height = 30, units = "cm")

# PC2
data_pc2 <- procesar_datos_heatmap(hpc2, pval2)
plot_pc2 <- crear_heatmap(data_pc2)
ggsave(paste0(base_dir, "plotpc2.jpeg"), plot = plot_pc2, dpi = 800, width = 30, height = 30, units = "cm")

# PC3
data_pc3 <- procesar_datos_heatmap(hpc3, pval3)
plot_pc3 <- crear_heatmap(data_pc3)
ggsave(paste0(base_dir, "plotpc3.jpeg"), plot = plot_pc3, dpi = 800, width = 30, height = 30, units = "cm")


