diff --git a/workflows/GC_MS_analysis/210903_primary_analysis_cmQTL_val_1_levene_RxNV.R b/workflows/GC_MS_analysis/210903_primary_analysis_cmQTL_val_1_levene_RxNV.R new file mode 100644 index 0000000000000000000000000000000000000000..291e54c653db34aefca7a7b2e4fafd9e7c2eba79 --- /dev/null +++ b/workflows/GC_MS_analysis/210903_primary_analysis_cmQTL_val_1_levene_RxNV.R @@ -0,0 +1,1841 @@ +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)