# This code was used to generate graphs in Nelson, Shoenhard, and Granato: 
# "Pharmacogenetic and whole-brain activity analyses uncover 
# integration of distinct molecular and circuit programs that drive learning"

#This is the latest version as of January 25, 2022



#### Load packages####

library(vroom)
library(janitor)
library(tidyverse)
library(readxl)
library(BBmisc) #contains normalize function
library(cluster) #for agnes
library(factoextra) #for visualizing distance matrices and clusters
library(ggrepel)
library(corrplot)


#### Data import and cleanup####

files <- list.files(path = "~/MATLAB/Jess_Sample_Data_211004/Data",
                    pattern = "*.csv", 
                    full.names = T)
#raw import
data <- vroom(files, 
              id = "path",
              .name_repair = janitor::make_clean_names)

#cleaning and tidying
data_wide <- data %>%
  separate(path, sep = "/", into = c("a","b","c","d","e","f","g","experiment")) %>%
  select(experiment, roi_name, signal_in_roi) %>%
  pivot_wider(names_from = "experiment",
              values_from = "signal_in_roi")

brain_regions <- read_xlsx("293_BrainRegions.xlsx")  #contains long and short brain region names

colnames(brain_regions) <- c("roi_name", "short_roi", "space_roi")

joined_wide <- left_join(brain_regions, data_wide) #adds brain regions that have no signal in any experiment

# converting empty cells to 0s (because that's what they are)
joined_wide_numbers <- joined_wide %>%
  select(-roi_name, -short_roi, -space_roi)

joined_wide_numbers[joined_wide_numbers==""]<-NA

joined_wide_numbers[is.na(joined_wide_numbers)] = 0

joined_wide_numbers$roi_name <- joined_wide$roi_name

joined_wide <- joined_wide_numbers

#write_csv(joined_wide, "RawData_MAP_Mapping_Cleaned.csv")

joined_wide_Neg <- joined_wide %>%
  select("roi_name",contains("Neg")) %>%
  mutate(roi_name = paste(roi_name, " Negative")) %>%
  pivot_longer(!roi_name, names_to = "experiment", values_to = "signal") %>%
  separate(experiment, sep = "_", into = c("treatment", "date", "stimulus","d","e","f","g","h")) %>%
  mutate(stimulus = recode(stimulus, "NoStimuli" = "NoStim")) %>%
  unite(experiment, c(treatment,date,stimulus)) %>%
  select(-d,-e,-f,-g,-h) %>%
  pivot_wider(names_from = "experiment",
              values_from = "signal")

joined_wide_Pos <- joined_wide %>%
  select("roi_name",contains("Pos"))  %>%
  mutate(roi_name = paste(roi_name, " Positive")) %>%
  pivot_longer(!roi_name, names_to = "experiment", values_to = "signal") %>%
  separate(experiment, sep = "_", into = c("treatment", "date", "stimulus","d","e","f","g","h")) %>%
  mutate(stimulus = recode(stimulus, "NoStimuli" = "NoStim")) %>%
  unite(experiment, c(treatment,date,stimulus)) %>%
  select(-d,-e,-f,-g,-h) %>%
  pivot_wider(names_from = "experiment",
              values_from = "signal")

stacked <- rbind(joined_wide_Neg, joined_wide_Pos)  # get negative and positive signal together

norm_try1 <- normalize(stacked, method = "range", range = c(0,1))  # normalize -- all values are positive, so "negative" and "positive" values are all rescaled together

#All columns (including positive and negative data) are now between 0 and 1

#converting to a matrix for clustering purposes
normed_t <- t(as.matrix(norm_try1))

colnames(normed_t) <- norm_try1$roi_name

normed_t_trim <- data.matrix(normed_t[-1,])

#### Canberra Method -- USE THIS ONE ####
distances_canberra <- dist(normed_t_trim, method = "canberra")
glimpse(distances_canberra)
fviz_dist(distances_canberra, 
          gradient = list(low = "blue", mid = NULL, high = "white"), 
          # gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"),
          show_labels = TRUE, 
          lab_size = 8)

?fviz_dist


m <- c( "average", "single", "complete", "ward")
names(m) <- c( "average", "single", "complete", "ward")

# code modified from https://uc-r.github.io/hc_clustering
#using "agnes" -- an alternative agglomerative method to hclust
ac <- function(x) {
  agnes(distances_canberra, method = x)$ac
}

map_dbl(m, ac)

# looks like the best method is Ward!
agnes_ward_clusters_canberra <- agnes(distances_canberra, method = "ward")
plot(agnes_ward_clusters_canberra)

# Major groups: ap2s1, But + MK801, Strych + hip14, kcna, and the Blob
sub_grp_agnes_5_canberra <- cutree(agnes_ward_clusters_canberra, k = 5)
fviz_cluster(list(cluster = sub_grp_agnes_5_canberra, data = distances_canberra, repel = TRUE, axes = 3)) +
  ggtitle("Canberra Distances, Agnes Clustering, 5 groups") +
  theme_bw()

# Separates out the Blob, particularly the MK801 + Butaclamol part
sub_grp_agnes_7_canberra <- cutree(agnes_ward_clusters_canberra, k = 7)
fviz_cluster(list(cluster = sub_grp_agnes_7_canberra, data = distances_canberra)) +
  ggtitle("Canberra Distances, Agnes Clustering, 7 groups") +
  theme_bw()

#Finally separates Strych and hip14
sub_grp_agnes_9_canberra <- cutree(agnes_ward_clusters_canberra, k = 9)
fviz_cluster(list(cluster = sub_grp_agnes_9_canberra, data = distances_canberra)) +
  ggtitle("Canberra Distances, Agnes Clustering, 9 groups") +
  theme_bw()

#### Clustering averages####
normed_t_tib <- as_tibble(normed_t_trim)
normed_t_tib$group <- rownames(normed_t_trim)

average_normed <- normed_t_tib %>%
  separate(group, sep = "_", into = c("treatment", "date", "stimulus")) %>%
  group_by(treatment) %>%
  mutate(across(!c(date,stimulus), as.numeric)) %>%
  summarise(across(!c(date,stimulus),mean))

average_normed_mat <- as.matrix(average_normed)
rownames(average_normed_mat) <- average_normed$treatment
average_normed_mat_trim <- average_normed_mat[,-1]

# Canberra method for distance matrix
dist_AVG_canberra <- dist(average_normed_mat_trim, method = "canberra") 
fviz_dist(dist_AVG_canberra, 
          gradient = list(low = "blue", mid = NULL, high = "white"), 
          show_labels = TRUE, 
          lab_size = 12)

agnes_ward_clusters_avg_canberra <- agnes(dist_AVG_canberra, method = "ward")
plot(agnes_ward_clusters_avg_canberra)

#### Assess correlations between pairs of treatments####

# norm_try1 is from when I normalized the positive and negative values TOGETHER from 0 to 1. EITHER the peak
# positive or the peak negative value should be 1, but not both

pos_only_normed <- norm_try1 %>%    # get positive only
  filter(grepl('Positive', roi_name)) %>%
  select(-roi_name)

#write_csv(pos_only_normed, "Normalized_MAP_Mapping_PosOnly.csv")

neg_only_normed <- norm_try1 %>%    # get negative only - note, these are positive values for now
  filter(grepl('Negative', roi_name)) %>%
  select(-roi_name)

#write_csv(neg_only_normed, "Normalized_MAP_Mapping_NegOnly.csv")

norm_neg_range_for_corr <- pos_only_normed - neg_only_normed  #subtract negative from positive for each region

norm_neg_range_for_corr_tib <- as_tibble(norm_neg_range_for_corr)

norm_neg_range_for_corr_tib$roi_name <- brain_regions$roi_name 

norm_neg_range_for_corr_tib_t <- t(as.matrix(norm_neg_range_for_corr_tib))  #transpose

colnames(norm_neg_range_for_corr_tib_t) <- brain_regions$space_roi

norm_neg_range_for_corr_tib_t <- data.matrix(norm_neg_range_for_corr_tib_t[-51,])

norm_neg_range_for_corr_tib_t_tibble <- as_tibble(norm_neg_range_for_corr_tib_t, 
                                                  rownames = "experiment", 
                                                  .name_repair = "minimal")

# Note: the "blob" is what I have nicknamed the large cluster into which pappa, cacna2d3, kcna1a, and 
# a few of the MK801 and butaclamol samples fall into. I averaged the butaclamol and MK801 samples
# that fell in vs. out of the "blob" and determined that brains without much negative signal landed
# in the "blob" while those with some negative signal generally avoided it. However, the positive signal
# between "blob" and "no blob" butaclamol brains and MK-801 brains was generally correlated. So,
# even though those brains wound up in different clusters, they actually aren't all that different.

for_corr_long <- norm_neg_range_for_corr_tib_t_tibble %>%     #average by condition
  pivot_longer(!experiment, names_to = "region", values_to = "signal") %>%
  pivot_wider(names_from = "experiment", values_from = "signal") %>%
  mutate("but_avg" = ((as.numeric(Butaclamol_210716_Hab) +
                         as.numeric(Butaclamol_210805_Taps) +
                         as.numeric(Butaclamol_210716_NoStim) +
                         as.numeric(Butaclamol_210722_Hab) +
                         as.numeric(Butaclamol_210722_Taps) +
                         as.numeric(Butaclamol_210722_NoStim)) / 6)) %>%
  mutate("but_noblob" = ((as.numeric(Butaclamol_210716_Hab) +
                            as.numeric(Butaclamol_210722_Taps) + 
                            as.numeric(Butaclamol_210722_NoStim)) / 3)) %>%
  mutate("but_blob" = ((as.numeric(Butaclamol_210805_Taps) +
                          as.numeric(Butaclamol_210716_NoStim) +
                          as.numeric(Butaclamol_210722_Hab)) / 3)) %>%                          
  mutate("ap2s1_avg" = ((as.numeric(ap2s1_210624_Hab) +
                           as.numeric(ap2s1_210708_Taps1) +
                           as.numeric(ap2s1_210624_NoStim) +
                           as.numeric(ap2s1_210708_Hab) +
                           as.numeric(ap2s1_210708_Taps2) +
                           as.numeric(ap2s1_210708_NoStim)) / 6)) %>%
  mutate("cacna_avg" = ((as.numeric(cacna2d3_210316_Hab) +
                           as.numeric(cacna2d3_210316_Taps) +
                           as.numeric(cacna2d3_210316_NoStim) +
                           as.numeric(cacna2d3_210603_Hab) +
                           as.numeric(cacna2d3_210617_Taps) +
                           as.numeric(cacna2d3_210603_NoStim)) / 6)) %>%
  mutate("hip14_avg" = ((as.numeric(hip14_210223_Hab) +
                           as.numeric(hip14_210223_Taps) +
                           as.numeric(hip14_210223_NoStim) +
                           as.numeric(hip14_210228_Taps) +
                           as.numeric(hip14_210228_NoStim) +
                           as.numeric(hip14_210527_Hab) +
                           as.numeric(hip14_210527_Taps) +
                           as.numeric(hip14_210527_NoStim)) / 8)) %>%
  mutate("kcna_avg" = ((as.numeric(kcna1a_210309_Hab) +
                          as.numeric(kcna1a_210309_Taps) +
                          as.numeric(kcna1a_210520_NoStim) +
                          as.numeric(kcna1a_210520_Hab) +
                          as.numeric(kcna1a_210520_Taps) +
                          as.numeric(kcna1a_210617_NoStim)) / 6)) %>%
  mutate("mk801_avg" = ((as.numeric(MK801_210701_Hab) +
                           as.numeric(MK801_210701_Taps) +
                           as.numeric(MK801_210701_NoStim) +
                           as.numeric(MK801_210715_Hab) +
                           as.numeric(MK801_210715_Taps) +
                           as.numeric(MK801_210715_NoStim)) / 6)) %>%
  mutate("mk801_blob" = ((as.numeric(MK801_210701_Hab) +
                            as.numeric(MK801_210701_NoStim) +
                            as.numeric(MK801_210715_Taps) +
                            as.numeric(MK801_210715_NoStim)) / 4)) %>%
  mutate("mk801_noblob" = ((as.numeric(MK801_210701_Taps) +
                              as.numeric(MK801_210715_Hab)) / 2)) %>%
  mutate("pappaa_avg" = ((as.numeric(pappaa_210416_Hab) +
                            as.numeric(pappaa_210416_Taps) +
                            as.numeric(pappaa_210416_NoStim) +
                            as.numeric(pappaa_210610_Hab) +
                            as.numeric(pappaa_210610_Taps) +
                            as.numeric(pappaa_210610_NoStim)) / 6)) %>%
  mutate("strych_avg" = ((as.numeric(Strychnine_210729_Hab) +
                            as.numeric(Strychnine_210729_Taps) +
                            as.numeric(Strychnine_210729_NoStim) +
                            as.numeric(Strychnine_210809_Hab) +
                            as.numeric(Strychnine_210809_Taps) +
                            as.numeric(Strychnine_210809_NoStim)) / 6))

#write_csv(for_corr_long, "Normalized_MAP_Mapping_Data_PosNegSUMMED.csv")

averages_only <- for_corr_long %>%
  select(region,
         ap2s1_avg, 
         but_avg,
         but_noblob,
         but_blob,
         cacna_avg, 
         hip14_avg, 
         kcna_avg, 
         mk801_avg,
         mk801_noblob,
         mk801_blob,
         pappaa_avg, 
         strych_avg)

averages_only$reg_cat <- averages_only$region   

averages_only1 <- averages_only  %>%
  separate(reg_cat, into = c("reg_big", "drop", sep = " "))

#### Graphs of correlations####

#### Butaclamol vs ap2s1####
but_ap2s1 <- averages_only1 %>%
  select(but_avg, ap2s1_avg, reg_big, region) %>%
  mutate("quadrant" = case_when(
    ((but_avg>0.05 & ap2s1_avg >0.05)~ "top_right"),
    ((but_avg< -0.05 & ap2s1_avg > 0.05) ~ "top_left"),
    ((but_avg < -0.05 & ap2s1_avg < -0.05) ~ "bottom_left"),
    ((but_avg > 0.05 & ap2s1_avg < -0.05) ~ "bottom_right"),
    (TRUE ~ "toosmall") 
  )
  ) %>%
  mutate("total_abs" = abs(but_avg) + abs(ap2s1_avg))
#  mutate("label_logical" = ifelse(
#    test = (abs(ap2s1_avg) + abs(but_avg) > .55), 
#   yes = as.character(region),
#    no = ""
# ))

but_ap2s1_extremes <- but_ap2s1 %>%
  filter(quadrant != "toosmall") %>%
  group_by(quadrant) %>%
  slice_max(order_by = total_abs, n = 5) %>%
  pull(region)

but_ap2s1_extremes

but_ap2s1_all <- but_ap2s1 %>%
  mutate(reg_big = ifelse(region %in% but_ap2s1_extremes, yes = "extreme", no = reg_big)) %>%
  mutate(reg_big = as.factor(reg_big))

ggplot(but_ap2s1_all, aes(x = but_avg, y = ap2s1_avg, color = reg_big)) +
  geom_point(size = 2, shape = 17) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_text_repel(aes(label=ifelse(region %in% but_ap2s1_extremes, as.character(region), "")),
                  max.overlaps = 32,   
                  hjust = 0,
                  vjust=0, 
                  show.legend = FALSE) +
  #  geom_text(aes(label = as.character(region))) +
  ggtitle("Ap2s1 vs. Butaclamol") +
  labs(x = "Butaclamol Average Signal", y = "Ap2s1 Average Signal", color = "Brain Region") +
  scale_color_manual(
    values = c("#F8766D",
               "#000000",
               "#C49A00",
               "#53B400",
               #  "#00C094",
               "#00B6EB",
               "#A58AFF",
               "#FB61D7"),
    labels = c("Diencephalon", 
               "Label", 
               "Ganglia", 
               "Mesencephalon", 
               "Rhombencephalon", 
               "Spinal Cord", 
               "Telencephalon")) +
  scale_x_continuous(limits = c(-1,1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-1,1), expand = c(0, 0)) +
  # scale_shape_manual(values = c(16)) +
  coord_fixed() +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        plot.title = element_text(size = 18))

#### Mk801 vs Butaclamol####
but_mk801 <- averages_only1 %>%
  select(but_avg, mk801_avg, reg_big, region) %>%
  mutate("quadrant" = case_when(
    ((but_avg>0.05 & mk801_avg >0.05)~ "top_right"),
    ((but_avg< -0.05 & mk801_avg > 0.05) ~ "top_left"),
    ((but_avg < -0.05 & mk801_avg < -0.05) ~ "bottom_left"),
    ((but_avg > 0.05 & mk801_avg < -0.05) ~ "bottom_right"),
    (TRUE ~ "toosmall") 
  )
  ) %>%
  mutate("total_abs" = abs(but_avg) + abs(mk801_avg))
#  mutate("label_logical" = ifelse(
#    test = (abs(ap2s1_avg) + abs(but_avg) > .55), 
#   yes = as.character(region),
#    no = ""
# ))

but_mk801_extremes <- but_mk801 %>%
  filter(quadrant != "toosmall") %>%
  group_by(quadrant) %>%
  slice_max(order_by = total_abs, n = 5) %>%
  pull(region)

but_mk801_extremes

but_mk801_all <- but_mk801 %>%
  mutate(reg_big = ifelse(region %in% but_mk801_extremes, yes = "extreme", no = reg_big)) %>%
  mutate(reg_big = as.factor(reg_big))

ggplot(but_mk801_all, aes(x = but_avg, y = mk801_avg, color = reg_big)) +
  geom_point(size = 2, shape = 17) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_text_repel(aes(label=ifelse(region %in% but_mk801_extremes, as.character(region), "")),
                  max.overlaps = 45,   
                  hjust = 0,
                  vjust=0, 
                  show.legend = FALSE) +
  #  geom_text(aes(label = as.character(region))) +
  ggtitle("Butaclamol vs. MK801") +
  labs(x = "Butaclamol Average Signal", y = "MK801 Average Signal", color = "Brain Region") +
  scale_color_manual(
    values = c("#F8766D",
               "#000000",
               "#C49A00",
               "#53B400",
               #  "#00C094",
               "#00B6EB",
               "#A58AFF",
               "#FB61D7"),
    labels = c("Diencephalon", 
               "Label", 
               "Ganglia", 
               "Mesencephalon", 
               "Rhombencephalon", 
               "Spinal Cord", 
               "Telencephalon")) +
  scale_x_continuous(limits = c(-1,1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-1,1), expand = c(0, 0)) +
  # scale_shape_manual(values = c(16)) +
  coord_fixed() +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        plot.title = element_text(size = 18))


#### Butaclamol vs Pappaa ####
but_pappaa <- averages_only1 %>%
  select(but_avg, pappaa_avg, reg_big, region) %>%
  mutate("quadrant" = case_when(
    ((but_avg>0.05 & pappaa_avg >0.05)~ "top_right"),
    ((but_avg< -0.05 & pappaa_avg > 0.05) ~ "top_left"),
    ((but_avg < -0.05 & pappaa_avg < -0.05) ~ "bottom_left"),
    ((but_avg > 0.05 & pappaa_avg < -0.05) ~ "bottom_right"),
    (TRUE ~ "toosmall") 
  )
  ) %>%
  mutate("total_abs" = abs(but_avg) + abs(pappaa_avg))
#  mutate("label_logical" = ifelse(
#    test = (abs(ap2s1_avg) + abs(but_avg) > .55), 
#   yes = as.character(region),
#    no = ""
# ))

but_pappaa_extremes <- but_pappaa %>%
  filter(quadrant != "toosmall") %>%
  group_by(quadrant) %>%
  slice_max(order_by = total_abs, n = 5) %>%
  pull(region)

but_pappaa_extremes

but_pappaa_all <- but_pappaa %>%
  mutate(reg_big = ifelse(region %in% but_pappaa_extremes, yes = "extreme", no = reg_big)) %>%
  mutate(reg_big = as.factor(reg_big))

ggplot(but_pappaa_all, aes(x = but_avg, y = pappaa_avg, color = reg_big)) +
  geom_point(size = 2, shape = 17) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_text_repel(aes(label=ifelse(region %in% but_pappaa_extremes, as.character(region), "")),
                  max.overlaps = 300,   
                  hjust = 0,
                  vjust=0, 
                  show.legend = FALSE) +
  #  geom_text(aes(label = as.character(region))) +
  ggtitle("Butaclamol vs. pappaa") +
  labs(x = "Butaclamol Average Signal", y = "pappaa Average Signal", color = "Brain Region") +
  scale_color_manual(
    values = c("#F8766D",
               "#000000",
               "#C49A00",
               "#53B400",
               #  "#00C094",
               "#00B6EB",
               "#A58AFF",
               "#FB61D7"),
    labels = c("Diencephalon", 
               "Label", 
               "Ganglia", 
               "Mesencephalon", 
               "Rhombencephalon", 
               "Spinal Cord", 
               "Telencephalon")) +
  scale_x_continuous(limits = c(-1,1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-1,1), expand = c(0, 0)) +
  # scale_shape_manual(values = c(16)) +
  coord_fixed() +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        plot.title = element_text(size = 18))

#### Hip14 vs Strychnine ####

hip14_strych <- averages_only1 %>%
  select(strych_avg, hip14_avg, reg_big, region) %>%
  mutate("quadrant" = case_when(
    ((strych_avg>0.05 & hip14_avg >0.05)~ "top_right"),
    ((strych_avg< -0.05 & hip14_avg > 0.05) ~ "top_left"),
    ((strych_avg < -0.05 & hip14_avg < -0.05) ~ "bottom_left"),
    ((strych_avg > 0.05 & hip14_avg < -0.05) ~ "bottom_right"),
    (TRUE ~ "toosmall") 
  )
  ) %>%
  mutate("total_abs" = abs(strych_avg) + abs(hip14_avg))
#  mutate("label_logical" = ifelse(
#    test = (abs(ap2s1_avg) + abs(but_avg) > .55), 
#   yes = as.character(region),
#    no = ""
# ))

hip14_strych_extremes <- hip14_strych %>%
  filter(quadrant != "toosmall") %>%
  group_by(quadrant) %>%
  slice_max(order_by = total_abs, n = 10) %>%
  pull(region)

hip14_strych_extremes

hip14_strych_all <- hip14_strych %>%
  mutate(reg_big = ifelse(region %in% hip14_strych_extremes, yes = "extreme", no = reg_big)) %>%
  mutate(reg_big = as.factor(reg_big))

ggplot(hip14_strych_all, aes(x = hip14_avg, y = strych_avg, color = reg_big)) +
  geom_point(size = 2, shape = 17) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_text_repel(aes(label=ifelse(region %in% hip14_strych_extremes, as.character(region), "")),
                  max.overlaps = 300,   
                  hjust = 0,
                  vjust=0, 
                  show.legend = FALSE) +
  #  geom_text(aes(label = as.character(region))) +
  ggtitle("hip14 vs. Strychnine") +
  labs(x = "hip14 Average Signal", y = "Strychnine Average Signal", color = "Brain Region") +
  scale_color_manual(
    values = c("#F8766D",
               "#000000",
               "#C49A00",
               "#53B400",
               #  "#00C094",
               "#00B6EB",
               "#A58AFF",
               "#FB61D7"),
    labels = c("Diencephalon", 
               "Label", 
               "Ganglia", 
               "Mesencephalon", 
               "Rhombencephalon", 
               "Spinal Cord", 
               "Telencephalon")) +
  scale_x_continuous(limits = c(-1,1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-1,1), expand = c(0, 0)) +
  # scale_shape_manual(values = c(16)) +
  coord_fixed() +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        plot.title = element_text(size = 18))

#### Hip14 vs kcna1a####

hip14_kcna <- averages_only1 %>%
  select(kcna_avg, hip14_avg, reg_big, region) %>%
  mutate("quadrant" = case_when(
    ((kcna_avg>0.05 & hip14_avg >0.05)~ "top_right"),
    ((kcna_avg< -0.05 & hip14_avg > 0.05) ~ "top_left"),
    ((kcna_avg < -0.05 & hip14_avg < -0.05) ~ "bottom_left"),
    ((kcna_avg > 0.05 & hip14_avg < -0.05) ~ "bottom_right"),
    (TRUE ~ "toosmall") 
  )
  ) %>%
  mutate("total_abs" = abs(kcna_avg) + abs(hip14_avg))
#  mutate("label_logical" = ifelse(
#    test = (abs(ap2s1_avg) + abs(but_avg) > .55), 
#   yes = as.character(region),
#    no = ""
# ))

hip14_kcna_extremes <- hip14_kcna %>%
  filter(quadrant != "toosmall") %>%
  group_by(quadrant) %>%
  slice_max(order_by = total_abs, n = 9) %>%
  pull(region)

hip14_kcna_extremes

hip14_kcna_all <- hip14_kcna %>%
  mutate(reg_big = ifelse(region %in% hip14_kcna_extremes, yes = "extreme", no = reg_big)) %>%
  mutate(reg_big = as.factor(reg_big))

ggplot(hip14_kcna_all, aes(x = hip14_avg, y = kcna_avg, color = reg_big)) +
  geom_point(size = 2, shape = 17) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_text_repel(aes(label=ifelse(region %in% hip14_kcna_extremes, as.character(region), "")),
                  max.overlaps = 300,   
                  hjust = 0,
                  vjust=0, 
                  show.legend = FALSE) +
  #  geom_text(aes(label = as.character(region))) +
  ggtitle("hip14 vs. kcna1a") +
  labs(x = "hip14 Average Signal", y = "kcna1a Average Signal", color = "Brain Region") +
  scale_color_manual(
    values = c("#F8766D",
               "#000000",
               "#C49A00",
               "#53B400",
               #  "#00C094",
               "#00B6EB",
               "#A58AFF",
               "#FB61D7"),
    labels = c("Diencephalon", 
               "Label", 
               "Ganglia", 
               "Mesencephalon", 
               "Rhombencephalon", 
               "Spinal Cord", 
               "Telencephalon")) +
  scale_x_continuous(limits = c(-1,1), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-1,1), expand = c(0, 0)) +
  # scale_shape_manual(values = c(16)) +
  coord_fixed() +
  theme_bw() +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 12),
        plot.title = element_text(size = 18))

#### Correlations####

strict_avg <- averages_only %>%
  select(-but_noblob, -but_blob, -mk801_noblob, -mk801_blob, -region, -reg_cat) 

strict_avg_num <- data.matrix(strict_avg)

rownames(strict_avg_num) <- averages_only$region

correlation_averages <- cor(strict_avg_num)
testRes = cor.mtest(strict_avg_num, conf.level = 0.95)

colnames(correlation_averages) <- c("ap2s1", "Butaclamol", "cacna2d3", "hip14", "kcna1a", "MK801", "pappaa", "Strychnine")
rownames(correlation_averages) <- c("ap2s1", "Butaclamol", "cacna2d3", "hip14", "kcna1a", "MK801", "pappaa", "Strychnine")

corrplot(correlation_averages, 
         method = 'number', 
         tl.col = "black", 
         order = "AOE", 
         p.mat = testRes$p, 
         sig.level = 0.05,
         insig = "blank")
corrplot(correlation_averages, method = 'ellipse', order = 'AOE', type = 'upper',tl.col = "black")
corrplot(correlation_averages, order = 'AOE', addCoef.col = 'black', tl.pos = 'd',
         cl.pos = 'n')

corrplot(correlation_averages, 
         addCoef.col ='black',
         type = "upper",
         tl.col = "black",
         method = "ellipse", 
         p.mat = testRes$p, sig.level = 0.05, 
         order = 'hclust', 
         insig = "blank")


#### Module identification and export for Matlab####

#MK801, Butaclamol, and Hip14 all upreg or all downreg
Module1 <- averages_only1 %>%
  select(but_avg, mk801_avg, hip14_avg, reg_big, region) %>%
  mutate("quadrant" = case_when(
    (((but_avg>0.05 & mk801_avg >0.05) & hip14_avg >0.05)~ "upreg"),
    (((but_avg < -0.05 & mk801_avg < -0.05) & hip14_avg < -0.05) ~ "downreg"),
    (TRUE ~ "toosmall") 
  )
  ) %>%
  mutate("total_abs" = abs(but_avg) + abs(mk801_avg) + abs(hip14_avg))# %>%
 # filter(quadrant != "toosmall")

write_csv(Module1, "Module_1_Regions.csv")

# Module 2: hip14 + kcna
Module2 <- averages_only1 %>%
  select(kcna_avg, hip14_avg, reg_big, region) %>%
  mutate("quadrant" = case_when(
    ((kcna_avg > 0.05 & hip14_avg > 0.05)~ "upreg"),
    ((kcna_avg < -0.05 & hip14_avg < -0.05) ~ "downreg"),
    (TRUE ~ "toosmall") 
  )
  ) %>%
  mutate("total_abs" = abs(kcna_avg) + abs(hip14_avg)) %>%
  filter(quadrant == "upreg")

write_csv(Module2, "Module_2_Regions.csv")

# Module 3: Opposite in Ap2s1+Pappaa and Butaclamol
Module3 <- averages_only1 %>%
  select(ap2s1_avg, pappaa_avg, but_avg, reg_big, region) %>%
  mutate("quadrant" = case_when(
    (((ap2s1_avg > 0.05) & but_avg < -0.05)~ "Ap2Up_ButDown"),
    (((pappaa_avg > 0.05) & but_avg < -0.05)~ "PapUp_ButDown"),
    (((ap2s1_avg < -0.05) & but_avg > 0.05) ~ "Ap2DownButUp"),
    (((pappaa_avg < -0.05) & but_avg > 0.05) ~ "PapDownButUp"),
    (TRUE ~ "toosmall") 
  )
  ) %>%
  mutate("total_abs" = abs(ap2s1_avg) + abs(pappaa_avg) + abs(but_avg)) %>%
  filter(quadrant != "toosmall")

Module3_top20 <- Module3 %>%
  slice_max(order_by = total_abs, n=20)

write_csv(Module3, "Module_3_Regions.csv")



Module1_ranked <- Module1_try2 %>%
  filter(quadrant != "toosmall") %>%
  group_by(quadrant) %>%
  slice_max(order_by = total_abs, n = 10) %>%
  pull(region, total_abs)

Module1_ranked