library(magrittr)
library(VennDiagram)
library(gridExtra)
library(grid)
library(plotrix)
library(limma)
library(dendextend)
library(edgeR)
library(eulerr)
library(stringr)
library(STRINGdb)
library(janitor)
library(rlang)
library(cowplot)
library(svglite)
library(ggalluvial)
library(parallel)
library(tidyverse)
library(ggrepel)
library(ggpubr)
library(emmeans)       
library(ggdendro)
library(AnnotationHub)
library(pheatmap)
library(RColorBrewer)
library(viridis)

# --- SETUP ---------------------------------------------------
Expt_Summary_All <- read.csv("Visits_df_22Oct25.csv")

Expt_Summary_All$Treat_Geno <- paste0(Expt_Summary_All$Genotype, '\n', Expt_Summary_All$Treatment_Index)
Expt_Summary_All$Treat_Geno <- factor(Expt_Summary_All$Treat_Geno,
                                       levels = c('WT\nVehicle', 'WT\nESI05', 'KO\nVehicle', 'KO\nESI05'))

# --- SHARED COLOR SCALE (defined once, reused throughout) ----
group_colors <- c('KO\nESI05'   = "#a31d47",
                  'KO\nVehicle' = "#ED2224",
                  'WT\nVehicle' = "grey0",
                  'WT\nESI05'   = "grey38")

sex_shapes <- c('F' = 1, 'M' = 0)

# ============================================================
# FUNCTION: SC_WHISKER_PLOTS_FIGURES
# ============================================================
SC_WHISKER_PLOTS_FIGURES <- function(DF, Expt_Type, Col2plt, ylabel, xcategory, ggtitle, Filename, ylimit) {
  p <- ggplot(DF %>% filter(Expt_Type == {{Expt_Type}}),
              aes(x = {{xcategory}}, y = {{Col2plt}}, color = {{xcategory}}, fill = {{xcategory}})) +
    stat_summary(fun = mean, geom = "point", size = 10, shape = 95) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, size = 0.5) +
    geom_jitter(aes(shape = factor(Sex, levels = c('M', 'F'))), fill = NA,
                position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.75),
                size = 2, alpha = 0.7) +
    theme_cowplot() +
    ylab(ylabel) +
    geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
    scale_color_manual(values = group_colors) +
    scale_fill_manual(values = group_colors) +
    scale_shape_manual(values = sex_shapes) +
    theme(legend.title = element_blank()) +
    theme(axis.title.x = element_blank(), legend.position = "none",
          plot.title = element_blank(), axis.title.y = element_blank(),
          axis.text.x = element_blank()) +
    ggtitle(ggtitle) +
    ylim(ylimit)

  print(p)
  ggsave(paste0(Filename, Sys.Date(), ".pdf"), dpi = 1000, width = 1.75, height = 1.5, units = "in")
}

# ============================================================
# FUNCTION: SC_WHISKER_PROPORTION_FIGURES
# ============================================================
SC_WHISKER_PROPORTION_FIGURES <- function(DF, Expt_Type, Col2plt, ylabel, xcategory, ggtitle, Filename, ylimit) {
  p <- ggplot(DF %>% filter(Expt_Type == {{Expt_Type}}),
              aes(x = {{xcategory}}, y = {{Col2plt}}, color = {{xcategory}}, fill = {{xcategory}})) +
    stat_summary(fun = mean, geom = "point", size = 10, shape = 95) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, size = 0.5) +
    geom_jitter(aes(shape = factor(Sex, levels = c('M', 'F'))), fill = NA,
                position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.75),
                size = 2, alpha = 0.7) +
    theme_cowplot() +
    ylab(ylabel) +
    ylim(ylimit) +
    geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
    scale_color_manual(values = group_colors) +
    scale_fill_manual(values = group_colors) +
    scale_shape_manual(values = sex_shapes) +
    theme(legend.title = element_blank()) +
    theme(axis.title.x = element_blank(), legend.position = "none",
          plot.title = element_blank(), axis.title.y = element_blank(),
          axis.text.x = element_blank()) +
    ggtitle(ggtitle)

  print(p)
  ggsave(paste0(Filename, Sys.Date(), ".pdf"),  dpi = 1000, width = 1.75, height = 1.5, units = "in")
  ggsave(paste0(Filename, Sys.Date(), ".jpeg"), dpi = 1000, width = 1.75, height = 1.5, units = "in")
}

# ============================================================
# FUNCTION: SC_WHISKER_PROPORTION_Sidebyside_FIGURES
# ============================================================
SC_WHISKER_PROPORTION_Sidebyside_FIGURES <- function(DF, Expt_Type, Col2plt, ylabel, xcategory, ggtitle, Filename, ylimit) {
  df_long <- DF %>% pivot_longer(cols = Col2plt, values_to = 'measure', names_to = 'colsID')

  p <- ggplot(df_long %>% filter(Expt_Type == {{Expt_Type}}),
              aes(x = interaction(colsID, Treat_Geno), y = measure, color = Treat_Geno, fill = Treat_Geno)) +
    stat_summary(fun = mean, geom = "point", size = 10, shape = 95) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.1, size = 0.5) +
    geom_jitter(aes(shape = factor(Sex, levels = c('M', 'F'))), fill = NA,
                position = position_jitterdodge(jitter.width = 0.5, dodge.width = 0.75),
                size = 2, alpha = 0.7) +
    theme_cowplot() +
    ylab(ylabel) +
    ylim(ylimit) +
    geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
    scale_color_manual(values = group_colors) +
    scale_fill_manual(values = group_colors) +
    scale_shape_manual(values = sex_shapes) +
    theme(legend.title = element_blank()) +
    theme(axis.title.x = element_blank(), legend.position = "none",
          plot.title = element_blank(), axis.title.y = element_blank(),
          axis.text.x = element_blank()) +
    ggtitle(ggtitle)

  print(p)
  ggsave(paste0(Filename, Sys.Date(), ".pdf"),  dpi = 1000, width = 2.75, height = 1.5, units = "in")
  ggsave(paste0(Filename, Sys.Date(), ".jpeg"), dpi = 1000, width = 2.75, height = 1.5, units = "in")
}

# ============================================================
# FUNCTION: Make_Dot_SEM_plot
# FIX: replaced deprecated aes_string() with .data[[]] pronoun
# FIX: added default value for ylimt to prevent errors on calls that omit it
# ============================================================
Make_Dot_SEM_plot <- function(df, what2filter, x_value, y_value, filename, ylimt = NULL) {
  plot <- ggplot(df %>% filter(Behavior == what2filter),
                 aes(x = .data[[x_value]], y = .data[[y_value]],
                     color = .data[[x_value]])) +
    geom_point(aes(shape = Sex), size = 2, alpha = 0.5, fill = 'white', stroke = 0.6,
               position = position_jitter(width = 0.05, height = 0)) +
    stat_summary(fun = mean, geom = "point", size = 10, shape = 95) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.25, size = 0.5) +
    theme_cowplot() +
    scale_color_manual(values = c('KO\nESI05'   = "#a31d47",
                                  'KO\nSaline'  = "#ED2224",
                                  'WT\nSaline'  = "grey0",
                                  'WT\nESI05'   = "grey38")) +
    scale_shape_manual(values = c(21, 22)) +
    theme(legend.position = "none",
          axis.title.x = element_blank(),
          axis.title.y = element_blank())

  if (!is.null(ylimt)) plot <- plot + ylim(ylimt)

  ggsave(paste0(filename, '.svg'), plot = plot, dpi = 900, width = 1.75, height = 2)
}

# ============================================================
# SECTION: PREFERENCE INDEX — CUP (Sociability & Familiarity)
# ============================================================

### Sociability
SC_WHISKER_PLOTS_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                         Col2plt = preference_Obj_time,
                         ylabel = 'Preference Index \n (time interacting with cup)',
                         xcategory = Treat_Geno,
                         ggtitle = 'Preference Index (with Novel mouse)',
                         Filename = '3SC_Plot_Preference_Index_Cup_Sociability_ILLUS',
                         ylimit = c(-0.75, 0.75))

Expt_Summary_All %>%
  filter(Expt_Type == 'Sociability') %>%
  group_by(Treat_Geno) %>%
  summarize(mean_val = mean(preference_Obj_time),
            sd_val   = sd(preference_Obj_time, na.rm = TRUE),
            n        = n(), .groups = 'drop')

anova_results <- aov(preference_Obj_time ~ Treatment_Index * Genotype,
                     data = Expt_Summary_All %>% filter(Expt_Type == 'Sociability'))
summary(anova_results)
TukeyHSD(anova_results)
emmeans(anova_results, pairwise ~ Treatment_Index * Genotype)

### Familiarity
SC_WHISKER_PLOTS_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Familiarity',
                         Col2plt = preference_Obj_time,
                         ylabel = 'Preference Index \n (time interacting with cup)',
                         xcategory = Treat_Geno,
                         ggtitle = 'Preference Index Cup (Habituation)',
                         Filename = '3SC_Plot_Preference_Index_FAMILIARITY_ILLUS',
                         ylimit = c(-0.75, 0.75))

Expt_Summary_All %>%
  filter(Expt_Type == 'Familiarity') %>%
  group_by(Treat_Geno) %>%
  summarize(mean_val = mean(preference_Obj_time),
            sd_val   = sd(preference_Obj_time, na.rm = TRUE),
            n        = n(), .groups = 'drop')

anova_results <- aov(preference_Obj_time ~ Treat_Geno * Sex,
                     data = Expt_Summary_All %>% filter(Expt_Type == 'Familiarity'))
summary(anova_results)
TukeyHSD(anova_results)

# ============================================================
# SECTION: PREFERENCE INDEX CHAMBER (Sociability & Familiarity)
# ============================================================

### Sociability
SC_WHISKER_PLOTS_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                         Col2plt = preference_Cham_time,
                         ylabel = 'Preference Index \n (time interacting in chamber)',
                         xcategory = Treat_Geno,
                         ggtitle = 'Preference Index',
                         Filename = '3SC_Plot_Preference_Index_CHAM_Sociability_ILLUS',
                         ylimit = c(-0.75, 0.75))

Expt_Summary_All %>%
  filter(Expt_Type == 'Sociability') %>%
  group_by(Treat_Geno) %>%
  summarize(number   = n(),
            mean     = mean(preference_Cham_time),
            sd_val   = sd(preference_Cham_time, na.rm = TRUE),
            SEM      = sd(preference_Cham_time) / sqrt(n()))

anova_results <- aov(preference_Cham_time ~ Treat_Geno * Sex,
                     data = Expt_Summary_All %>% filter(Expt_Type == 'Sociability'))
summary(anova_results)
TukeyHSD(anova_results)

### Familiarity
SC_WHISKER_PLOTS_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Familiarity',
                         Col2plt = preference_Cham_time,
                         ylabel = 'Preference Index \n (time interacting in chamber)',
                         xcategory = Treat_Geno,
                         ggtitle = 'Preference Index',
                         Filename = '3SC_Plot_Preference_Index_TIME_CHAM_FAMILIARITY_ILLUS',
                         ylimit = c(-0.75, 0.75))

# ============================================================
# SECTION: PREFERENCE INDEX VISITS TO CUPS (Sociability & Familiarity)
# ============================================================

### Sociability
SC_WHISKER_PLOTS_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                         Col2plt = preference_Obj_visit,
                         ylabel = 'Preference Index \n (unique approaches to the cup)',
                         xcategory = Treat_Geno,
                         ggtitle = 'Preference Index visits (Sociability)',
                         Filename = '3SC_Plot_Preference_Index_Visits_Cup_Sociability_ILLUS',
                         ylimit = c(-0.75, 0.75))

Expt_Summary_All %>%
  filter(Expt_Type == 'Sociability') %>%
  group_by(Treat_Geno) %>%
  summarize(mean_val = mean(preference_Obj_visit),
            sd_val   = sd(preference_Obj_visit, na.rm = TRUE),   # FIX: was sd(preference_Cham_time)
            n        = n(), .groups = 'drop')

anova_results <- aov(preference_Obj_visit ~ Treatment_Index * Genotype,  # FIX: was preference_Cham_time
                     data = Expt_Summary_All %>% filter(Expt_Type == 'Sociability'))
summary(anova_results)
TukeyHSD(anova_results)

### Familiarity
SC_WHISKER_PLOTS_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Familiarity',
                         Col2plt = preference_Obj_visit,
                         ylabel = 'Preference Index \n (unique approaches to the cup)',
                         xcategory = Treat_Geno,
                         ggtitle = 'Preference Index visits (Familiarity)',
                         Filename = '3SC_Plot_Preference_Index_Visits_Cup_Familiarity_FIGURE',
                         ylimit = c(-0.75, 0.75))

# ============================================================
# SECTION: CORRELATION Sociability vs Familiarity
# ============================================================

groups <- c('WT\nESI05', 'WT\nVehicle', 'KO\nESI05', 'KO\nVehicle', 'All')

plot_correlation_plot <- function(DF, groups, what2plot, filename) {
  for (i in groups) {
    if (i == 'All') {
      data <- DF %>% filter(Expt_Type %in% c('Sociability', 'Familiarity'))
    } else {
      data <- DF %>% filter(Expt_Type %in% c('Sociability', 'Familiarity'), Treat_Geno == i)
    }

    p <- data %>%
      pivot_wider(names_from = Expt_Type, values_from = !!sym(what2plot),
                  id_cols = c(MouseID, Treat_Geno, Sex)) %>%
      ggplot(aes(x = Sociability, y = Familiarity)) +
      geom_point(aes(color = Treat_Geno, shape = Sex), size = 5, alpha = 0.7) +
      geom_abline(intercept = 0, slope = 1, linetype = 2, color = 'grey') +
      geom_hline(yintercept = 0, linetype = 2, color = 'grey') +
      geom_vline(xintercept = 0, linetype = 2, color = 'grey') +
      geom_smooth(method = 'lm', formula = y ~ x + 0, se = TRUE, color = 'blue', fill = 'grey') +
      theme_cowplot() +
      xlab(paste('Preference Index (Sociability) -', what2plot)) +
      ylab(paste('Preference Index (Familiarity) -', what2plot)) +
      scale_color_manual(values = group_colors) +
      scale_shape_manual(values = sex_shapes) +
      theme(legend.position = "right") +
      xlim(-1, 1) +   # FIX: was xlim(1,-1) — reversed axes
      ylim(-1, 1)     # FIX: was ylim(1,-1)

    print(p)
    filename_use <- paste0(filename, "_", gsub('\n', '_', i), "_", Sys.Date(), ".pdf")
    ggsave(filename_use, plot = p, dpi = 1000, width = 5, height = 5, units = "in")
  }
}

plot_correlation_plot(Expt_Summary_All, groups, 'preference_Obj_time',  'Correlation_Obj_Time')
plot_correlation_plot(Expt_Summary_All, groups, 'preference_Cham_time', 'Correlation_Cham_Time')

# ============================================================
# SECTION: PROPORTION TIME Side-by-side (Cup & Chamber)
# ============================================================

### Cups: FC vs NC
SC_WHISKER_PROPORTION_Sidebyside_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                                          Col2plt = c("prop_frames_FC", "prop_frames_NC"),
                                          ylabel = 'Proportion of total time',
                                          xcategory = Treat_Geno,
                                          ggtitle = 'Proportion (Novel Object)',
                                          Filename = '3SC_Plot_Proportion_Object_FCvsNC_Sociability_Figure',
                                          ylimit = c(0, 0.4))
SC_WHISKER_PROPORTION_Sidebyside_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Familiarity',
                                          Col2plt = c("prop_frames_FC", "prop_frames_NC"),
                                          ylabel = 'Proportion of total time',
                                          xcategory = Treat_Geno,
                                          ggtitle = 'Proportion (Novel Object)',
                                          Filename = '3SC_Plot_Proportion_Object_FCvsNC_Familiarity_Figure',
                                          ylimit = c(0, 0.4))

df_long_prob_cup <- Expt_Summary_All %>%
  pivot_longer(cols = c("prop_frames_FC", "prop_frames_NC"), values_to = 'measure', names_to = 'colsID')
anova_results <- aov(measure ~ Treat_Geno * colsID,
                     data = df_long_prob_cup %>% filter(Expt_Type == 'Sociability'))
summary(anova_results)
em <- emmeans(anova_results, pairwise ~ Treat_Geno * colsID)
contrast(em, method = "revpairwise", by = "Treat_Geno")

### Chambers: FO vs NO
SC_WHISKER_PROPORTION_Sidebyside_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                                          Col2plt = c("prop_frames_FO", "prop_frames_NO"),
                                          ylabel = 'Proportion of total time',
                                          xcategory = Treat_Geno,
                                          ggtitle = 'Proportion (Novel Chamber)',
                                          Filename = '3SC_Plot_Proportion_Object_FOvsNO_Sociability_Figure',
                                          ylimit = c(0, 1))
SC_WHISKER_PROPORTION_Sidebyside_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Familiarity',
                                          Col2plt = c("prop_frames_FO", "prop_frames_NO"),
                                          ylabel = 'Proportion of total time',
                                          xcategory = Treat_Geno,
                                          ggtitle = 'Proportion (Novel Chamber)',
                                          Filename = '3SC_Plot_Proportion_Object_FOvsNO_Familiarity_Figure',
                                          ylimit = c(0, 1))

df_long_prob_chamber <- Expt_Summary_All %>%
  pivot_longer(cols = c("prop_frames_FO", "prop_frames_NO"), values_to = 'measure', names_to = 'colsID')
anova_results <- aov(measure ~ Treat_Geno * colsID,
                     data = df_long_prob_chamber %>% filter(Expt_Type == 'Sociability'))
summary(anova_results)
em <- emmeans(anova_results, pairwise ~ Treat_Geno * colsID)
contrast(em, method = "revpairwise", by = "Treat_Geno")

anova_results <- aov(measure ~ Treat_Geno * colsID,
                     data = df_long_prob_chamber %>% filter(Expt_Type == 'Familiarity'))
summary(anova_results)
emmeans(anova_results, pairwise ~ Treat_Geno * colsID)

### Cup visits: FC vs NC
SC_WHISKER_PROPORTION_Sidebyside_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                                          Col2plt = c("n_visits_FC", "n_visits_NC"),
                                          ylabel = 'Number of visits',
                                          xcategory = Treat_Geno,
                                          ggtitle = 'Visits (Novel Cup)',
                                          Filename = '3SC_Plot_VISITS_CUP_Sociability_FCvsNC_Figure',
                                          ylimit = c(0, 200))
SC_WHISKER_PROPORTION_Sidebyside_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Familiarity',
                                          Col2plt = c("n_visits_FO", "n_visits_NO"),
                                          ylabel = 'Number of visits',
                                          xcategory = Treat_Geno,
                                          ggtitle = 'Visits (Novel Cup)',
                                          Filename = '3SC_Plot_VISITS_CUP_Familiarity_FCvsNC_Figure',
                                          ylimit = c(0, 100))
SC_WHISKER_PROPORTION_Sidebyside_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                                          Col2plt = c("n_visits_FO", "n_visits_NO"),
                                          ylabel = 'Number of visits',
                                          xcategory = Treat_Geno,
                                          ggtitle = 'Visits (Novel Chamber)',
                                          Filename = '3SC_Plot_VISITS_CHAMBER_Sociability_FOvsNO_Figure',
                                          ylimit = c(0, 100))

df_long_visit_cup <- Expt_Summary_All %>%
  pivot_longer(cols = c("n_visits_FC", "n_visits_NC"), values_to = 'measure', names_to = 'colsID')
anova_results <- aov(measure ~ Treat_Geno * colsID,
                     data = df_long_visit_cup %>% filter(Expt_Type == 'Sociability'))
summary(anova_results)
em <- emmeans(anova_results, ~ Treat_Geno * colsID)
contrast(em, method = "revpairwise", by = "Treat_Geno")

# ============================================================
# SECTION: PROPORTION TIME Individual columns
# ============================================================

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = prop_frames_NC,
                               ylabel = 'Proportion of total time',
                               xcategory = Treat_Geno,
                               ggtitle = 'Proportion (Novel Mouse)',
                               Filename = '3SC_Plot_Proportion_Novel_Object_Sociability_Figure',
                               ylimit = c(0, 0.25))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = prop_frames_FO,
                               ylabel = 'Proportion of total time',
                               xcategory = Treat_Geno,
                               ggtitle = 'Proportion time in chamber \n (Novel Object)',
                               Filename = '3SC_Plot_Proportion_Familiar_Chamber_Sociability_Figure',
                               ylimit = c(0, 1))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = prop_frames_NO,
                               ylabel = 'Proportion of total time',
                               xcategory = Treat_Geno,
                               ggtitle = 'Proportion time chamber \n (Novel Chamber)',
                               Filename = '3SC_Plot_Proportion_Novel_Chamber_Sociability_Figure',
                               ylimit = c(0, 1))

# ============================================================
# SECTION: MEAN TIME PER VISIT
# ============================================================

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = mean_time_FC_sec, ylabel = 'Mean time per visit',
                               xcategory = Treat_Geno,
                               ggtitle = 'Mean time interacting with chamber \n (Novel Object)',
                               Filename = '3SC_Plot_Mean_Novel_Cup_Object_Sociability_Figure',
                               ylimit = c(0, 1.5))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = mean_time_NC_sec, ylabel = 'Mean time per visit',
                               xcategory = Treat_Geno,
                               ggtitle = 'Mean time interacting with chamber \n (Novel Mouse)',
                               Filename = '3SC_Plot_Mean_Novel_Cup_Mouse_Sociability_Figure',
                               ylimit = c(0, 1.5))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = mean_time_NO_sec, ylabel = 'Mean time per visit',
                               xcategory = Treat_Geno,
                               ggtitle = 'Mean time spent in novel mouse chamber',
                               Filename = '3SC_Plot_Mean_Novel_Chamber_Mouse_Sociability_Figure',
                               ylimit = c(0, 50))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = mean_time_FO_sec, ylabel = 'Mean time per visit',
                               xcategory = Treat_Geno,
                               ggtitle = 'Mean time spent in novel object chamber',
                               Filename = '3SC_Plot_Mean_Novel_Chamber_Object_Sociability_Figure',
                               ylimit = c(0, 50))

# ============================================================
# SECTION: VISIT COUNTS
# ============================================================

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = n_visits_FC, ylabel = 'Number of visits',
                               xcategory = Treat_Geno,
                               ggtitle = 'Number of visits interacting with \n Novel Object',
                               Filename = '3SC_Plot_NumberVisits_Novel_Object_Sociability_ILLUS',
                               ylimit = c(0, 200))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = n_visits_NC, ylabel = 'Number of visits',
                               xcategory = Treat_Geno,
                               ggtitle = 'Number of visits interacting with \n Novel Mouse',
                               Filename = '3SC_Plot_NumberVisits_Novel_Mouse_Sociability_Figure',
                               ylimit = c(0, 250))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = n_visits_NO, ylabel = 'Number of visits',
                               xcategory = Treat_Geno,
                               ggtitle = 'Number of visits interacting with Novel mouse chamber',
                               Filename = '3SC_Plot_NumberVisits_Novel_Mouse_Chamber_Sociability_Figure',
                               ylimit = c(0, 100))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = n_visits_FO, ylabel = 'Number of visits',
                               xcategory = Treat_Geno,
                               ggtitle = 'Number of visits interacting with Novel Object chamber',
                               Filename = '3SC_Plot_NumberVisits_Novel_Object_Chamber_Sociability_Figure',
                               ylimit = c(0, 100))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Familiarity',
                               Col2plt = n_visits_FC, ylabel = 'Number of visits',
                               xcategory = Treat_Geno,
                               ggtitle = 'Number of visits interacting with \n Novel Object',
                               Filename = '3SC_Plot_NumberVisits_Novel_Object_Familiarity_Figure',
                               ylimit = c(0, 200))

# ============================================================
# SECTION: DISTANCE
# ============================================================

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Sociability',
                               Col2plt = dist_sum_px, ylabel = 'Distance travelled',
                               xcategory = Treat_Geno,
                               ggtitle = 'Distance traveled Sociability',
                               Filename = '3SC_Plot_DISTANCE_Novel_Object_Sociability_ILLUS_',
                               ylimit = c(25000, 125000))

SC_WHISKER_PROPORTION_FIGURES(DF = Expt_Summary_All, Expt_Type = 'Familiarity',
                               Col2plt = dist_sum_px, ylabel = 'Distance travelled',
                               xcategory = Treat_Geno,
                               ggtitle = 'Distance traveled Habituation',
                               Filename = '3SC_Plot_DISTANCE_Novel_Object_Familiarity_ILLUS_',
                               ylimit = c(0, 150000))

anova_results <- aov(dist_sum_px ~ Treatment_Index * Genotype,
                     data = Expt_Summary_All %>% filter(Expt_Type == 'Sociability'))
summary(anova_results)
TukeyHSD(anova_results)
Expt_Summary_All %>%
  filter(Expt_Type == 'Sociability') %>%
  group_by(Treat_Geno) %>%
  summarize(mean_val = mean(dist_sum_px), sd_val = sd(dist_sum_px, na.rm = TRUE), n = n(), .groups = 'drop')

# ============================================================
# SECTION: PCA
# ============================================================

df_pca <- Expt_Summary_All %>%
  dplyr::select(!c(FileName, MouseID, Treatment, Treatment_Index, Treat_Geno, Batch,
                   Batch_Index, Expt_Type, ExptGroup, Cage.ID, Cage, Genotype, DOB, Sex, Index))

pca_result <- prcomp(df_pca, scale = FALSE, center = TRUE)
plot(pca_result, type = "l")
summary(pca_result)
biplot(pca_result)

pc_scores <- as.data.frame(pca_result$x)
pc_scores$Treatment_Geno <- Expt_Summary_All$Treat_Geno

ggplot(pc_scores, aes(x = PC1, y = PC2, color = Treatment_Geno)) +
  geom_point(size = 3) +
  theme_cowplot()

# ============================================================
# SECTION: HIERARCHICAL CLUSTERING
# ============================================================
# FIX: removed duplicate scale() call (df_pca_scaled); using data_scaled consistently
# FIX: df$Cluster replaced with Expt_Summary_All$Cluster

data_scaled <- scale(df_pca)
hc <- hclust(dist(data_scaled), method = "ward.D2")
plot(hc, main = "Hierarchical Clustering Dendrogram", xlab = "Samples", ylab = "Distance")

clusters <- cutree(hc, k = 5)
Expt_Summary_All$Cluster <- clusters   # FIX: was df$Cluster (df undefined)

hc_dend <- as.dendrogram(hc)
geno_treat_colors <- Expt_Summary_All$Treat_Geno[order.dendrogram(hc_dend)]
color_map <- c('WT\nVehicle' = "grey0", 'WT\nESI05' = "grey38",
               'KO\nVehicle' = "#ED2224", 'KO\nESI05' = "#a31d47")
labels_colors(hc_dend) <- color_map[geno_treat_colors]

plot(hc_dend, main = "Hierarchical Clustering by Genotype/Treatment")
legend("topright", legend = names(color_map), col = color_map, pch = 15, cex = 1)

# ============================================================
# SECTION: SUMMARY BARPLOTS — Group counts
# ============================================================

ggplot(Expt_Summary_All %>%
         filter(Expt_Type == 'Sociability') %>%
         mutate(interaction_var = fct_relevel(interaction(Genotype, Treatment_Index),
                                              "WT.Vehicle", "WT.ESI05", "KO.Vehicle", "KO.ESI05")),
       aes(x = interaction_var, fill = Treat_Geno)) +
  geom_bar(alpha = 0.75) +
  theme_cowplot() +
  xlab("") + ylab("Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  geom_text(stat = 'count', aes(label = after_stat(count)), vjust = -0.5) +
  scale_fill_manual(values = group_colors)

ggsave(paste0('3SC_Sociability_Counts_Groups_Barplot', Sys.Date(), ".pdf"),
       dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_Sociability_Counts_Groups_Barplot', Sys.Date(), ".jpeg"),
       dpi = 1000, width = 6, height = 6, units = "in")

ggplot(Expt_Summary_All %>%
         filter(Expt_Type == 'Sociability') %>%
         mutate(interaction_var = fct_relevel(interaction(Genotype, Sex, Treatment_Index),
                                              "WT.M.Vehicle", "WT.F.Vehicle", "WT.M.ESI05", "WT.F.ESI05",
                                              "KO.M.Vehicle", "KO.F.Vehicle", "KO.M.ESI05", "KO.F.ESI05")),
       aes(x = interaction_var, fill = Treat_Geno)) +
  geom_bar(alpha = 0.75) +
  theme_cowplot() +
  xlab("") + ylab("Count") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  geom_text(stat = 'count', aes(label = after_stat(count)), vjust = -0.5) +
  scale_fill_manual(values = group_colors)

ggsave(paste0('3SC_Sociability_Counts_SEX_Groups_Barplot', Sys.Date(), ".pdf"),
       dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_Sociability_Counts_SEX_Groups_Barplot', Sys.Date(), ".jpeg"),
       dpi = 1000, width = 6, height = 6, units = "in")

# ============================================================
# SECTION: PREFERENCE INDEX BOXPLOTS (Sociability)
# ============================================================



ggplot(Expt_Summary_All %>% filter(Expt_Type == 'Sociability',
                                    !(MouseID %in% c('MO1')),
                                    Treat_Geno != 'ESI050'),
       aes(x = Treat_Geno, y = preference_Obj_time, color = Treat_Geno, fill = Treat_Geno)) +
  geom_boxplot(position = position_dodge(width = 0.75), alpha = 0.25) +
  geom_jitter(aes(shape = Sex), position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.75),
              size = 5, alpha = 0.7) +
  theme_cowplot() +
  ylab('Preference Index \n (time interacting with cup)') + xlab("") +
  scale_color_manual(values = group_colors) +
  geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
  scale_fill_manual(values = group_colors) +
  scale_shape_manual(values = sex_shapes) +
  theme(axis.title.x = element_blank())

ggsave(paste0('3SC_Preference_Index_Boxplot', Sys.Date(), ".pdf"),  dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_Preference_Index_Boxplot', Sys.Date(), ".jpeg"), dpi = 1000, width = 6, height = 6, units = "in")

### Breakdown by Sex
ggplot(Expt_Summary_All %>% filter(Expt_Type == 'Sociability',
                                    !(MouseID %in% c('MO1')),
                                    Treat_Geno != 'ESI050'),
       aes(x = Treat_Geno, y = preference_Obj_time, color = Treat_Geno, fill = Treat_Geno, shape = Sex)) +
  geom_boxplot(position = position_dodge(width = 0.75), alpha = 0.25) +
  geom_jitter(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.75), size = 5, alpha = 0.7) +
  theme_cowplot() +
  ylab('Preference Index \n (time interacting with cup)') +
  scale_color_manual(values = group_colors) +
  geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
  scale_fill_manual(values = group_colors) +
  scale_shape_manual(values = sex_shapes) +
  theme(axis.title.x = element_blank())

ggsave(paste0('3SC_Preference_Index_PerSEX_Boxplot', Sys.Date(), ".pdf"),  dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_Preference_Index_PerSEX_Boxplot', Sys.Date(), ".jpeg"), dpi = 1000, width = 6, height = 6, units = "in")

# ============================================================
# SECTION: PREFERENCE INDEX BOXPLOTS (Familiarity)
# ============================================================

ggplot(Expt_Summary_All %>% filter(Expt_Type == 'Familiarity',
                                    !(MouseID %in% c('MO1')),
                                    Treat_Geno != 'ESI050'),
       aes(x = Treat_Geno, y = preference_Obj_time, color = Treat_Geno, fill = Treat_Geno)) +
  geom_boxplot(position = position_dodge(width = 0.75), alpha = 0.25) +
  geom_jitter(aes(shape = Sex), position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.75),
              size = 5, alpha = 0.7) +
  theme_cowplot() +
  ylab('Preference Index \n (time interacting with cup)') + xlab("") +
  scale_color_manual(values = group_colors) +
  geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
  scale_fill_manual(values = group_colors) +
  scale_shape_manual(values = sex_shapes) +
  theme(axis.title.x = element_blank())

ggsave(paste0('3SC_Preference_FAMILIARITY_Index_Boxplot', Sys.Date(), ".pdf"),  dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_Preference_FAMILIARITY_Index_Boxplot', Sys.Date(), ".jpeg"), dpi = 1000, width = 6, height = 6, units = "in")

### Breakdown by Sex
ggplot(Expt_Summary_All %>% filter(Expt_Type == 'Familiarity',
                                    !(MouseID %in% c('MO1')),
                                    Treat_Geno != 'ESI050'),
       aes(x = Treat_Geno, y = preference_Obj_time, color = Treat_Geno, fill = Treat_Geno, shape = Sex)) +
  geom_boxplot(position = position_dodge(width = 0.75), alpha = 0.25) +
  geom_jitter(position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.75), size = 5, alpha = 0.7) +
  theme_cowplot() +
  ylab('Preference Index \n (time interacting with cup)') +
  scale_color_manual(values = group_colors) +
  geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
  scale_fill_manual(values = group_colors) +
  scale_shape_manual(values = sex_shapes) +
  theme(axis.title.x = element_blank())

ggsave(paste0('3SC_Preference_Index_FAMILIARITY_PerSEX_Boxplot', Sys.Date(), ".pdf"),  dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_Preference_Index_FAMILIARITY_PerSEX_Boxplot', Sys.Date(), ".jpeg"), dpi = 1000, width = 6, height = 6, units = "in")

# ============================================================
# SECTION: TIME INTERACTING (FO and NO) — Sociability
# ============================================================

ggplot(Expt_Summary_All %>% filter(Expt_Type == 'Sociability',
                                    !(MouseID %in% c('MO1')),
                                    Treat_Geno != 'ESI050'),
       aes(x = Treat_Geno, y = time_FO_sec, color = Treat_Geno, fill = Treat_Geno)) +
  geom_boxplot(position = position_dodge(width = 0.75), alpha = 0.25) +
  geom_jitter(aes(shape = Sex), position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.75),
              size = 5, alpha = 0.7) +
  theme_cowplot() +
  ylab('Total time interacting with novel object (seconds)') + xlab("") +
  scale_color_manual(values = group_colors) +
  geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
  scale_fill_manual(values = group_colors) +
  scale_shape_manual(values = sex_shapes) +
  theme(axis.title.x = element_blank())

ggsave(paste0('3SC_TIMEINTERACTING_WITH_FO_PerSEX_Boxplot', Sys.Date(), ".pdf"),  dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_TIMEINTERACTING_WITH_FO_PerSEX_Boxplot', Sys.Date(), ".jpeg"), dpi = 1000, width = 6, height = 6, units = "in")

ggplot(Expt_Summary_All %>% filter(Expt_Type == 'Sociability',
                                    !(MouseID %in% c('MO1')),
                                    Treat_Geno != 'ESI050'),
       aes(x = Treat_Geno, y = time_NO_sec, color = Treat_Geno, fill = Treat_Geno)) +
  geom_boxplot(position = position_dodge(width = 0.75), alpha = 0.25) +
  geom_jitter(aes(shape = Sex), position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.75),
              size = 5, alpha = 0.7) +
  theme_cowplot() +
  ylab('Total time interacting with novel mouse (seconds)') + xlab("") +
  scale_color_manual(values = group_colors) +
  geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
  scale_fill_manual(values = group_colors) +
  scale_shape_manual(values = sex_shapes) +
  theme(axis.title.x = element_blank())

ggsave(paste0('3SC_TIMEINTERACTING_WITH_NO_PerSEX_Boxplot', Sys.Date(), ".pdf"),  dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_TIMEINTERACTING_WITH_NO_PerSEX_Boxplot', Sys.Date(), ".jpeg"), dpi = 1000, width = 6, height = 6, units = "in")

# ============================================================
# SECTION: TIME INTERACTING (FO and NO) Familiarity
# ============================================================

ggplot(Expt_Summary_All %>% filter(Expt_Type == 'Familiarity',
                                    !(MouseID %in% c('MO1')),
                                    Treat_Geno != 'ESI050'),
       aes(x = Treat_Geno, y = time_FO_sec, color = Treat_Geno, fill = Treat_Geno)) +
  geom_boxplot(position = position_dodge(width = 0.75), alpha = 0.25) +
  geom_jitter(aes(shape = Sex), position = position_jitterdodge(jitter.width = 0.2, dodge.width = 0.75),
              size = 5, alpha = 0.7) +
  theme_cowplot() +
  ylab('Total time interacting with novel object (seconds)') + xlab("") +
  scale_color_manual(values = group_colors) +
  geom_hline(yintercept = 0, linetype = 2, color = 'blue') +
  scale_fill_manual(values = group_colors) +
  scale_shape_manual(values = sex_shapes) +
  theme(axis.title.x = element_blank())

ggsave(paste0('3SC_TIMEINTERACTING_WITH_FO_FAMILIARITY_PerSEX_Boxplot', Sys.Date(), ".pdf"),  dpi = 1000, width = 6, height = 6, units = "in")
ggsave(paste0('3SC_TIMEINTERACTING_WITH_FO_FAMILIARITY_PerSEX_Boxplot', Sys.Date(), ".jpeg"), dpi = 1000, width = 6, height = 6, units = "in")




