library(ggridges)
library(ggplot2)
library(dplyr)
library(tidyr)
library(forcats)
library(tibble)

cholesterol_ihc <- read_excel("NSS_IHC_data_processed.xlsx") # Make a spreadsheet from the raw counts data that provide proportions of the counted neurons that are each score
cholesterol_ihc$Group <- ifelse(cholesterol_ihc$Braak < 3, "HC","AD")

cholesterol_ihc$SN_score <- (cholesterol_ihc$SN_p0)*0 + (cholesterol_ihc$SN_p1)*1 + (cholesterol_ihc$SN_p2)*2 + (cholesterol_ihc$SN_p3)*3
cholesterol_ihc$LC_score <- (cholesterol_ihc$LC_p0)*0 + (cholesterol_ihc$LC_p1)*1 + (cholesterol_ihc$LC_p2)*2 + (cholesterol_ihc$LC_p3)*3

cholesterol_ihc$Marker <- factor(cholesterol_ihc$Marker, levels=c("LDLR","MYLIP","SREBP2","ABCA1","HMGCS1"), ordered = TRUE)

cholesterol_ihc_HC <- cholesterol_ihc[cholesterol_ihc$Group=="HC",]

cholesterol_ihc_AD <- cholesterol_ihc[cholesterol_ihc$Group=="AD",]

# Convert to long format
long_data <- cholesterol_ihc %>%
  pivot_longer(cols = c(SN_score, LC_score),
               names_to = "Region",
               values_to = "AvgScore") %>%
  mutate(Region = gsub("_score", "", Region))

long_data_HC <- long_data[long_data$Group=="HC",]
long_data_AD <- long_data[long_data$Group=="AD",]


# Perform Wilcoxon rank-sum test for each marker comparing SN and LC
significance_results_HC <- cholesterol_ihc_HC %>%
  group_by(Marker) %>%
  summarise(p_value = wilcox.test(SN_score, LC_score, paired=TRUE, exact=FALSE)$p.value)

significance_results_AD <- cholesterol_ihc_AD %>%
  group_by(Marker) %>%
  summarise(p_value = wilcox.test(SN_score, LC_score, paired=TRUE, exact=FALSE)$p.value)

# Convert p-values to significance labels
significance_results_HC <- significance_results_HC %>%
  mutate(significance = case_when(
    p_value < 0.001 ~ "***",
    p_value < 0.01 ~ "**",
    p_value < 0.05 ~ "*",
    TRUE ~ "ns"))

significance_results_AD <- significance_results_AD %>%
  mutate(significance = case_when(
    p_value < 0.001 ~ "***",
    p_value < 0.01 ~ "**",
    p_value < 0.05 ~ "*",
    TRUE ~ "ns"))

# Create a named vector for custom facet labels
marker_labels_HC <- significance_results_HC %>%
  mutate(label = paste0(Marker, " (", significance, ")")) %>%
  select(Marker, label) %>%
  deframe()

marker_labels_AD <- significance_results_AD %>%
  mutate(label = paste0(Marker, " (", significance, ")")) %>%
  select(Marker, label) %>%
  deframe()


region_labels <- c("Locus coeruleus","Substantia nigra")
names(region_labels) <- c("LC","SN")

# Plot with custom facet labels showing significance
ggplot(long_data_HC, aes(x = AvgScore, fill = Region)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 12) +
  geom_density(alpha = 0.9, adjust=2) +
  facet_grid(Region ~ Marker, labeller = labeller(Marker = marker_labels_HC),switch="y") +  # Custom facet labels
  labs(title = "Average Neuron Scores by Region and Protein (Braak stage < III)",
       x = "Average Score",
       y = "Frequency",
       caption = "Significance (Wilcox): *** p < 0.001, ** p < 0.01, * p < 0.05, ns = not significant") + 
  scale_x_continuous(breaks=c(0,1,2,3),labels=c("0","+","++","+++")) +
  scale_fill_manual(values=c("#5857FF","#FFAC0F"))+
  theme_linedraw()+
  theme(legend.position = "none", strip.placement = "outside", strip.text = element_text(face = "bold"))

ggplot(long_data_AD, aes(x = AvgScore, fill = Region)) +
  geom_histogram(position = "identity", alpha = 0.5, bins = 12) +
  geom_density(alpha = 0.9, adjust=2) +
  facet_grid(Region ~ Marker, labeller = labeller(Marker = marker_labels_AD),switch="y") +  # Custom facet labels
  labs(title = "Average Neuron Scores by Region and Protein (Braak stage VI)",
       x = "Average Score",
       y = "Frequency",
       caption = "Significance (Wilcox): *** p < 0.001, ** p < 0.01, * p < 0.05, ns = not significant") + 
  scale_x_continuous(breaks=c(0,1,2,3),labels=c("0","+","++","+++")) +
  scale_fill_manual(values=c("#5857FF","#FFAC0F"))+
  theme_linedraw()+
  theme(legend.position = "none", strip.placement = "outside", strip.text = element_text(face = "bold"))