added GC analysis script for implementation

workflows/GC_MS_analysis/210903_primary_analysis_cmQTL_val_1_levene_RxNV.R

+rm(list = ls())
+# Set working directory to source file location ---------------------------
+
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+getwd()
+
+library(openxlsx)
+library(tidyverse)
+library(car)
+library(pheatmap)
+library(broom)
+library(ggpubr)
+library(viridisLite)
+library(ggtext)
+library(glue)
+
+current <- getwd()
+source <- str_c(current,"/..")
+
+cur_date <- str_c(str_replace_all(Sys.Date(),"^.{2}|-",""))
+
+out <- str_c(cur_date, "analysis", sep = "_")
+
+if (file.exists(out)) {
+  cat("The folder already exists")
+} else {
+  dir.create(out)
+}
+
+out_dir <- str_c(current, out, sep = "/")
+
+# Data loading ------------------------------------------------------------
+
+exclude <- read_csv("compound_rts_exclusion.csv")
+
+latest <- str_sort(str_extract(list.files(pattern = "^\\d{6}_normalization$"),
+                               pattern = "^\\d{6}_normalization"),
+                   decreasing = T)[[1]]
+
+latest_norm <- str_c(current, "/", latest)
+
+setwd(latest_norm)
+
+latest_norm_date <- str_extract(latest, pattern = "^\\d{6}")
+
+sam_dat <- read_csv(str_c(latest_norm_date, "_cmQTL_val_1_2_sam_dat_GC.csv"))
+met_dat <- read_csv(str_c(latest_norm_date, "_cmQTL_val_1_2_met_dat_GC.csv"), col_types = "ccfdf") %>% 
+  filter(str_detect(Compound_Name,"ribitol", negate = T))
+feat_dat <- read_csv(str_c(latest_norm_date, "_cmQTL_val_1_2_feat_dat_GC.csv"), col_types = "fcffffTiiffidfffcfdddd")
+
+setwd(out_dir)
+
+# Data combination --------------------------------------------------------
+
+ met_dat<- met_dat %>% 
+  left_join(exclude) %>% 
+  filter(!exclude)
+
+mets <- met_dat$met
+
+genotypes <- sam_dat %>% distinct(genotype,.keep_all = T) %>% 
+  select(genotype, alias)
+
+GC_classes <- feat_dat %>% 
+  distinct(Compound_Class, Compound_Name, met)
+
+GC_long <- feat_dat %>% 
+  left_join(met_dat) %>% 
+  left_join(exclude) %>% 
+  filter(exp == 1, class == "sample", !exclude)
+
+
+
+# Means -------------------------------------------------------------------
+
+means_tec_rep <- GC_long %>% 
+  group_by(plantline, alias, genotype, treatment, met, tissue, LIMS_ID) %>% 
+  summarise(mean_tec_rep = mean(rescaled)) %>% 
+  ungroup()
+
+means <- means_tec_rep %>% 
+  group_by(met, tissue, treatment, genotype, alias) %>% 
+  summarise(mean = mean(mean_tec_rep),
+            sd = sd(mean_tec_rep),
+            n = n()) %>% 
+  ungroup() %>% 
+  mutate(se = sd/sqrt(n))
+
+miss_per_treat <- GC_long %>% 
+  group_by(met, tissue, treatment) %>% 
+  summarise(na = sum(is.na(area)),
+            n_treat = n()) %>% 
+  mutate(percent_na = na/n_treat*100) %>% 
+  ungroup()
+
+
+# Significance analysis ---------------------------------------------------
+
+
+
+GC_tidy <- means_tec_rep %>% 
+  pivot_wider(id_cols = c( genotype, alias, treatment, tissue, LIMS_ID),
+              names_from = met,
+              values_from = mean_tec_rep)
+
+GC_tidy_numeric <- GC_tidy %>% 
+  select(all_of(mets))
+
+GC_aov <- map(.x = GC_tidy_numeric, .f = ~aov(.x ~ alias * treatment*tissue, data = GC_tidy))
+
+
+GC_tuk <- map(.x = GC_aov, .f = ~TukeyHSD(.x)) %>% 
+  map(.f = tidy) %>% 
+  map2(.y = names(GC_aov), .f = ~.x %>% mutate(var = .y)) %>% 
+  purrr::reduce(bind_rows)
+
+sig_GC_groups <- GC_tuk %>% 
+  filter(term == "alias:treatment:tissue") %>% 
+  separate(col = contrast, into = c("group1", "group2"), sep = "-")
+
+sig_GC <- sig_GC_groups %>%  
+  separate(group1, into = c("alias1", "treatment1", "tissue1"), sep = ":") %>% 
+  separate(group2, into = c("alias2", "treatment2", "tissue2"), sep = ":") %>% 
+  filter(treatment1 == treatment2, tissue1 == tissue2) %>% 
+  filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT") %>% 
+  ungroup() %>% 
+  mutate(p.signif = if_else(adj.p.value <= 0.05, "*", "ns")) %>% 
+  left_join(means, by = c("var" = "met", "treatment1" = "treatment", "alias1" = "alias", "tissue1" = "tissue")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1 = mean, se1 = se) %>%
+  left_join(means, by = c("var" = "met", "treatment2" = "treatment", "alias2" = "alias", "tissue2" = "tissue")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1, mean2 = mean, se1, se2 = se) %>% 
+  group_by(treatment1, tissue1, var) %>% 
+  mutate(tot_val1 = mean1 + se1,
+         tot_val2 = mean2 + se2,
+         y.position = 1.1*(max(tot_val1, tot_val2))) %>% 
+  rename(met = var,
+         group1 = alias1,
+         group2 = alias2,
+         treatment = treatment1,
+         tissue = tissue1) %>% 
+  mutate(met = as_factor(met),
+         group1 = as_factor(group1),
+         group2 = as_factor(group2),
+         treatment = as_factor(treatment),
+         tissue = as_factor(tissue)) %>% 
+  ungroup()
+
+sig_mets <- sig_GC %>% 
+  filter(p.signif == "*") %>% 
+  distinct(met) %>% 
+  mutate(sig = T)
+
+
+# Significance with t-test ------------------------------------------------
+
+GC_tidy <- means_tec_rep %>% 
+  pivot_wider(id_cols = c( genotype, alias, treatment, tissue, LIMS_ID),
+              names_from = met,
+              values_from = mean_tec_rep) %>% 
+  mutate(group = str_c(tissue, treatment, alias, sep = "_"))
+
+GC_tidy_numeric <- GC_tidy %>% 
+  select(all_of(mets))
+
+GC_t <- map(.x = GC_tidy_numeric,
+             .f = ~pairwise.t.test(x = .x ,
+                                   g = GC_tidy$group,
+                                   p.adjust.method = "none",
+                                   pool.sd = F)) %>% 
+  map(.f = tidy)
+
+GC_t_tidy <- GC_t %>% 
+  map2(.y = names(GC_t), .f = ~.x %>% mutate(var = .y)) %>% 
+  purrr::reduce(bind_rows)
+
+sig_GC<- GC_t_tidy %>% 
+  separate(group1, into = c("tissue1", "treatment1", "alias1"), sep = "_") %>% 
+  separate(group2, into = c("tissue2", "treatment2", "alias2"), sep = "_") %>% 
+  filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT",
+         tissue1 == tissue2, treatment1 == treatment2) %>% 
+  group_by(var) %>% 
+  mutate(adj.p.value = p.adjust(p.value)) %>% 
+  #mutate(adj.p.value = p.value * 121) %>% 
+  ungroup() %>% 
+  mutate(p.signif = if_else(p.value <= 0.05, "*", "ns")) %>% 
+  left_join(means, by = c("var" = "met", "treatment1" = "treatment", "alias1" = "alias", "tissue1" = "tissue")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1 = mean, se1 = se, p.value, adj.p.value) %>%
+  left_join(means, by = c("var" = "met", "treatment2" = "treatment", "alias2" = "alias", "tissue2" = "tissue")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1, mean2 = mean, se1, se2 = se, p.value, adj.p.value) %>% 
+  group_by(treatment1, tissue1, var) %>% 
+  mutate(tot_val1 = mean1 + se1,
+         tot_val2 = mean2 + se2,
+         y.position = 1.1*(max(tot_val1, tot_val2))) %>% 
+  ungroup() %>% 
+  rename(met = var,
+         group1 = alias1,
+         group2 = alias2,
+         treatment = treatment1,
+         tissue = tissue1) %>% 
+  left_join(miss_per_treat) %>% 
+  mutate(met = as_factor(met),
+         group1 = as_factor(group1),
+         group2 = as_factor(group2),
+         treatment = as_factor(treatment),
+         tissue = as_factor(tissue),
+         p.signif = if_else(percent_na >= 60, "ns", p.signif))
+  
+
+sig_mets <- sig_GC %>% 
+  filter(p.signif == "*") %>% 
+  distinct(met) %>% 
+  mutate(sig = T)
+
+
+
+# Heatmap scaled all-----------------------------------------------------------------
+
+heat.GC <- means %>% 
+  group_by(met) %>% 
+  mutate(se = sd/sqrt(n),
+         total_norm = mean/mean(mean),
+         log_norm = log2(total_norm),
+         log_norm = if_else(is.na(log_norm)| is.infinite(log_norm), 0, log_norm),
+         log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm)),
+         log_norm_level = if_else(is.na(log_norm_level), 0, log_norm_level)) %>% 
+  left_join(sig_mets) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>% 
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = log_norm_level) %>% 
+  left_join(met_dat) %>%
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>% 
+  select(contains("leaves"), contains("fruits")) %>% as.matrix()
+
+annotation_row <- heat.GC %>% 
+  select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- sam_dat %>% 
+  filter(exp == 1, class == "sample") %>%
+  distinct(tissue, treatment, genotype) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  as.data.frame()
+
+rownames(annotation_col) <- annotation_col$group
+
+annotation_col <- annotation_col %>% 
+  select(-group, treatment, tissue, genotype)
+
+heat.GC_signif <- means %>% 
+  distinct(tissue, treatment, alias, genotype, met) %>% 
+  #filter(genotype != "MoneyMaker") %>% 
+  left_join(sig_GC, by = c("met", "tissue", "treatment", "alias" = "group2")) %>% 
+  left_join(sig_mets) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")),
+         signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = signif) %>% 
+  left_join(met_dat) %>% 
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+heat_cols <- colnames(mat.heat.GC)
+
+rownames(heat.GC_signif) <- heat.GC_signif$met
+
+mat.heat.GC_signif <- heat.GC_signif %>% 
+  select(all_of(heat_cols)) %>% as.matrix()
+
+ann_colors = list(
+  tissue = c(fruits = "red",leaves = "darkgreen"),
+  genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+  treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+
+pheatmap.GC <- pheatmap(mat.heat.GC,
+                         #color = plasma(14),
+                         #cellwidth = 9,
+                         #cellheight = 6,
+                         #breaks = c(-6.5,-5.5,-4.5,-3.5,-2.5,-1.5, -0.5 ,0.5,1.5,2.5,3.5,4.5,5.5,6.5),
+                         cluster_rows = T,
+                         cluster_cols = T,
+                         annotation_names_row = F,
+                         show_rownames = F,
+                         annotation_row = annotation_row,
+                         annotation_col = annotation_col,
+                         display_numbers = mat.heat.GC_signif,
+                         number_color = "black",
+                         fontsize_number = 6,
+                         annotation_colors = ann_colors,
+                        angle_col = 45,
+                        fontsize_col = 6,
+                         filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+                                          "cmQTL_val1_heatmap_rel_all.jpg",
+                                          sep = "_")
+)
+
+# Heatmap scaled per tissue-----------------------------------------------------------------
+
+heat.GC <- means%>% 
+  group_by(tissue, met) %>% 
+  mutate(total_norm = mean/mean(mean),
+         log_norm = log2(total_norm),
+         log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+         log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>% 
+  ungroup() %>% 
+  left_join(sig_mets) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>% 
+  pivot_wider(id_cols = c( met),
+              names_from = group,
+              values_from = log_norm_level) %>% 
+  left_join(met_dat) %>%
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>% 
+  select(contains("leaves"), contains("fruits")) %>% as.matrix()
+
+annotation_row <- heat.GC %>% 
+  select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- sam_dat %>% 
+  filter(exp == 1, class == "sample") %>%
+  distinct(tissue, treatment, genotype) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  as.data.frame()
+
+rownames(annotation_col) <- annotation_col$group
+
+annotation_col <- annotation_col %>% 
+  select(-group, treatment, tissue, genotype)
+
+heat.GC_signif <- means %>% 
+  distinct(tissue, treatment, alias, genotype, met) %>% 
+  #filter(genotype != "MoneyMaker") %>% 
+  left_join(sig_GC, by = c("met", "tissue", "treatment", "alias" = "group2")) %>% 
+  left_join(sig_mets) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")),
+         signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+  pivot_wider(id_cols = c( met),
+              names_from = group,
+              values_from = signif) %>% 
+  left_join(met_dat) %>% 
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+heat_cols <- colnames(mat.heat.GC)
+
+rownames(heat.GC_signif) <- heat.GC_signif$met
+
+mat.heat.GC_signif <- heat.GC_signif %>% 
+  select(all_of(heat_cols)) %>% as.matrix()
+
+ann_colors = list(
+  tissue = c(fruits = "red",leaves = "darkgreen"),
+  genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+  treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+
+
+pheatmap.GC <- pheatmap(mat.heat.GC,
+                         #color = plasma(14),
+                         #cellwidth = 16,
+                         #cellheight = 4,
+                         #breaks = c(-6.5,-5.5,-4.5,-3.5,-2.5,-1.5, -0.5 ,0.5,1.5,2.5,3.5,4.5,5.5,6.5),
+                         cluster_rows = T,
+                         cluster_cols = T,
+                         annotation_names_row = F,
+                         show_rownames = F,
+                         annotation_row = annotation_row,
+                         annotation_col = annotation_col,
+                         display_numbers = mat.heat.GC_signif,
+                         number_color = "black",
+                         fontsize_number = 6,
+                        angle_col = 45,
+                        fontsize_col = 6,
+                         annotation_colors = ann_colors,
+                         filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+                                          "cmQTL_val1_heatmap_rel_tissue.jpg",
+                                          sep = "_")
+)
+
+
+# Heatmap scaled per tissue and treatment-----------------------------------------------------------------
+
+heat.GC <- means%>% 
+  group_by(tissue, met, treatment) %>% 
+  mutate(total_norm = mean/mean(mean),
+         log_norm = log2(total_norm),
+         log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+         log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>% 
+  ungroup() %>% 
+  left_join(sig_mets) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>% 
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = log_norm_level) %>% 
+  left_join(met_dat) %>%
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>% 
+  select(contains("leaves"), contains("fruits")) %>% as.matrix()
+
+annotation_row <- heat.GC %>% 
+  select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- sam_dat %>% 
+  filter(exp == 1, class == "sample") %>%
+  distinct(tissue, treatment, genotype) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  as.data.frame()
+
+rownames(annotation_col) <- annotation_col$group
+
+annotation_col <- annotation_col %>% 
+  select(-group, treatment, tissue, genotype)
+
+heat.GC_signif <- means %>% 
+  distinct(tissue, treatment, alias, genotype, met) %>% 
+  #filter(genotype != "MoneyMaker") %>% 
+  left_join(sig_GC, by = c("met", "tissue", "treatment", "alias" = "group2")) %>% 
+  left_join(sig_mets) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")),
+         signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+  pivot_wider(id_cols = c( met),
+              names_from = group,
+              values_from = signif) %>% 
+  left_join(met_dat) %>% 
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+heat_cols <- colnames(mat.heat.GC)
+
+rownames(heat.GC_signif) <- heat.GC_signif$met
+
+mat.heat.GC_signif <- heat.GC_signif %>% 
+  select(all_of(heat_cols)) %>% as.matrix()
+
+ann_colors = list(
+  tissue = c(fruits = "red",leaves = "darkgreen"),
+  genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+  treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+
+
+pheatmap.GC <- pheatmap(mat.heat.GC,
+                         #color = plasma(14),
+                         #cellwidth = 16,
+                         #cellheight = 4,
+                         #breaks = c(-6.5,-5.5,-4.5,-3.5,-2.5,-1.5, -0.5 ,0.5,1.5,2.5,3.5,4.5,5.5,6.5),
+                         cluster_rows = T,
+                         cluster_cols = T,
+                         annotation_names_row = F,
+                         show_rownames = F,
+                         annotation_row = annotation_row,
+                         annotation_col = annotation_col,
+                         display_numbers = mat.heat.GC_signif,
+                         number_color = "black",
+                         fontsize_number = 6,
+                        angle_col = 45,
+                        fontsize_col = 6,
+                         annotation_colors = ann_colors,
+                         filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+                                          "cmQTL_val1_heatmap_rel_tissue_treatment.jpg",
+                                          sep = "_")
+)
+
+
+
+# Per metabolite comparisons ----------------------------------------------
+
+norm_MM <- means %>% 
+  filter(genotype == "MoneyMaker") %>% 
+  select(tissue, treatment, met, MM_mean = mean)
+
+
+
+fc <- means_tec_rep %>% 
+  left_join(norm_MM) %>% 
+  mutate(fc = mean_tec_rep/MM_mean) %>% 
+  group_by(tissue, treatment, alias, genotype, met) %>% 
+  summarise(mean_fc = mean(fc),
+            sd = sd(fc),
+            n = n()) %>% 
+  mutate(se = sd/sqrt(n),
+         group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  left_join(sig_GC, by = c("tissue", "treatment", "met", "alias" = "group2")) %>% 
+  group_by(met) %>% 
+  filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>% 
+  left_join(met_dat) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  ungroup() %>% 
+  mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+         genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+plotmets <- fc %>% distinct(met) %>% as_vector()
+plottissues <- fc %>% distinct(tissue) %>% as_vector() %>% as.character
+
+labelnames <- plotmets %>% as_tibble() %>% 
+  rename(met = value) %>% 
+  left_join(met_dat) %>% left_join(GC_classes) %>%
+  select(Compound_Name) %>% as_vector()
+
+
+plot_out <- vector("list", length = length(plotmets))
+per_comp_y <- fc %>%
+  group_by(tissue, treatment, met) %>% 
+  summarise(mean = max(mean_fc),
+            se = max(se))
+
+for(tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    
+    
+    sig_bar <- fc %>% 
+      group_by(tissue, treatment, met) %>% 
+      mutate(tot_val = max(mean_fc + se)) %>%
+      mutate(y.position = tot_val + 0.25*max(tot_val)) %>%
+      ungroup() %>% 
+      rename(genotype1 = genotype) %>% 
+      left_join(genotypes, by = c("group1" = "alias")) %>% 
+      select(-group1) %>% 
+      mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*")))) %>% 
+      rename(group1 = genotype1, group2 = genotype) %>% 
+      filter(!is.na(p.signif)) %>% 
+      filter(met == plotmets[[meta]], tissue2 == plottissues[[tiss]])
+    
+    
+    plot_out[[tiss]][[meta]] <- fc %>% 
+      filter(met == plotmets[[meta]], tissue == plottissues[[tiss]]) %>% 
+      mutate(treatment = fct_relevel(treatment, c("0.4", "0.6", "0.8", "1"))) %>% 
+      ggplot(aes(x = genotype, y = mean_fc)) +
+      geom_col(position = "dodge", aes(fill = genotype), color = "black") +
+      geom_errorbar(aes(ymin = (mean_fc-se), ymax = (mean_fc + se)), position = position_dodge(0.9), width = 0.25, size = 0.75)+
+      stat_pvalue_manual(sig_bar, label = "p.signif", y.position = "y.position",
+                         step.increase = 0.07,
+                         hide.ns = T) +
+      facet_grid(cols = vars(treatment)) +
+      theme(axis.text.x = element_markdown(angle = 45, hjust = 1),
+            panel.background = element_rect(fill = "white"),
+            panel.border = element_rect(color = "black",fill = NA),
+            text = element_text(size = 14),
+            legend.text = element_markdown()) +
+      ylab("Mean fold-change")+
+      ggtitle(label = str_c(labelnames[[meta]], "in", plottissues[[tiss]], sep = " ")) +
+      scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+  }
+  
+}
+
+plot_out[[1]][[2]]
+
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_GC_col_plots.pdf"),
+    width = 9,
+    height = 6)
+
+for (tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    print(plot_out[[tiss]][[meta]])  
+  }
+  
+}
+
+dev.off()
+
+# Per metabolite comparisons unscaled----------------------------------------------
+
+norm_MM_1 <- means %>% 
+  filter(genotype == "MoneyMaker", treatment == 1) %>% 
+  select(tissue, met, MM_mean = mean)
+
+fc_1_ind <- means_tec_rep %>% 
+  left_join(norm_MM_1) %>% 
+  mutate(fc = mean_tec_rep/MM_mean)
+
+fc_1 <- means_tec_rep %>% 
+  left_join(norm_MM_1) %>% 
+  mutate(fc = mean_tec_rep/MM_mean) %>% 
+  group_by(tissue, treatment, alias, genotype, met) %>% 
+  summarise(mean_fc = mean(fc),
+            sd = sd(fc),
+            n = n()) %>% 
+  mutate(se = sd/sqrt(n),
+         group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  left_join(sig_GC, by = c("tissue", "treatment", "met", "alias" = "group2")) %>% 
+  group_by(met) %>% 
+  #filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>% 
+  left_join(met_dat) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  ungroup() %>% 
+  mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+         genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+
+plotmets <- fc_1 %>% distinct(met) %>% as_vector()
+plottissues <- fc_1 %>% distinct(tissue) %>% as_vector() %>% as.character
+
+labelnames <- plotmets %>% as_tibble() %>% 
+  rename(met = value) %>% 
+  left_join(met_dat) %>% left_join(GC_classes) %>%
+  select(Compound_Name) %>% as_vector()
+
+
+plot_out <- vector("list", length = length(plotmets))
+per_comp_y <- fc_1 %>%
+  group_by(tissue, treatment, met) %>% 
+  summarise(mean = max(mean_fc),
+            se = max(se))
+
+for(tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    
+    
+    sig_bar <- fc_1 %>% 
+      group_by(tissue, treatment, met) %>% 
+      mutate(tot_val = max(mean_fc + se)) %>%
+      mutate(y.position = tot_val + 0.25*max(tot_val)) %>%
+      ungroup() %>% 
+      rename(genotype1 = genotype) %>% 
+      left_join(genotypes, by = c("group1" = "alias")) %>% 
+      select(-group1) %>% 
+      mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*")))) %>% 
+      rename(group1 = genotype1, group2 = genotype) %>% 
+      filter(!is.na(p.signif)) %>% 
+      filter(met == plotmets[[meta]], tissue2 == plottissues[[tiss]])
+    
+    
+    plot_out[[tiss]][[meta]] <- fc_1 %>% 
+      filter(met == plotmets[[meta]], tissue == plottissues[[tiss]]) %>% 
+      mutate(treatment = fct_relevel(treatment, c("0.4", "0.6", "0.8", "1"))) %>% 
+      ggplot(aes(x = genotype, y = mean_fc)) +
+      geom_col(position = "dodge", aes(fill = genotype), color = "black") +
+      geom_errorbar(aes(ymin = (mean_fc-se), ymax = (mean_fc + se)), position = position_dodge(0.9), width = 0.25, size = 0.75)+
+      stat_pvalue_manual(sig_bar, label = "p.signif", y.position = "y.position",
+                         step.increase = 0.07,
+                         hide.ns = T) +
+      facet_grid(cols = vars(treatment)) +
+      theme(axis.text.x = element_markdown(angle = 45, hjust = 1),
+            panel.background = element_rect(fill = "white"),
+            panel.border = element_rect(color = "black",fill = NA),
+            text = element_text(size = 14),
+            legend.text = element_markdown()) +
+      ylab("Mean fold-change")+
+      ggtitle(label = str_c(labelnames[[meta]], "in", plottissues[[tiss]], sep = " ")) +
+      scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+  }
+  
+}
+
+plot_out[[tiss]][[meta]]
+plot_out[[1]][[2]]
+
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_scaled1_GC_col_plots.pdf"),
+    width = 15.8/2.54,
+    height = 8/2.54)
+
+for (tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    print(plot_out[[tiss]][[meta]])  
+  }
+  
+}
+
+dev.off()
+
+
+# Heatmap scaled Wildtype-----------------------------------------------------------------
+
+heat.GC <- fc %>% 
+  group_by(Compound_Name, met) %>% 
+  mutate(log_norm = log2(mean_fc),
+         log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+         log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>% 
+  ungroup() %>% 
+  #mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>% 
+  pivot_wider(id_cols = c(Compound_Name, met),
+              names_from = group,
+              values_from = log_norm) %>% 
+  left_join(GC_classes) %>%  
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>% 
+  select(contains("leaves"), contains("fruits")) %>% as.matrix()
+
+annotation_row <- heat.GC %>% 
+  select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- sam_dat %>% 
+  filter(exp == 1, class == "sample") %>%
+  distinct(tissue, treatment, genotype) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  as.data.frame()
+
+rownames(annotation_col) <- annotation_col$group
+
+annotation_col <- annotation_col %>% 
+  select(-group, treatment, tissue, genotype)
+
+heat.GC_signif <- means %>% 
+  distinct(tissue, treatment, alias, genotype, met) %>% 
+  #filter(genotype != "MoneyMaker") %>% 
+  left_join(sig_GC, by = c("met", "tissue", "treatment", "alias" = "group2")) %>% 
+  left_join(sig_mets) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")),
+         signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = signif) %>% 
+  left_join(met_dat) %>% 
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+heat_cols <- colnames(mat.heat.GC)
+
+rownames(heat.GC_signif) <- heat.GC_signif$met
+
+mat.heat.GC_signif <- heat.GC_signif %>% 
+  select(all_of(heat_cols)) %>% as.matrix()
+
+ann_colors = list(
+  tissue = c(fruits = "red",leaves = "darkgreen"),
+  genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+  treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+
+pheatmap.GC <- pheatmap(mat.heat.GC,
+                         #color = RColorBrewer::brewer.pal(7,"BuRd"),
+                         #cellwidth = 16,
+                         #cellheight = 4,
+                         #breaks = c(-2.5,-1.5, -0.5 ,0.5,1.5,2.5),
+                         cluster_rows = T,
+                         cluster_cols = T,
+                         annotation_names_row = F,
+                         show_rownames = F,
+                         annotation_row = annotation_row,
+                         annotation_col = annotation_col,
+                         display_numbers = mat.heat.GC_signif,
+                         number_color = "black",
+                         fontsize_number = 6,
+                        angle_col = 45,
+                        fontsize_col = 6,
+                         annotation_colors = ann_colors,
+                         filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+                                          "cmQTL_val_1_heatmap_rel_tissue_wt.jpg",
+                                          sep = "_")
+)
+
+
+
+# Levene analysis simple anova-------------------------------------------------------------
+
+miss_per_clust <- GC_long %>% 
+  group_by(met, tissue) %>% 
+  summarise(na = sum(is.na(area)),
+            n = n()) %>% 
+  ungroup() %>% 
+  mutate(percent_na_clust = na/n*100) %>% 
+  select(met, tissue, percent_na_clust)
+
+miss_any_treat <- miss_per_treat %>% 
+  group_by(met, tissue) %>% 
+  summarise(miss_comp_treat = if_else(any(percent_na == 100), T,F))
+
+means_tec_rep_lt <- means_tec_rep %>% 
+  group_by(met, tissue, treatment, genotype, alias) %>% 
+  mutate(lev_t = abs(log10(mean_tec_rep) - log10(median(mean_tec_rep)))) %>% 
+  ungroup()
+  
+
+means_lt <- means_tec_rep_lt %>% 
+  ungroup() %>% 
+  group_by(met, tissue, treatment, genotype, alias) %>% 
+  summarise(mean = mean(lev_t),
+            sd = sd(lev_t),
+            n = n()) %>% 
+  ungroup() %>% 
+  mutate(se = sd/sqrt(n))
+
+lt_tidy <- means_tec_rep_lt %>% 
+  filter(tissue == "fruits") %>% 
+  ungroup() %>% 
+  pivot_wider(id_cols = c(tissue, treatment, genotype, alias, LIMS_ID),
+                          names_from = met,
+              values_from = lev_t)
+
+lt_tidy_numeric <- lt_tidy %>% 
+  select(all_of(mets))
+
+lt_aov <- map(.x = lt_tidy_numeric, .f = ~aov(.x ~ alias*treatment, data = lt_tidy))
+
+lt_aov_tidy <- lt_aov %>% 
+  map(tidy) %>% 
+  map2(.y = names(lt_aov), .f = ~.x %>% mutate(lt_vars = .y)) %>% 
+  purrr::reduce(bind_rows)
+
+lt_tuk <- map(.x = lt_aov, .f = ~TukeyHSD(.x)) %>% 
+  map(.f = tidy) %>% 
+  map2(.y = names(lt_aov), .f = ~.x %>% mutate(var = .y)) %>% 
+  purrr::reduce(bind_rows)
+
+sig_GC_lt_groups <- lt_tuk %>% 
+  filter(term == "alias:treatment") %>% 
+  separate(col = contrast, into = c("group1", "group2"), sep = "-")
+
+sig_GC_lt <- sig_GC_lt_groups %>%  
+  separate(group1, into = c("alias1", "treatment1"), sep = ":") %>% 
+  separate(group2, into = c("alias2", "treatment2"), sep = ":") %>% 
+  filter(treatment1 == treatment2) %>% 
+  filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT") %>% 
+  ungroup() %>% 
+  mutate(p.signif = if_else(adj.p.value <= 0.05, "*", "ns")) %>% 
+  left_join(means_lt, by = c("var" = "met", "treatment1" = "treatment", "alias1" = "alias")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2,  p.signif, mean1 = mean, se1 = se) %>%
+  left_join(means_lt, by = c("var" = "met", "treatment2" = "treatment", "alias2" = "alias")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2,  p.signif, mean1, mean2 = mean, se1, se2 = se) %>% 
+  group_by(treatment1,  var) %>% 
+  mutate(tot_val1 = mean1 + se1,
+         tot_val2 = mean2 + se2,
+         y.position = 1.1*(max(tot_val1, tot_val2))) %>% 
+  rename(met = var,
+         group1 = alias1,
+         group2 = alias2,
+         treatment = treatment1) %>% 
+  mutate(met = as_factor(met),
+         group1 = as_factor(group1),
+         group2 = as_factor(group2),
+         treatment = as_factor(treatment)) %>% 
+  ungroup()
+
+sig_mets <- sig_GC_lt %>% 
+  filter(p.signif == "*") %>% 
+  distinct(met) %>% 
+  mutate(sig = T)
+
+# Levene with t-test ------------------------------------------------
+skip <-  F
+
+if(skip == T) {
+  print("Levene t-test skipped") 
+} else {
+GC_tidy <- means_tec_rep_lt %>% 
+  pivot_wider(id_cols = c( genotype, alias, treatment, tissue, LIMS_ID),
+              names_from = met,
+              values_from = lev_t) %>% 
+  mutate(group = str_c(tissue, treatment, alias, sep = "_"))
+
+GC_tidy_numeric <- GC_tidy %>% 
+  select(all_of(mets))
+
+GC_t <- map(.x = GC_tidy_numeric,
+            .f = ~pairwise.t.test(x = .x ,
+                                  g = GC_tidy$group,
+                                  p.adjust.method = "none",
+                                  pool.sd = F)) %>% 
+  map(.f = tidy)
+
+GC_t_tidy <- GC_t %>% 
+  map2(.y = names(GC_t), .f = ~.x %>% mutate(var = .y)) %>% 
+  purrr::reduce(bind_rows)
+
+sig_GC_lt<- GC_t_tidy %>% 
+  separate(group1, into = c("tissue1", "treatment1", "alias1"), sep = "_") %>% 
+  separate(group2, into = c("tissue2", "treatment2", "alias2"), sep = "_") %>% 
+  filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT",
+         tissue1 == tissue2, treatment1 == treatment2) %>% 
+  group_by(var) %>% 
+  mutate(adj.p.value = p.adjust(p.value)) %>% 
+  #mutate(adj.p.value = p.value * 121) %>% 
+  ungroup() %>% 
+  mutate(p.signif = if_else(p.value <= 0.05, "*", "ns")) %>% 
+  left_join(means_lt, by = c("var" = "met", "treatment1" = "treatment", "alias1" = "alias", "tissue1" = "tissue")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1 = mean, se1 = se, p.value, adj.p.value) %>%
+  left_join(means_lt, by = c("var" = "met", "treatment2" = "treatment", "alias2" = "alias", "tissue2" = "tissue")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1, mean2 = mean, se1, se2 = se, p.value, adj.p.value) %>% 
+  group_by(treatment1, tissue1, var) %>% 
+  mutate(tot_val1 = mean1 + se1,
+         tot_val2 = mean2 + se2,
+         y.position = 1.1*(max(tot_val1, tot_val2))) %>% 
+  ungroup() %>% 
+  rename(met = var,
+         group1 = alias1,
+         group2 = alias2,
+         treatment = treatment1,
+         tissue = tissue1) %>% 
+  left_join(miss_per_treat) %>% 
+  mutate(met = as_factor(met),
+         group1 = as_factor(group1),
+         group2 = as_factor(group2),
+         treatment = as_factor(treatment),
+         tissue = as_factor(tissue),
+         p.signif = if_else(percent_na >= 60, "ns", p.signif))
+
+
+sig_mets <- sig_GC_lt %>% 
+  filter(p.signif == "*") %>% 
+  distinct(met) %>% 
+  mutate(sig = T)
+
+}
+# Per metabolite comparisons unscaled levene----------------------------------------------
+
+norm_MM_1_lt <- means_lt %>% 
+  filter(genotype == "MoneyMaker", treatment == 1) %>% 
+  select(tissue, met, MM_mean = mean)
+
+fc_1_lt <- means_tec_rep_lt %>% 
+  #left_join(norm_MM_1) %>% 
+  #mutate(fc = mean_tec_rep/MM_mean) %>% 
+  group_by(tissue, treatment, alias, genotype, met) %>% 
+  summarise(mean_fc = mean(lev_t),
+            sd = sd(lev_t),
+            n = n()) %>% 
+  mutate(se = sd/sqrt(n),
+         group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  left_join(sig_GC_lt, by = c("tissue", "treatment", "met", "alias" = "group2")) %>% 
+  group_by(met) %>% 
+  #filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>% 
+  left_join(met_dat) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  ungroup() %>% 
+  mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+         genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+
+plotmets <- fc_1_lt %>% distinct(met) %>% as_vector()
+plottissues <- fc_1_lt %>% distinct(tissue) %>% as_vector() %>% as.character
+
+labelnames <- plotmets %>% as_tibble() %>% 
+  rename(met = value) %>% 
+  left_join(met_dat) %>% left_join(GC_classes) %>%
+  select(Compound_Name) %>% as_vector()
+
+
+plot_out <- vector("list", length = length(plotmets))
+per_comp_y <- fc_1_lt %>%
+  group_by(tissue, treatment, met) %>% 
+  summarise(mean = max(mean_fc),
+            se = max(se))
+
+for(tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    
+    
+    sig_bar <- fc_1_lt %>% 
+      group_by(tissue, treatment, met) %>% 
+      mutate(tot_val = max(mean_fc + se)) %>%
+      mutate(y.position = tot_val + 0.25*max(tot_val)) %>%
+      ungroup() %>% 
+      rename(genotype1 = genotype) %>% 
+      left_join(genotypes, by = c("group1" = "alias")) %>% 
+      select(-group1) %>% 
+      mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*")))) %>% 
+      rename(group1 = genotype1, group2 = genotype) %>% 
+      filter(!is.na(p.signif)) %>% 
+      filter(met == plotmets[[meta]], tissue2 == plottissues[[tiss]])
+    
+    
+    plot_out[[tiss]][[meta]] <- fc_1_lt %>% 
+      filter(met == plotmets[[meta]], tissue == plottissues[[tiss]]) %>% 
+      mutate(treatment = fct_relevel(treatment, c("0.4", "0.6", "0.8", "1"))) %>% 
+      ggplot(aes(x = genotype, y = mean_fc)) +
+      geom_col(position = "dodge", aes(fill = genotype), color = "black") +
+      geom_errorbar(aes(ymin = (mean_fc-se), ymax = (mean_fc + se)), position = position_dodge(0.9), width = 0.25, size = 0.75)+
+      stat_pvalue_manual(sig_bar, label = "p.signif", y.position = "y.position",
+                         step.increase = 0.07,
+                         hide.ns = T) +
+      facet_grid(cols = vars(treatment)) +
+      theme(axis.text.x = element_markdown(angle = 45, hjust = 1),
+            panel.background = element_rect(fill = "white"),
+            panel.border = element_rect(color = "black",fill = NA),
+            text = element_text(size = 14),
+            legend.text = element_markdown()) +
+      ylab("Mean fold-change")+
+      ggtitle(label = str_c(labelnames[[meta]], "in", plottissues[[tiss]], sep = " ")) +
+      scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+  }
+  
+}
+
+plot_out[[tiss]][[meta]]
+plot_out[[1]][[2]]
+
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_scaled1_GC_lev_t_col_plots.pdf"),
+    width = 15.8/2.54,
+    height = 8/2.54)
+
+for (tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    print(plot_out[[tiss]][[meta]])  
+  }
+  
+}
+
+dev.off()
+
+# CV individuals analysis -------------------------------------------------------------
+library(bootstrap)
+
+theta <- function(x){
+  sd(x)/mean(x)
+}
+
+cv <- means %>% 
+  ungroup() %>% 
+  group_by(genotype, tissue, alias, met) %>% 
+  summarise(grand_mean = mean(mean),
+            grand_sd = sd(mean),
+            n = n()) %>% 
+  mutate(cv = grand_sd/grand_mean)
+
+cv_jack_ind <- means_tec_rep %>% 
+  group_by(genotype, tissue, treatment, alias, met) %>% 
+  summarise(cv = jackknife(mean_tec_rep, theta)$jack.values) %>% 
+  mutate(jack_rep = row_number()) %>% 
+  ungroup()
+
+cv_jack_ind_mean <- cv_jack_ind %>% 
+  group_by(genotype,tissue, treatment, alias, met) %>% 
+  summarise(mean_cv = mean(cv),
+            sd_cv = sd(cv),
+            n = n()) %>% 
+  ungroup() %>% 
+  mutate(se_cv = sd_cv/sqrt(n))
+
+cv_jack_ind_wide <- cv_jack_ind %>% 
+  pivot_wider(id_cols = c(genotype, tissue, alias, treatment, jack_rep),
+              names_from = met,
+              values_from = cv) %>% 
+  mutate(group = str_c(tissue, treatment, alias, sep = "_"))
+
+cv_jack_ind_tidy_numeric <- cv_jack_ind_wide %>% 
+  select(all_of(mets))
+
+cv_jack_ind_t <- map(.x = cv_jack_ind_tidy_numeric,
+            .f = ~pairwise.t.test(x = .x ,
+                                  g = GC_tidy$group,
+                                  p.adjust.method = "none",
+                                  pool.sd = F)) %>% 
+  map(.f = tidy)
+
+cv_jack_ind_t_tidy <- cv_jack_ind_t %>% 
+  map2(.y = names(cv_jack_ind_t), .f = ~.x %>% mutate(var = .y)) %>% 
+  purrr::reduce(bind_rows)
+
+sig_cv_jack_ind<- cv_jack_ind_t_tidy %>% 
+  separate(group1, into = c("tissue1", "treatment1", "alias1"), sep = "_") %>% 
+  separate(group2, into = c("tissue2", "treatment2", "alias2"), sep = "_") %>% 
+  filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT",
+         tissue1 == tissue2, treatment1 == treatment2) %>% 
+  group_by(var) %>% 
+  mutate(adj.p.value = p.adjust(p.value)) %>% 
+  #mutate(adj.p.value = p.value * 121) %>% 
+  ungroup() %>% 
+  mutate(p.signif = if_else(adj.p.value <= 0.05, "*", "ns")) %>% 
+  left_join(means, by = c("var" = "met", "treatment1" = "treatment", "alias1" = "alias", "tissue1" = "tissue")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1 = mean, se1 = se, p.value, adj.p.value) %>%
+  left_join(means, by = c("var" = "met", "treatment2" = "treatment", "alias2" = "alias", "tissue2" = "tissue")) %>% 
+  select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1, mean2 = mean, se1, se2 = se, p.value, adj.p.value) %>% 
+  group_by(treatment1, tissue1, var) %>% 
+  mutate(tot_val1 = mean1 + se1,
+         tot_val2 = mean2 + se2,
+         y.position = 1.1*(max(tot_val1, tot_val2))) %>% 
+  ungroup() %>% 
+  rename(met = var,
+         group1 = alias1,
+         group2 = alias2,
+         treatment = treatment1,
+         tissue = tissue1) %>% 
+  left_join(miss_per_treat) %>% 
+  mutate(met = as_factor(met),
+         group1 = as_factor(group1),
+         group2 = as_factor(group2),
+         treatment = as_factor(treatment),
+         tissue = as_factor(tissue),
+         p.signif = if_else(percent_na >= 60, "ns", p.signif))
+
+
+sig_mets <- sig_cv_jack_ind %>% 
+  filter(p.signif == "*") %>% 
+  distinct(met) %>% 
+  mutate(sig = T)
+
+
+# Per metabolite comparisons cv ind unscaled----------------------------------------------
+
+norm_MM_1 <- cv_jack_ind %>% 
+  filter(genotype == "MoneyMaker", treatment == 1) %>% 
+  select(tissue, met, MM_mean = cv)
+
+fc_1_ind_cv_ind <- cv_jack_ind %>% 
+  left_join(norm_MM_1) %>% 
+  mutate(fc = cv/MM_mean)
+
+fc_1_cv_ind <- cv_jack_ind %>% 
+  left_join(norm_MM_1) %>% 
+  mutate(fc = cv/MM_mean) %>% 
+  group_by(tissue, treatment, alias, genotype, met) %>% 
+  summarise(mean_fc = mean(fc),
+            sd = sd(fc),
+            n = n()) %>% 
+  mutate(se = sd/sqrt(n),
+         group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  left_join(sig_cv_jack_ind, by = c("tissue", "treatment", "met", "alias" = "group2")) %>% 
+  group_by(met) %>% 
+  #filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>% 
+  left_join(met_dat) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  ungroup() %>% 
+  mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+         genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+
+plotmets <- fc_1_cv_ind %>% distinct(met) %>% as_vector()
+plottissues <- fc_1_cv_ind %>% distinct(tissue) %>% as_vector() %>% as.character
+
+labelnames <- plotmets %>% as_tibble() %>% 
+  rename(met = value) %>% 
+  left_join(met_dat) %>% left_join(GC_classes) %>%
+  select(Compound_Name) %>% as_vector()
+
+
+plot_out <- vector("list", length = length(plotmets))
+per_comp_y <- fc_1_cv_ind %>%
+  group_by(tissue, treatment, met) %>% 
+  summarise(mean = max(mean_fc),
+            se = max(se))
+
+for(tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    
+    
+    sig_bar <- fc_1_cv_ind %>% 
+      group_by(tissue, treatment, met) %>% 
+      mutate(tot_val = max(mean_fc + se)) %>%
+      mutate(y.position = tot_val + 0.25*max(tot_val)) %>%
+      ungroup() %>% 
+      rename(genotype1 = genotype) %>% 
+      left_join(genotypes, by = c("group1" = "alias")) %>% 
+      select(-group1) %>% 
+      mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*")))) %>% 
+      rename(group1 = genotype1, group2 = genotype) %>% 
+      filter(!is.na(p.signif)) %>% 
+      filter(met == plotmets[[meta]], tissue2 == plottissues[[tiss]])
+    
+    
+    plot_out[[tiss]][[meta]] <- fc_1_cv_ind %>% 
+      filter(met == plotmets[[meta]], tissue == plottissues[[tiss]]) %>% 
+      mutate(treatment = fct_relevel(treatment, c("0.4", "0.6", "0.8", "1"))) %>% 
+      ggplot(aes(x = genotype, y = mean_fc)) +
+      geom_col(position = "dodge", aes(fill = genotype), color = "black") +
+      geom_errorbar(aes(ymin = (mean_fc-se), ymax = (mean_fc + se)), position = position_dodge(0.9), width = 0.25, size = 0.75)+
+      stat_pvalue_manual(sig_bar, label = "p.signif", y.position = "y.position",
+                         step.increase = 0.07,
+                         hide.ns = T) +
+      facet_grid(cols = vars(treatment)) +
+      theme(axis.text.x = element_markdown(angle = 45, hjust = 1),
+            panel.background = element_rect(fill = "white"),
+            panel.border = element_rect(color = "black",fill = NA),
+            text = element_text(size = 14),
+            legend.text = element_markdown()) +
+      ylab("Mean fold-change")+
+      ggtitle(label = str_c(labelnames[[meta]], "in", plottissues[[tiss]], sep = " ")) +
+      scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+  }
+  
+}
+
+plot_out[[tiss]][[meta]]
+plot_out[[1]][[2]]
+
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_scaled1_cv_ind_GC_col_plots.pdf"),
+    width = 15.8/2.54,
+    height = 8/2.54)
+
+for (tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    print(plot_out[[tiss]][[meta]])  
+  }
+  
+}
+
+dev.off()
+
+
+# Heatmap scaled Wildtype-----------------------------------------------------------------
+
+norm_MM <- cv_jack_ind_mean %>% 
+  filter(genotype == "MoneyMaker") %>% 
+  select(tissue, treatment, met, MM_mean = mean_cv)
+
+fc_cv_ind <- cv_jack_ind %>% 
+  left_join(norm_MM) %>% 
+  mutate(fc = cv/MM_mean) %>% 
+  group_by(tissue, treatment, alias, genotype, met) %>% 
+  summarise(mean_fc = mean(fc),
+            sd = sd(fc),
+            n = n()) %>% 
+  mutate(se = sd/sqrt(n),
+         group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  left_join(sig_cv_jack_ind, by = c("tissue", "treatment", "met", "alias" = "group2")) %>% 
+  group_by(met) %>% 
+  filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>% 
+  left_join(met_dat) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  ungroup() %>% 
+  mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+         genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+
+heat.GC <- fc_cv_ind %>% 
+  group_by(Compound_Name, met) %>% 
+  mutate(log_norm = log2(mean_fc),
+         log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+         log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>% 
+  ungroup() %>% 
+  #mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>% 
+  pivot_wider(id_cols = c( Compound_Name, met),
+              names_from = group,
+              values_from = log_norm) %>% 
+  left_join(GC_classes) %>%  
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>% 
+  select(contains("leaves"), contains("fruits")) %>% as.matrix()
+
+annotation_row <- heat.GC %>% 
+  select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- sam_dat %>% 
+  filter(exp == 1, class == "sample") %>%
+  distinct(tissue, treatment, genotype) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%  
+  as.data.frame()
+
+rownames(annotation_col) <- annotation_col$group
+
+annotation_col <- annotation_col %>% 
+  select(-group, treatment, tissue, genotype)
+
+heat.GC_signif <- means %>% 
+  distinct(tissue, treatment, alias, genotype, met) %>% 
+  #filter(genotype != "MoneyMaker") %>% 
+  left_join(sig_cv_jack_ind, by = c("met", "tissue", "treatment", "alias" = "group2")) %>% 
+  left_join(sig_mets) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")),
+         signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = signif) %>% 
+  left_join(met_dat) %>% 
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+heat_cols <- colnames(mat.heat.GC)
+
+rownames(heat.GC_signif) <- heat.GC_signif$met
+
+mat.heat.GC_signif <- heat.GC_signif %>% 
+  select(all_of(heat_cols)) %>% as.matrix()
+
+ann_colors = list(
+  tissue = c(fruits = "red",leaves = "darkgreen"),
+  genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+  treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+
+pheatmap.GC <- pheatmap(mat.heat.GC,
+                        #color = colorRampPalette(c("blue","white", "red"))(9),
+                        #cellwidth = 16,
+                        #cellheight = 4,
+                        #breaks = seq(-2.25, 2.25, 0.5),
+                        cluster_rows = T,
+                        cluster_cols = T,
+                        annotation_names_row = F,
+                        show_rownames = F,
+                        annotation_row = annotation_row,
+                        annotation_col = annotation_col,
+                        display_numbers = mat.heat.GC_signif,
+                        number_color = "black",
+                        fontsize_number = 6,
+                        angle_col = 45,
+                        fontsize_col = 6,
+                        annotation_colors = ann_colors,
+                        filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+                                         "cmQTL_val_1_cv_ind_heatmap_rel_tissue_wt.jpg",
+                                         sep = "_")
+)
+
+
+
+
+
+# CV analysis -------------------------------------------------------------
+library(bootstrap)
+
+theta <- function(x){
+  sd(x)/mean(x)
+}
+
+cv <- means %>% 
+  ungroup() %>% 
+  group_by(genotype, tissue, alias, met) %>% 
+  summarise(grand_mean = mean(mean),
+            grand_sd = sd(mean),
+            n = n()) %>% 
+  mutate(cv = grand_sd/grand_mean)
+
+cv_jack <- means %>% 
+  group_by(genotype, tissue, alias, met) %>% 
+  summarise(cv = jackknife(mean, theta)$jack.values) %>% 
+  mutate(jack_rep = row_number()) %>% 
+  ungroup()
+
+cv_jack_mean <- cv_jack %>% 
+  group_by(genotype,tissue, alias, met) %>% 
+  summarise(mean_cv = mean(cv),
+            sd_cv = sd(cv),
+            n = n()) %>% 
+  ungroup() %>% 
+  mutate(se_cv = sd_cv/sqrt(n))
+
+cv_jack_wide <- cv_jack %>% 
+  pivot_wider(id_cols = c(genotype, tissue, alias, jack_rep),
+              names_from = met,
+              values_from = cv) %>% 
+  mutate(group = str_c(tissue, alias, sep = "_"))
+
+cv_jack_numeric <- cv_jack_wide %>% 
+  select(-c(genotype, tissue, alias, jack_rep, group))
+
+GC_cv_jack <- map(.x = cv_jack_numeric, .f = ~aov(.x ~ alias, data = cv_jack_wide))
+
+
+GC_jack_t <- map(.x = cv_jack_numeric, .f = ~pairwise.t.test(x = .x , g = cv_jack_wide$group,                                                              p.adjust.method = "none")) %>% 
+  map(.f = tidy)
+
+GC_jack_t_tidy <- GC_jack_t %>% 
+  map2(.y = names(GC_jack_t), .f = ~.x %>% mutate(var = .y)) %>% 
+  purrr::reduce(bind_rows)
+
+sig_GC_cv_groups <- GC_jack_t_tidy %>% 
+  separate(group1, into = c("tissue1", "alias1"), sep = "_") %>% 
+  separate(group2, into = c("tissue2", "alias2"), sep = "_") %>% 
+  filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT",
+         tissue1 == tissue2) %>% 
+  group_by(var) %>% 
+  mutate(adj.p.value = p.adjust(p.value)) %>% 
+  ungroup() %>% 
+  mutate(p.signif = if_else(adj.p.value <= 0.05, "*", "ns")) %>% 
+  left_join(cv_jack_mean, by = c("var" = "met", "alias1" = "alias", "tissue1" = "tissue")) %>% 
+  select(var,  alias1, alias2, p.signif, mean_cv1 = mean_cv, se_cv1 = se_cv, tissue1, tissue2, p.value, adj.p.value) %>%
+  left_join(cv_jack_mean, by = c("var" = "met", "alias2" = "alias", "tissue2" = "tissue")) %>% 
+  select(var, alias1, alias2, p.signif, mean_cv1, mean_cv2 = mean_cv, se_cv1, se_cv2 = se_cv, tissue1, tissue2, p.value, adj.p.value) %>% 
+  group_by(var, tissue1) %>% 
+  mutate(tot_val1 = mean_cv1 + se_cv1,
+         tot_val2 = mean_cv2 + se_cv2,
+         y.position = 1.1*(max(tot_val1, tot_val2))) %>%
+  ungroup() %>% 
+  left_join(genotypes, by = c("alias1" = "alias")) %>% 
+  rename(genotype1 = genotype) %>% 
+  left_join(genotypes, by = c("alias2" = "alias")) %>% 
+  rename(genotype2 = genotype) %>% 
+  rename(met = var,
+         group1 = genotype1,
+         group2 = genotype2,
+         tissue = tissue1) %>% 
+  left_join(miss_any_treat) %>% 
+  left_join(miss_per_clust) %>% 
+  group_by(met, tissue) %>% 
+  mutate(p.signif = if_else(miss_comp_treat == T| percent_na_clust >= 50, "ns", p.signif),
+         met = as_factor(met),
+         group1 = as_factor(group1),
+         group2 = as_factor(group2)) %>% 
+  ungroup()
+
+
+sig_mets_cv <- sig_GC_cv_groups %>% 
+  filter(p.signif == "*") %>% 
+  distinct(met, .keep_all = T) %>% 
+  mutate(sig = T) %>% 
+  select(met, sig)
+
+
+
+# CV Heatmap scaled per tissue-----------------------------------------------------------------
+
+heat.GC <- cv_jack_mean%>% 
+  group_by(tissue, met) %>% 
+  mutate(total_norm = mean_cv/mean(mean_cv),
+         log_norm = log2(total_norm),
+         log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+         log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>% 
+  ungroup() %>% 
+  left_join(sig_mets_cv) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_"))) %>% 
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = log_norm_level) %>% 
+  left_join(met_dat) %>%
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>% 
+  select(contains("leaves"), contains("fruits")) %>% as.matrix()
+
+annotation_row <- heat.GC %>% 
+  select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- sam_dat %>% 
+  filter(exp == 1, class == "sample") %>%
+  distinct(tissue, genotype) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_"))) %>%  
+  as.data.frame()
+
+rownames(annotation_col) <- annotation_col$group
+
+annotation_col <- annotation_col %>% 
+  select(-group, tissue, genotype)
+
+heat.GC_signif <- cv_jack_mean %>% 
+  distinct(tissue, alias, genotype, met) %>% 
+  #filter(genotype != "MoneyMaker") %>% 
+  left_join(sig_GC_cv_groups, by = c("met", "tissue" = "tissue2", "alias" = "alias2")) %>% 
+  left_join(sig_mets_cv) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_")),
+         signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = signif) %>% 
+  left_join(met_dat) %>% 
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+heat_cols <- colnames(mat.heat.GC)
+
+rownames(heat.GC_signif) <- heat.GC_signif$met
+
+mat.heat.GC_signif <- heat.GC_signif %>% 
+  select(all_of(heat_cols)) %>% as.matrix()
+
+ann_colors = list(
+  tissue = c(fruits = "red",leaves = "darkgreen"),
+  genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+  treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+
+
+pheatmap.GC <- pheatmap(mat.heat.GC,
+                         #color = plasma(14),
+                         #cellwidth = 16,
+                         #cellheight = 4,
+                         #breaks = c(-6.5,-5.5,-4.5,-3.5,-2.5,-1.5, -0.5 ,0.5,1.5,2.5,3.5,4.5,5.5,6.5),
+                         cluster_rows = T,
+                         cluster_cols = T,
+                         annotation_names_row = F,
+                         show_rownames = F,
+                         annotation_row = annotation_row,
+                         annotation_col = annotation_col,
+                         display_numbers = mat.heat.GC_signif,
+                         number_color = "black",
+                         fontsize_number = 6,
+                        angle_col = 45,
+                        fontsize_col = 6,
+                         annotation_colors = ann_colors,
+                         filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+                                          "cmQTL_val1_cv_heatmap_rel_tissue.jpg",
+                                          sep = "_")
+)
+
+
+
+# CV Heatmap unscaled-----------------------------------------------------------------
+
+heat.GC <- cv_jack_mean%>% 
+  group_by(tissue, met) %>% 
+  mutate(total_norm = mean_cv/mean(mean_cv),
+         log_norm = log2(total_norm),
+         log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+         log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>% 
+  ungroup() %>% 
+  left_join(sig_mets_cv) %>% 
+  #filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_"))) %>% 
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = mean_cv) %>% 
+  left_join(met_dat) %>%
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>% 
+  select(contains("leaves"), contains("fruits")) %>% as.matrix()
+
+annotation_row <- heat.GC %>% 
+  select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- sam_dat %>% 
+  filter(exp == 1, class == "sample") %>%
+  distinct(tissue, genotype) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_"))) %>%  
+  as.data.frame()
+
+rownames(annotation_col) <- annotation_col$group
+
+annotation_col <- annotation_col %>% 
+  select(-group, tissue, genotype)
+
+heat.GC_signif <- cv_jack_mean %>% 
+  distinct(tissue, alias, genotype, met) %>% 
+  #filter(genotype != "MoneyMaker") %>% 
+  left_join(sig_GC_cv_groups, by = c("met", "tissue" = "tissue2", "alias" = "alias2")) %>% 
+  left_join(sig_mets_cv) %>% 
+  #filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_")),
+         signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = signif) %>% 
+  left_join(met_dat) %>% 
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+heat_cols <- colnames(mat.heat.GC)
+
+rownames(heat.GC_signif) <- heat.GC_signif$met
+
+mat.heat.GC_signif <- heat.GC_signif %>% 
+  select(all_of(heat_cols)) %>% as.matrix()
+
+ann_colors = list(
+  tissue = c(fruits = "red",leaves = "darkgreen"),
+  genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+  treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+
+
+pheatmap.GC <- pheatmap(mat.heat.GC,
+                        #color = plasma(14),
+                        #cellwidth = 16,
+                        #cellheight = 4,
+                        #breaks = c(-6.5,-5.5,-4.5,-3.5,-2.5,-1.5, -0.5 ,0.5,1.5,2.5,3.5,4.5,5.5,6.5),
+                        cluster_rows = T,
+                        cluster_cols = T,
+                        annotation_names_row = F,
+                        show_rownames = F,
+                        annotation_row = annotation_row,
+                        annotation_col = annotation_col,
+                        #display_numbers = mat.heat.GC_signif,
+                        number_color = "black",
+                        fontsize_number = 6,
+                        angle_col = 45,
+                        fontsize_col = 6,
+                        annotation_colors = ann_colors,
+                        filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+                                         "cmQTL_val1_cv_heatmap_unscaled.jpg",
+                                         sep = "_")
+)
+
+
+dev.off()
+
+# CV dotplots -------------------------------------------------------------
+
+cv %>% 
+  filter(genotype == "MoneyMaker") %>% 
+  mutate(`cv > 1` = if_else(cv > 1, "yes", "no")) %>% 
+  ggplot(aes(x = tissue, y = cv)) +
+  geom_dotplot(aes(fill = `cv > 1`), stackdir = "center", binaxis = "y", binwidth = 0.1) +
+  geom_hline(aes(yintercept = 1), color = "red") +
+  theme(panel.background = element_rect(fill = "white"),
+        panel.border = element_rect(color = "black",fill = NA),
+        legend.position = "bottom",
+        legend.text = element_text(size = 8),
+        text = element_text(size = 10)) +
+  ylim(c(-0.15, 2))
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+             "cv_dotplot.jpg",
+             sep = "_"), units = "cm",
+       width = 15.9, height = 6)
+
+
+# CV Heatmap wildtype-----------------------------------------------------------------
+
+cv_jack_mean_wildtype <- cv_jack_mean %>% 
+  filter(genotype == "MoneyMaker") %>% 
+  select(tissue, met, mean_cv_wt = mean_cv)
+
+heat.GC <- cv_jack_mean%>% 
+  left_join(cv_jack_mean_wildtype) %>% 
+  mutate(fc_cv = mean_cv/mean_cv_wt) %>% 
+  group_by(tissue, met) %>% 
+  mutate(total_norm = mean_cv/mean(mean_cv),
+         log_norm = log2(fc_cv),
+         log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+         log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>% 
+  ungroup() %>% 
+  left_join(sig_mets_cv) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_"))) %>% 
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = log_norm) %>% 
+  left_join(met_dat) %>%
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>% 
+  select(contains("leaves"), contains("fruits")) %>% as.matrix()
+
+annotation_row <- heat.GC %>% 
+  select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- sam_dat %>% 
+  filter(exp == 1, class == "sample") %>%
+  distinct(tissue, genotype) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_"))) %>%  
+  as.data.frame()
+
+rownames(annotation_col) <- annotation_col$group
+
+annotation_col <- annotation_col %>% 
+  select(-group, tissue, genotype)
+
+heat.GC_signif <- cv_jack_mean %>% 
+  distinct(tissue, alias, genotype, met) %>% 
+  #filter(genotype != "MoneyMaker") %>% 
+  left_join(sig_GC_cv_groups, by = c("met", "tissue" = "tissue2", "alias" = "alias2")) %>% 
+  left_join(sig_mets_cv) %>% 
+  filter(sig == T) %>% 
+  mutate(group = as_factor(str_c(tissue, genotype, sep = "_")),
+         signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+  pivot_wider(id_cols = c(met),
+              names_from = group,
+              values_from = signif) %>% 
+  left_join(met_dat) %>% 
+  left_join(GC_classes) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  as.data.frame()
+
+heat_cols <- colnames(mat.heat.GC)
+
+rownames(heat.GC_signif) <- heat.GC_signif$met
+
+mat.heat.GC_signif <- heat.GC_signif %>% 
+  select(all_of(heat_cols)) %>% as.matrix()
+
+ann_colors = list(
+  tissue = c(fruits = "red",leaves = "darkgreen"),
+  genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+  treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+
+
+pheatmap.GC <- pheatmap(mat.heat.GC,
+                         #color = plasma(14),
+                         #cellwidth = 16,
+                         #cellheight = 4,
+                         #breaks = c(-6.5,-5.5,-4.5,-3.5,-2.5,-1.5, -0.5 ,0.5,1.5,2.5,3.5,4.5,5.5,6.5),
+                         cluster_rows = T,
+                         cluster_cols = T,
+                         annotation_names_row = F,
+                         show_rownames = F,
+                         annotation_row = annotation_row,
+                         annotation_col = annotation_col,
+                         display_numbers = mat.heat.GC_signif,
+                         number_color = "black",
+                         fontsize_number = 6,
+                        angle_col = 45,
+                        fontsize_col = 6,
+                         annotation_colors = ann_colors,
+                         filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+                                          "cmQTL_val1_cv_heatmap_rel_wt.jpg",
+                                          sep = "_")
+)
+
+# Per metabolite comparisons cv----------------------------------------------
+
+norm_MM <- means %>% 
+  filter(genotype == "MoneyMaker") %>% 
+  select(tissue, treatment, met, MM_mean = mean)
+
+fc_cv <- cv_jack%>% 
+  left_join(cv_jack_mean_wildtype) %>% 
+  mutate(fc_cv = cv/mean_cv_wt) %>% 
+  group_by(tissue, alias, genotype, met) %>% 
+  summarise(mean_fc = mean(fc_cv),
+            sd = sd(fc_cv),
+            n = n()) %>% 
+  mutate(se = sd/sqrt(n),
+         group = as_factor(str_c(tissue, genotype, sep = "_"))) %>%  
+  left_join(sig_GC_cv_groups, by = c("tissue", "met", "alias" = "alias2")) %>% 
+  group_by(met) %>% 
+  #filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>% 
+  left_join(met_dat) %>% 
+  arrange(Compound_Class, Compound_Name) %>% 
+  ungroup() %>% 
+  mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+         genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+
+plotmets <- fc_cv %>% distinct(met) %>% as_vector()
+plottissues <- fc_cv %>% distinct(tissue) %>% as_vector() %>% as.character
+
+labelnames <- plotmets %>% as_tibble() %>% 
+  rename(met = value) %>% 
+  left_join(met_dat) %>% left_join(GC_classes) %>%
+  select(Compound_Name) %>% as_vector()
+
+
+plot_out <- vector("list", length = length(plotmets))
+
+
+
+for(tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    
+    
+    sig_bar <- fc_cv %>% 
+      group_by(tissue, met) %>% 
+      mutate(tot_val = max(mean_fc + se)) %>%
+      mutate(y.position = tot_val + 0.25*max(tot_val)) %>%
+      ungroup() %>% 
+      filter(!is.na(p.signif)) %>% 
+      mutate(group2 = as_factor(if_else(group2 == "MoneyMaker", glue("{group2}"),glue("*{group2}*"))),
+             group2 = fct_relevel(group2, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*"))) %>% 
+      filter(met == plotmets[[meta]], tissue2 == plottissues[[tiss]])
+    
+    
+    plot_out[[tiss]][[meta]] <- fc_cv %>% 
+      filter(met == plotmets[[meta]], tissue == plottissues[[tiss]]) %>% 
+      ggplot(aes(x = genotype, y = mean_fc)) +
+      geom_col(position = "dodge", aes(fill = genotype), color = "black") +
+      geom_errorbar(aes(ymin = (mean_fc-se), ymax = (mean_fc + se)), position = position_dodge(0.9), width = 0.25, size = 0.75)+
+      stat_pvalue_manual(sig_bar, label = "p.signif", y.position = "y.position",
+                         step.increase = 0.07,
+                         hide.ns = T) +
+      theme(axis.text.x = element_markdown(angle = 45, hjust = 1),
+            panel.background = element_rect(fill = "white"),
+            panel.border = element_rect(color = "black",fill = NA),
+            text = element_text(size = 14),
+            legend.text = element_markdown()) +
+      ylab("Mean fold-change")+
+      ggtitle(label = str_c(labelnames[[meta]], "in", plottissues[[tiss]], sep = " ")) +
+      scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+  }
+  
+}
+
+plot_out[[1]][[2]]
+
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_cv_GC_col_plots.pdf"),
+    width = 15.8/2.54,
+    height = 8/2.54)
+
+for (tiss in seq_along(plottissues)) {
+  for (meta in seq_along(plotmets)) {
+    print(plot_out[[tiss]][[meta]])  
+  }
+  
+}
+
+dev.off()
+
+# Files to remake figures -------------------------------------------------
+
+write_csv(fc_1_ind, "individual_values.csv")
+write_csv(fc_1, "mean_values_se_n.csv")
+write_csv(sig_GC, "p_values.csv")
+
+write_csv(means_tec_rep_lt, "individual_values_levene.csv")
+write_csv(fc_1_lt, "mean_values_se_n_levene.csv")
+write_csv(sig_GC_lt, "p_values_levene.csv")
+
+write_csv(cv_jack_ind, "individual_values_cv.csv")
+write_csv(cv_jack_ind_mean, "mean_values_se_n_cv.csv")
+write_csv(sig_cv_jack_ind, "p_values_cv.csv")
+
+write_csv(cv_jack, "individual_values_cv.csv")
+write_csv(cv_jack_mean, "mean_values_se_n_cv.csv")
+write_csv(sig_GC_cv_groups, "p_values_cv.csv")
+
+# Log used code ------------------------------------------------------------
+
+file_name <- sys.frame(1)$ofile
+
+file.copy(file_name,
+          to = str_c(out_dir, str_remove(file_name, current), "_", str_replace_all(Sys.Date(),"^.{2}|-",""), ".R"),overwrite = T)