# ======================================================
# 1. PAQUETES
# ======================================================
library(readxl)
library(dplyr)
library(stringr)
library(ggplot2)
library(tidyr)

# ======================================================
# 2. ARCHIVO
# ======================================================
archivo <- "C:/Users/Lenovo/Downloads/Cuantificación LD Mitofection m-hMito D7.xlsx"

# ======================================================
# 3. SELECCIONAR PESTAÑAS LD
# ======================================================
pestanas <- excel_sheets(archivo)
pestanas_LD <- pestanas[str_detect(pestanas, "^LD")]

# ======================================================
# 4. LECTURA + LIMPIEZA DE DATOS
# ======================================================
datos_LD <- lapply(pestanas_LD, function(p) {
  
  read_excel(archivo, sheet = p, col_types = "text") %>%
    select(ImageNumber, ObjectNumber, AreaShape_Area) %>%
    mutate(
      Condicion = str_remove(str_trim(p), "^LD\\s*"),
      Area_pixels = suppressWarnings(
        as.numeric(str_replace(AreaShape_Area, ",", "."))
      ),
      Area_um2 = Area_pixels * 0.001752
    ) %>%
    filter(!is.na(Area_um2), Area_um2 > 0)
  
}) %>%
  bind_rows()

# ======================================================
# 5. DEFINIR RANGOS (BINS)
# ======================================================
rangos <- c(0, 2, 5, 10, 15, 20, 25, 30, 40)

labels_rangos <- c(
  "0–2", "2–5", "5–10", "10–15",
  "15–20", "20–25", "25–30", ">30"
)

datos_LD <- datos_LD %>%
  filter(Area_um2 < Inf) %>%
  mutate(
    Rango_area = cut(
      Area_um2,
      breaks = rangos,
      labels = labels_rangos,
      include.lowest = TRUE,
      right = FALSE
    )
  ) %>%
  filter(!is.na(Rango_area))

# ======================================================
# 6. FRECUENCIA RELATIVA POR CONDICIÓN
# ======================================================
frecuencia_rangos <- datos_LD %>%
  group_by(Condicion, Rango_area) %>%
  summarise(n_gotas = n(), .groups = "drop") %>%
  group_by(Condicion) %>%
  mutate(freq_rel = n_gotas / sum(n_gotas) * 100) %>%
  ungroup()

# ======================================================
# 7. CHI-CUADRADO Y RESIDUOS
# ======================================================
tabla_dist <- table(datos_LD$Condicion, datos_LD$Rango_area)
chisq_res <- suppressWarnings(chisq.test(tabla_dist))

residuos <- as.data.frame(as.table(chisq_res$stdres))
colnames(residuos) <- c("Condicion", "Rango_area", "residuo")

frecuencia_rangos <- frecuencia_rangos %>%
  left_join(
    residuos %>%
      mutate(
        Condicion = as.character(Condicion),
        Rango_area = as.character(Rango_area)
      ),
    by = c("Condicion", "Rango_area")
  ) %>%
  
  mutate(
    residuo_cat = case_when(
      residuo >  2 ~ "Positive",
      residuo < -2 ~ "Negative",
      TRUE         ~ NA_character_
    )
  )

frecuencia_rangos <- frecuencia_rangos %>%
  mutate(
    sig_label = case_when(
      abs(residuo) >= 4 ~ "***",
      abs(residuo) >= 3 ~ "**",
      abs(residuo) >= 2 ~ "*",
      TRUE              ~ NA_character_
    )
  )

sig_points <- frecuencia_rangos %>%
  filter(!is.na(sig_label), Condicion != "Non-transferred")


# ======================================================
# 8. FACTORES + POSICIÓN DE LOS PUNTOS
# ======================================================
frecuencia_rangos <- frecuencia_rangos %>%
  mutate(
    Condicion = recode(Condicion, "C-" = "Non-transferred"),
    Condicion = factor(Condicion,
                       levels = c("Non-transferred", "mMito", "hMito")
    ),
    Rango_area = factor(Rango_area,
                        levels = labels_rangos,
                        ordered = TRUE
    ),
    residuo_cat = factor(
      residuo_cat,
      levels = c("Positive", "Negative")
    ),
    y_point = freq_rel + 1     # <<< CLAVE
  )

frecuencia_rangos <- frecuencia_rangos %>%
  mutate(
    x_barra = interaction(Rango_area, Condicion, drop = TRUE),
    y_point = freq_rel + 1
  )

sig_points <- frecuencia_rangos %>%
  filter(!is.na(residuo_cat), Condicion != "Non-transferred")

dodge <- position_dodge2(width = 0.9, preserve = "single")

# ======================================================
# 9. GRÁFICO FINAL
# ======================================================
ggplot(frecuencia_rangos,
       aes(x = Condicion, y = freq_rel, fill = Condicion)) +
  
  geom_col(
    color = "black",
    linewidth = 0.2,
    width = 0.8
  ) +
  
  geom_text(
    data = sig_points,
    aes(
      y = y_point,
      label = sig_label,
      color = residuo_cat
    ),
    size = 5,
    fontface = "bold"
  ) +
  
  facet_grid(~ Rango_area, switch = "x") +
  
  
  scale_color_manual(
    name = "χ² residual",
    values = c(
      "Positive" = "black",
      "Negative" = "black"
    ),
    labels = c(
      "Positive" = "Positive residual",
      "Negative" = "Negative residual"
    )
  ) +
  
  scale_fill_manual(
    values = c(
      "Non-transferred" = "black",
      "mMito"           = "#FF8000",
      "hMito"           = "#800000"
      )
  ) +
  
  labs(
    x = "Lipid droplet area (µm²)",
    y = "Relative frequency (%)",
    title = "Lipid droplet size distribution"
  ) +
 
  theme_classic() +
  theme(
    axis.text.x      = element_blank(),
    axis.ticks.x     = element_blank(),
    axis.line.x      = element_line(color = "black", linewidth = 1),
    axis.title.x = element_text(size = 12), # nombre eje X
    axis.title.y = element_text(size = 12), # nombre eje Y
    legend.title     = element_blank()
  )