diff --git a/assays/cmQTL_val1_GH_2020_GC_MS/isa.assay.xlsx b/assays/cmQTL_val1_GH_2020_GC_MS/isa.assay.xlsx
index b8164b2d3765b2022354f54487473f2f9c124569..fec9ae80414f7a938612073d165363cbe9515ac4 100644
Binary files a/assays/cmQTL_val1_GH_2020_GC_MS/isa.assay.xlsx and b/assays/cmQTL_val1_GH_2020_GC_MS/isa.assay.xlsx differ
diff --git a/assays/cmQTL_val1_GH_2020_apolar_LC_MS/isa.assay.xlsx b/assays/cmQTL_val1_GH_2020_apolar_LC_MS/isa.assay.xlsx
index b66f56c12dd50fc8c4a220757376c0f1b7297ba5..d71ec871b3103f0e3903031114936aa8e1d2d5fb 100644
Binary files a/assays/cmQTL_val1_GH_2020_apolar_LC_MS/isa.assay.xlsx and b/assays/cmQTL_val1_GH_2020_apolar_LC_MS/isa.assay.xlsx differ
diff --git a/assays/cmQTL_val1_GH_2020_phenotype/README.md b/assays/cmQTL_val1_GH_2020_phenotype/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/assays/cmQTL_val1_GH_2020_phenotype/dataset/.gitkeep b/assays/cmQTL_val1_GH_2020_phenotype/dataset/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/assays/cmQTL_val1_GH_2020_phenotype/dataset/Phenotype_data.xlsx b/assays/cmQTL_val1_GH_2020_phenotype/dataset/Phenotype_data.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..49f0b056b692c0273bc06d1ae9649155755754d6
--- /dev/null
+++ b/assays/cmQTL_val1_GH_2020_phenotype/dataset/Phenotype_data.xlsx
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:560263fb9eaa085bd2f4bfca6a34994218faf8df0692c832630a26e398feded6
+size 26244
diff --git a/assays/cmQTL_val1_GH_2020_phenotype/isa.assay.xlsx b/assays/cmQTL_val1_GH_2020_phenotype/isa.assay.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..3e851c783b11d9dcd637fa384d777bfa5d58a342
Binary files /dev/null and b/assays/cmQTL_val1_GH_2020_phenotype/isa.assay.xlsx differ
diff --git a/assays/cmQTL_val1_GH_2020_phenotype/protocols/.gitkeep b/assays/cmQTL_val1_GH_2020_phenotype/protocols/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/assays/cmQTL_val1_GH_2020_polar_LC_MS/isa.assay.xlsx b/assays/cmQTL_val1_GH_2020_polar_LC_MS/isa.assay.xlsx
index 9b65eeee7b50883b797e546f4b7d8fe8554ac16b..13d3c2676ea954e9e12aae167a664e6330ef4120 100644
Binary files a/assays/cmQTL_val1_GH_2020_polar_LC_MS/isa.assay.xlsx and b/assays/cmQTL_val1_GH_2020_polar_LC_MS/isa.assay.xlsx differ
diff --git a/isa.investigation.xlsx b/isa.investigation.xlsx
index 8b6b9a0464054e9414d7c47c89590b25c330e900..0a206237fde46208957ac428eefd7c2a4218536d 100644
Binary files a/isa.investigation.xlsx and b/isa.investigation.xlsx differ
diff --git a/studies/cmQTL_val1_GH_2020_phenotype/README.md b/studies/cmQTL_val1_GH_2020_phenotype/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/studies/cmQTL_val1_GH_2020_phenotype/isa.study.xlsx b/studies/cmQTL_val1_GH_2020_phenotype/isa.study.xlsx
new file mode 100644
index 0000000000000000000000000000000000000000..dfb8b7f0bcac5823bcae6bbd2812959d40e6077f
Binary files /dev/null and b/studies/cmQTL_val1_GH_2020_phenotype/isa.study.xlsx differ
diff --git a/studies/cmQTL_val1_GH_2020_phenotype/protocols/.gitkeep b/studies/cmQTL_val1_GH_2020_phenotype/protocols/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/studies/cmQTL_val1_GH_2020_phenotype/resources/.gitkeep b/studies/cmQTL_val1_GH_2020_phenotype/resources/.gitkeep
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/workflows/GC_MS_correlation/220531_correlation_matrices.R b/workflows/GC_MS_correlation/220531_correlation_matrices.R
new file mode 100644
index 0000000000000000000000000000000000000000..1f44b872b59c794cd5ff9fa6c6309c9280b26bef
--- /dev/null
+++ b/workflows/GC_MS_correlation/220531_correlation_matrices.R
@@ -0,0 +1,281 @@
+library(tidyverse)
+library(glue)
+library(here)
+
+
+# Directory setting -------------------------------------------------------
+
+current <- getwd()
+source <- str_c(current,"/..")
+
+cur_date <- str_c(str_replace_all(Sys.Date(),"^.{2}|-",""))
+
+out <- str_c(cur_date, "correlation", sep = "_")
+
+if (file.exists(out)) {
+ cat("The folder already exists")
+} else {
+ dir.create(out)
+}
+
+out_dir <- str_c(current, out, sep = "/")
+
+latest <- str_sort(str_extract(list.files(pattern = "^\\d{6}_analysis$"),
+ pattern = "^\\d{6}_analysis"),
+ decreasing = T)[[1]]
+
+latest_analysis <- str_c(current, "/", latest)
+
+setwd(latest_analysis)
+
+latest_analysis_date <- str_extract(latest, pattern = "^\\d{6}")
+
+# File loading ------------------------------------------------------------
+
+fc_1 <- read_csv("mean_values_se_n.csv") %>%
+ mutate(genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+genotypes <- fc_1 %>%
+ distinct(alias, genotype) %>%
+ mutate(genotype_label = str_remove_all(genotype, "\\*"))
+
+fc_1_ind <- read_csv("individual_values.csv") %>%
+ select(-genotype) %>%
+ left_join(genotypes)
+
+sig_GC <- read_csv("p_values.csv")
+
+fc_1_lt <- read_csv("mean_values_se_n_levene.csv") %>%
+ mutate(genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+fc_1_ind_lt <- read_csv("individual_values_levene.csv") %>%
+ select(-genotype) %>%
+ left_join(genotypes)
+
+sig_GC_lt <- read_csv("p_values_levene.csv")
+
+fc_1_cv <- read_csv("mean_values_se_n_cv.csv") %>%
+ mutate(genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+fc_1_ind_cv <- read_csv("individual_values_cv.csv") %>%
+ select(-genotype) %>%
+ left_join(genotypes)
+
+sig_GC_cv <- read_csv("p_values_cv.csv") %>%
+ mutate(group1 = as_factor(if_else(group1 == "MoneyMaker", glue("{group1}"),glue("*{group1}*"))),
+ group1 = fct_relevel(group1, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")),
+ group2 = as_factor(if_else(group2 == "MoneyMaker", glue("{group2}"),glue("*{group2}*"))),
+ group2 = fct_relevel(group2, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+setwd(out_dir)
+
+
+# Data combination --------------------------------------------------------
+
+
+met_dat <- fc_1 %>%
+ distinct(met, Compound_Name, Compound_Class)
+
+
+per_comp_y <- fc_1 %>%
+ group_by(tissue, treatment, met) %>%
+ summarise(tot_val = max(mean_fc + se))
+
+cb_scale <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442",
+ "#0072B2", "#D55E00","#000000")
+bw_scale <- c("black", "black", "black", "black", "black", "black", "black")
+
+
+# Correlation network -------------------------------------------------------------
+library(jmuOutlier)
+
+g_vec <- genotypes$genotype %>%
+ sort() %>%
+ as.character()
+
+g_vec_lab <- genotypes %>%
+ arrange(genotype) %>%
+ select(genotype_label) %>%
+ as_vector()
+
+
+## Loop start --------------------------------------------------------------
+
+
+for(g in seq_along(genotypes$genotype)) {
+ print(str_c("starting with", g_vec_lab[[g]], sep = " "))
+
+cor_base <- fc_1_ind %>%
+ filter(#Compound_Class %in% c("Phospholipid", "Phosphatidylcholine"),
+ genotype %in% c(g_vec[[g]]),
+ tissue == "fruits")
+
+cor.GC <- cor_base %>%
+ group_by(met) %>%
+ mutate(log_norm = log2(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(genotype, treatment, LIMS_ID),
+ names_from = met,
+ values_from = fc) %>%
+ #arrange(Compound_Class, Compound_Name) %>%
+ as.data.frame()
+
+cor.GC_numeric <- cor.GC %>%
+ select(starts_with("m_"))
+
+corr_mat <- cor(cor.GC_numeric, method = "spearman")
+
+x <- vector(mode = "list", length = length(colnames(cor.GC_numeric)))
+names(x) <- colnames(cor.GC_numeric)
+
+start <- Sys.time()
+for (i in seq_along(colnames(cor.GC_numeric))) {
+ #names(x[[i]]) <- colnames(cor.GC_numeric)
+ for (n in seq_along(colnames(cor.GC_numeric))) {
+ if(i > n) {
+ x[[i]][[n]] <- NA
+ } else {
+ x[[i]][[n]] <- perm.cor.test(cor.GC_numeric[[i]], cor.GC_numeric[[n]],
+ method = "spearman",
+ num.sim = 20000
+ )$p.value
+ }
+ }
+ print(i)
+}
+
+stop <- Sys.time()
+stop-start
+
+for (i in seq_along(colnames(cor.GC_numeric))) {
+ names(x[[i]]) <- colnames(cor.GC_numeric)
+}
+
+corr_p <- x %>%
+ reduce(bind_rows)
+
+corr_p_long <- corr_p %>%
+ mutate(met1 = colnames(cor.GC_numeric)) %>%
+ pivot_longer(cols = starts_with("m_"),
+ values_to = "cor_p_val",
+ names_to = "met2") %>%
+ mutate(met1_num = as.numeric(str_extract(met1,"(?<=m_)\\d+")),
+ met2_num = as.numeric(str_extract(met2,"(?<=m_)\\d+"))) %>%
+ mutate(small = if_else(met1_num < met2_num, met1_num, met2_num),
+ big = if_else(met1_num > met2_num, met1_num, met2_num),
+ sorted = str_c(small, big, sep = "_")) %>%
+ filter(!is.na(cor_p_val)) %>%
+ mutate(p_adj = p.adjust(cor_p_val, method = "fdr")) %>%
+ filter(!met1 == met2)
+
+corr_mat_long <- corr_mat %>%
+ as_tibble(rownames = "met1") %>%
+ pivot_longer(cols = starts_with("m_"),
+ values_to = "cor_val",
+ names_to = "met2") %>%
+ mutate(met1_num = as.numeric(str_extract(met1,"(?<=m_)\\d+")),
+ met2_num = as.numeric(str_extract(met2,"(?<=m_)\\d+"))) %>%
+ mutate(small = if_else(met1_num < met2_num, met1_num, met2_num),
+ big = if_else(met1_num > met2_num, met1_num, met2_num),
+ sorted = str_c(small, big, sep = "_")) %>%
+ distinct(sorted, .keep_all = T) %>%
+ filter(!met1 == met2) %>%
+ mutate(pos_neg = as_factor(if_else(cor_val > 0, "pos", "neg")),
+ weight = cor_val)
+
+corr_mat_comb <- corr_mat_long %>%
+ left_join(corr_p_long)
+
+write_csv(corr_mat_comb,
+ str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "_corr_mat_",
+ g_vec_lab[[g]],
+ ".csv"))
+
+}
+
+## Loop end ----------------------------------------------------------------
+
+stop()
+
+library(viridis)
+library(igraph)
+library(ggnetwork)
+
+library(ggraph)
+
+classes_tidy <- met_dat %>%
+ mutate(`ChEBI simplified` = if_else(str_detect(Compound_Class,"carbohydrate"),"carbohydrate or derivative",
+ if_else(str_detect(Compound_Class,"amino_acid"), "amino acid or derivative",
+ if_else(str_detect(Compound_Class,"carboxylic_acid"), "carboxylic acid",
+ "other"))))
+
+annotation_col_corr <- classes_tidy %>%
+ select(`ChEBI simplified`) %>%
+ as.data.frame()
+
+rownames(annotation_col_corr) <- classes_tidy$met
+
+vir_scale <- plasma(4, begin = 0.2, end = 1)
+ann_colors_corr = list(
+ `ChEBI simplified`= c("amino acid or derivative"= vir_scale[1],
+ "carbohydrate or derivative" = vir_scale[2],
+ "carboxylic acid" = vir_scale[3],
+ "other" = vir_scale[4])
+)
+
+
+
+ann_colors_graph = tibble(`ChEBI simplified`= names(ann_colors_corr$`ChEBI simplified`),
+ color = ann_colors_corr$`ChEBI simplified`)
+
+node_cols <- tibble(met = colnames(cor.GC_numeric)) %>%
+ left_join(classes_tidy) %>%
+ left_join(ann_colors_graph)
+
+node_col_vec <- node_cols$color
+
+
+net_raw <- graph_from_data_frame(corr_mat_comb, vertices = node_cols, directed = F)
+
+
+set.seed(200)
+plot.igraph(net_raw, vertex.color = node_col_vec, vertex.label = "")
+
+set.seed(121)
+layout <- layout_with_graphopt(net_raw , charge = 0.2,
+ spring.length = 1,
+ spring.constant = 200,
+ mass = 50)
+
+ggraph(net_raw, layout = layout) +
+ geom_edge_arc0(aes(color = pos_neg), strength = 0.1) +
+ geom_node_point(aes(color = `ChEBI simplified`), size = 3)
+# geom_node_text(aes(label = name))
+
+
+ggnetwork(net_raw, layout = layout) %>%
+ ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
+ geom_edges(aes(color = pos_neg),curvature = 0.1) +
+ geom_nodes(aes(fill = `ChEBI simplified`), color = "black", shape = 21, size = 3) +
+ theme_blank() +
+ scale_fill_manual(values = vir_scale, aesthetics = "fill")
+
+#other methods
+
+ corr_mat_plot <- corr_mat_long %>%
+ filter(abs(cor_val) >= 0.75) %>%
+ select(met1, met2) %>%
+ as.matrix()
+
+ corr_mat_plot <- corr_mat
+ corr_mat_plot[abs(corr_mat_plot)<0.5] <- 0
+
+ net_raw <- graph_from_edgelist(corr_mat_plot, directed = F)
+ net_raw <- graph_from_adjacency_matrix(corr_mat_plot,
+ weighted=T,
+ mode="undirected",
+ diag=F)
diff --git a/workflows/GC_MS_figures/Figure maker_cmQTL1_paper.R b/workflows/GC_MS_figures/Figure maker_cmQTL1_paper.R
new file mode 100644
index 0000000000000000000000000000000000000000..4a3cd4c337d894f93511e8688116bd2f598bfda0
--- /dev/null
+++ b/workflows/GC_MS_figures/Figure maker_cmQTL1_paper.R
@@ -0,0 +1,2362 @@
+rm(list = ls())
+library(tidyverse)
+library(ggpubr)
+library(glue)
+library(ggtext)
+library(cowplot)
+library(ggbeeswarm)
+library(extrafont)
+library(viridis)
+
+# Directory setting -------------------------------------------------------
+
+
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+getwd()
+
+current <- getwd()
+source <- str_c(current,"/..")
+
+cur_date <- str_c(str_replace_all(Sys.Date(),"^.{2}|-",""))
+
+out <- str_c(cur_date, "Figures", sep = "_")
+
+if (file.exists(out)) {
+ cat("The folder already exists")
+} else {
+ dir.create(out)
+}
+
+out_dir <- str_c(current, out, sep = "/")
+
+latest <- str_sort(str_extract(list.files(pattern = "^\\d{6}_analysis$"),
+ pattern = "^\\d{6}_analysis"),
+ decreasing = T)[[1]]
+
+latest_analysis <- str_c(current, "/", latest)
+
+setwd(latest_analysis)
+
+latest_analysis_date <- str_extract(latest, pattern = "^\\d{6}")
+
+# File loading ------------------------------------------------------------
+
+fc_1 <- read_csv("mean_values_se_n.csv") %>%
+ mutate(genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+genotypes <- fc_1 %>%
+ distinct(alias, genotype) %>%
+ mutate(genotype_label = str_remove_all(genotype, "\\*"))
+
+fc_1_ind <- read_csv("individual_values.csv") %>%
+ select(-genotype) %>%
+ left_join(genotypes)
+
+sig_GC <- read_csv("p_values.csv")
+
+fc_1_lt <- read_csv("mean_values_se_n_levene.csv") %>%
+ mutate(genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+fc_1_ind_lt <- read_csv("individual_values_levene.csv") %>%
+ select(-genotype) %>%
+ left_join(genotypes)
+
+sig_GC_lt <- read_csv("p_values_levene.csv")
+
+fc_1_cv <- read_csv("mean_values_se_n_cv.csv") %>%
+ mutate(genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+fc_1_ind_cv <- read_csv("individual_values_cv.csv") %>%
+ select(-genotype) %>%
+ left_join(genotypes)
+
+sig_GC_cv <- read_csv("p_values_cv.csv") %>%
+ mutate(group1 = as_factor(if_else(group1 == "MoneyMaker", glue("{group1}"),glue("*{group1}*"))),
+ group1 = fct_relevel(group1, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")),
+ group2 = as_factor(if_else(group2 == "MoneyMaker", glue("{group2}"),glue("*{group2}*"))),
+ group2 = fct_relevel(group2, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+
+setwd(out_dir)
+
+
+# Data combination --------------------------------------------------------
+
+
+met_dat <- fc_1 %>%
+ distinct(met, Compound_Name, Compound_Class)
+
+classes_tidy <- met_dat %>%
+ mutate(`ChEBI simplified` = if_else(str_detect(Compound_Class,"carbohydrate"),"carbohydrate or derivative",
+ if_else(str_detect(Compound_Class,"amino_acid"), "amino acid or derivative",
+ if_else(str_detect(Compound_Class,"carboxylic_acid"), "carboxylic acid",
+ "other")))) %>%
+ mutate(ChEBI = if_else(str_detect(Compound_Class,"carbohydrate"),"CH",
+ if_else(str_detect(Compound_Class,"amino_acid"), "AA",
+ if_else(str_detect(Compound_Class,"carboxylic_acid"), "CA",
+ "other"))),
+ ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other")))
+
+genotypes <- fc_1 %>%
+ distinct(alias, genotype) %>%
+ mutate(genotype_label = str_remove_all(genotype, "\\*"))
+
+per_comp_y <- fc_1 %>%
+ group_by(tissue, treatment, met) %>%
+ summarise(tot_val = max(mean_fc + se))
+
+cb_scale <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442",
+ "#0072B2", "#D55E00", "#CC79A7","#000000")
+bw_scale <- c("black", "black", "black", "black", "black", "black", "black")
+vir_scale <- plasma(4, begin = 0.2, end = 1)
+
+
+# Used plots --------------------------------------------------------------
+
+com_theme <- theme(axis.text.x = element_markdown(angle = 45, hjust = 1, size = 6, family = "sans"),
+ axis.text.y = element_text(size = 6, family = "sans"),
+ axis.title.x = element_blank(),
+ axis.title.y = element_text(size = 6, family = "sans"),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ strip.text = element_text(size = 8, family = "sans", margin = margin(t = 1, r = 1, b = 1, l = 1 , unit = "pt")),
+ text = element_text(size = 6, family = "sans"),
+ legend.title = element_blank(),
+ legend.text = element_markdown(size = 6),
+ plot.margin = unit(c(1,1,1,1), "mm"),
+ legend.margin = margin(t = 0, r = 2, b = 0, l = 2 , unit = "mm"))
+
+make_box_dot_plot <- function(plot_met, plot_tissue, plot_label, plot_legend, plot_genotypes, plot_fill){
+
+ per_comp_y <- fc_1_ind %>%
+ filter(genotype %in% plot_genotypes) %>%
+ group_by(tissue, treatment, met) %>%
+ summarise(min_y = min(fc),
+ max_y = max(fc)) %>%
+ mutate(tot_val = 1.1*max_y) %>%
+ ungroup()
+
+ sig_bar <- sig_GC %>%
+ left_join(genotypes, by = c("group1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("group2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y) %>%
+ group_by(tissue, treatment, met) %>%
+ mutate(y.position = 1.1 * tot_val,
+ p.signif = if_else(p.value <= 0.05, "*", "ns")) %>%
+ ungroup() %>%
+ filter(met == plot_met, tissue == plot_tissue,
+ group1 %in% plot_genotypes, group2 %in% plot_genotypes)
+
+ binwidth <- sig_bar %>%
+ summarise(min_y_comp = min(min_y),
+ max_y_comp = max(max_y)) %>%
+ mutate(binwidth = (max_y_comp - min_y_comp)/50) %>%
+ select(binwidth) %>%
+ as_vector() %>%
+ as.numeric()
+
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ plot_out <- fc_1_ind %>%
+ filter(met == plot_met, tissue == plot_tissue, genotype %in% plot_genotypes) %>%
+ ggplot(aes(x = genotype, y = fc)) +
+ geom_boxplot(position = "dodge", aes(fill = genotype), color = "black") +
+ geom_dotplot(aes(fill = genotype, color = genotype), binaxis = "y", stackdir = "center", binwidth = binwidth, dotsize = 1.5) +
+ 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(legend.position = plot_legend) +
+ com_theme +
+ ylab("Mean fold-change")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ # ggtitle(label = "") +
+ scale_fill_manual(values = cb_scale[plot_fill], aesthetics = "fill") +
+ scale_color_manual(values = bw_scale, aesthetics = "color")
+
+ plot_out
+
+}
+
+make_box_dot_plot_levene <- function(plot_met, plot_tissue, plot_label, plot_legend, plot_genotypes, plot_fill){
+
+ per_comp_y <- fc_1_ind_lt %>%
+ filter(genotype %in% plot_genotypes) %>%
+ group_by(tissue, treatment, met) %>%
+ summarise(min_y = min(lev_t),
+ max_y = max(lev_t)) %>%
+ mutate(tot_val = 1.1*max_y) %>%
+ ungroup()
+
+ sig_bar <- sig_GC_lt %>%
+ left_join(genotypes, by = c("group1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("group2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y) %>%
+ group_by(tissue, treatment, met) %>%
+ mutate(y.position = 1.1 * tot_val,
+ p.signif = if_else(p.value <= 0.05, "*", "ns")) %>%
+ ungroup() %>%
+ filter(met == plot_met, tissue == plot_tissue,
+ group1 %in% plot_genotypes, group2 %in% plot_genotypes)
+
+ binwidth <- sig_bar %>%
+ summarise(min_y_comp = min(min_y),
+ max_y_comp = max(max_y)) %>%
+ mutate(binwidth = (max_y_comp - min_y_comp)/50) %>%
+ select(binwidth) %>%
+ as_vector() %>%
+ as.numeric()
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ plot_out <- fc_1_ind_lt %>%
+ filter(met == plot_met, tissue == plot_tissue, genotype %in% plot_genotypes) %>%
+ ggplot(aes(x = genotype, y = lev_t)) +
+ geom_boxplot(position = "dodge", aes(fill = genotype), color = "black") +
+ geom_dotplot(aes(fill = genotype, color = genotype), binaxis = "y", stackdir = "center", binwidth = binwidth, dotsize = 2) +
+ 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(legend.position = plot_legend) +
+ com_theme +
+ ylab("Levene's transformed value")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ #ggtitle(label = plot_label) +
+ scale_fill_manual(values = cb_scale[plot_fill], aesthetics = "fill") +
+ scale_color_manual(values = bw_scale, aesthetics = "color")
+
+ plot_out
+
+}
+
+
+make_box_dot_plot_cv <- function(plot_met, plot_tissue, plot_label, plot_legend, plot_genotypes, plot_fill){
+
+ per_comp_y_cv <- fc_1_ind_cv %>%
+ filter(genotype %in% plot_genotypes) %>%
+ group_by(tissue, met) %>%
+ summarise(min_y = min(cv),
+ max_y = max(cv)) %>%
+ mutate(tot_val = 1.1*max_y) %>%
+ ungroup()
+
+ sig_bar <- sig_GC_cv %>%
+ left_join(per_comp_y_cv) %>%
+ group_by(tissue, met) %>%
+ mutate(y.position = 1.1 * tot_val,
+ p.signif = if_else(p.value <= 0.05, "*", "ns")) %>%
+ ungroup() %>%
+ filter(met == plot_met, tissue == plot_tissue,
+ group1 %in% plot_genotypes, group2 %in% plot_genotypes)
+
+ binwidth <- sig_bar %>%
+ summarise(min_y_comp = min(min_y),
+ max_y_comp = max(max_y)) %>%
+ mutate(binwidth = (max_y_comp - min_y_comp)/50) %>%
+ select(binwidth) %>%
+ as_vector() %>%
+ as.numeric()
+
+ plot_out <- fc_1_ind_cv %>%
+ left_join(met_dat) %>%
+ filter(met == plot_met, tissue == plot_tissue, genotype %in% plot_genotypes) %>%
+ ggplot(aes(x = genotype, y = cv)) +
+ geom_boxplot(position = "dodge", aes(fill = genotype), color = "black") +
+ geom_dotplot(aes(fill = genotype, color = genotype), binaxis = "y", stackdir = "center", binwidth = binwidth, dotsize = 2) +
+ stat_pvalue_manual(sig_bar, label = "p.signif", y.position = "y.position",
+ step.increase = 0.07,
+ hide.ns = T) +
+ #facet_grid(cols = vars(treatment)) +
+ com_theme +
+ theme(legend.position = plot_legend) +
+ ylab("CV")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_manual(values = cb_scale[plot_fill], aesthetics = "fill") +
+ scale_color_manual(values = bw_scale, aesthetics = "color")
+
+ plot_out
+
+}
+
+
+# Malic acid fruits PanK4-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_70", "fruits", "Malic acid nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+
+saveRDS(last_plot(), "Malic_acid_nominal_pank4-1_fruits.RDS")
+
+p2 <- make_box_dot_plot_cv("m_70", "fruits", "Malic acid CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+saveRDS(last_plot(), "Malic_acid_cv_pank4-1_fruits.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_70", "fruits", "Malic acid", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+saveRDS(leg, "Malic_acid_legend_pank4-1_fruits.RDS")
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Malic_acid_fruits_pank4-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+p2 <- make_box_dot_plot_levene("m_70", "fruits", "Malic acid levene transformed", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+ggsave("Malic_acid_fruits_pank4-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+saveRDS(last_plot(), "Malic_acid_lt_pank4-1_fruits.RDS")
+
+# Malic acid leaves PanK4-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_70", "leaves", "Malic acid nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+
+saveRDS(last_plot(), "Malic_acid_nominal_pank4-1_leaves.RDS")
+
+p2 <- make_box_dot_plot_cv("m_70", "leaves", "Malic acid CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+
+saveRDS(last_plot(), "Malic_acid_cv_pank4-1_leaves.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_70", "leaves", "Malic acid", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+
+saveRDS(leg, "Malic_acid_legend_pank4-1_leaves.RDS")
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Malic_acid_leaves_panK4-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p2 <- make_box_dot_plot_levene("m_70", "leaves", "Malic acid", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+ggsave("Malic_acid_leaves-panK4-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+saveRDS(last_plot(), "Malic_acid_lt_pank4-1_leaves.RDS")
+
+
+# Malic acid fruits and leaves pank4-1 ------------------------------------
+
+p1 <- make_box_dot_plot("m_70", "fruits", "Malic acid nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+
+p2 <- make_box_dot_plot_cv("m_70", "fruits", "Malic acid CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+
+p3 <- make_box_dot_plot("m_70", "leaves", "Malic acid nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p3
+
+p4 <- make_box_dot_plot_cv("m_70", "leaves", "Malic acid CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p4
+
+leg <- get_legend(make_box_dot_plot("m_70", "leaves", "Malic acid", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Malic_acid_fruits_leaves_panK4-1_combined.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+
+# Phenylalanine fruits log2-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_81", "fruits", "Phenylalanine nominal", "none", c("MoneyMaker","*log2-1*"), c(1,3))
+p1
+
+saveRDS(last_plot(), "Phenylalanine_nominal_log2-1_fruits.RDS")
+
+p2 <- make_box_dot_plot_cv("m_81", "fruits", "Phenylalanine CV", "none", c("MoneyMaker","*log2-1*"), c(1,3))
+p2
+
+saveRDS(last_plot(), "Phenylalanine_CV_log2-1_fruits.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_81", "fruits", "Phenylalanine", "bottom", c("MoneyMaker","*log2-1*"), c(1,3)))
+leg
+
+saveRDS(leg, "Phenylalanine_leg_log2-1_fruits.RDS")
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Phenylalanine_fruits_log2-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p2 <- make_box_dot_plot_levene("m_81", "fruits", "Phenylalanine levene transformed", "bottom", c("MoneyMaker","*log2-1*"), c(1,3))
+p2
+ggsave("Phenylalanine_fruits_log2-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+saveRDS(last_plot(), "Phenylalanine_lt_log2-1_fruits.RDS")
+
+
+# Phenylalanine leaves log2-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_81", "leaves", "Phenylalanine nominal", "none", c("MoneyMaker","*log2-1*"), c(1,3))
+p1
+saveRDS(last_plot(), "Phenylalanine_nominal_log2-1_leaves.RDS")
+
+p2 <- make_box_dot_plot_cv("m_81", "leaves", "Phenylalanine CV", "none", c("MoneyMaker","*log2-1*"), c(1,3))
+p2
+saveRDS(last_plot(), "Phenylalanine_cv_log2-1_leaves.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_81", "leaves", "Phenylalanine", "bottom", c("MoneyMaker","*log2-1*"), c(1,3)))
+leg
+saveRDS(leg, "Phenylalanine_leg_log2-1_leaves.RDS")
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Phenylalanine_leaves_log2-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p2 <- make_box_dot_plot_levene("m_81", "leaves", "Phenylalanine levene transformed", "bottom", c("MoneyMaker","*log2-1*"), c(1,3))
+p2
+ggsave("Phenylalanine_leaves_log-2-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+saveRDS(last_plot(), "Phenylalanine_lt_log2-1_leaves.RDS")
+# Phenylalanine fruits and leaves log2-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_81", "fruits", "Phenylalanine nominal", "none", c("MoneyMaker","*log2-1*"), c(1,3))
+p1
+
+p2 <- make_box_dot_plot_cv("m_81", "fruits", "Phenylalanine CV", "none", c("MoneyMaker","*log2-1*"), c(1,3))
+p2
+
+p3 <- make_box_dot_plot("m_81", "leaves", "Phenylalanine nominal", "none", c("MoneyMaker","*log2-1*"), c(1,3))
+p3
+
+p4 <- make_box_dot_plot_cv("m_81", "leaves", "Phenylalanine CV", "none", c("MoneyMaker","*log2-1*"), c(1,3))
+p4
+
+leg <- get_legend(make_box_dot_plot("m_81", "leaves", "Phenylalanine", "bottom", c("MoneyMaker","*log2-1*"), c(1,3)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Phenylalanine_fruits_leaves_log2-1_combined.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+# Erythritol fruits PanK4-1-------------------------------------------------------
+
+
+p1 <- make_box_dot_plot("m_24", "fruits", "Erythritol nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+saveRDS(last_plot(), "Erythritol_nominal_pank4-1_fruits.RDS")
+
+p2 <- make_box_dot_plot_cv("m_24", "fruits", "Erythritol CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+saveRDS(last_plot(), "Erythritol_cv_pank4-1_fruits.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_24", "fruits", "Erythritol", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+saveRDS(leg, "Erythritol_leg_pank4-1_fruits.RDS")
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Erythritol_fruits_panK4-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+p2 <- make_box_dot_plot_levene("m_24", "fruits", "Erythritol levene transformed", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+ggsave("Erythritol_fruits_panK4-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+saveRDS(last_plot(), "Erythritol_lt_pank4-1_fruits.RDS")
+
+# Erythritol leaves PanK4-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_24", "leaves", "Erythritol nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+saveRDS(last_plot(), "Erythritol_nominal_pank4-1_leaves.RDS")
+
+p2 <- make_box_dot_plot_cv("m_24", "leaves", "Erythritol CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+saveRDS(last_plot(), "Erythritol_cv_pank4-1_leaves.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_24", "leaves", "Erythritol", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+saveRDS(leg, "Erythritol_leg_pank4-1_leaves.RDS")
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Erythritol_leaves_panK4-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p2 <- make_box_dot_plot_levene("m_24", "leaves", "Erythritol levene transformed", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+ggsave("Erythritol_leaves_panK4-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+saveRDS(last_plot(), "Erythritol_lt_pank4-1_leaves.RDS")
+
+
+# Isocitrate fruits PanK4-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_63", "fruits", "Isocitrate nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+saveRDS(last_plot(), "Isocit_nominal_pank4-1_fruits.RDS")
+
+p2 <- make_box_dot_plot_cv("m_63", "fruits", "Isocitrate CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+saveRDS(last_plot(), "Isocit_cv_pank4-1_fruits.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_63", "fruits", "Isocitrate", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+saveRDS(leg, "Isocit_leg_pank4-1_fruits.RDS")
+
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Isocitrate_fruits_panK4-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+p2 <- make_box_dot_plot_levene("m_63", "fruits", "Isocitrate levene transformed", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+ggsave("Isocitrate_fruits_panK4-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+saveRDS(last_plot(), "Isocit_lt_pank4-1_fruits.RDS")
+
+# Isocitrate leaves PanK4-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_63", "leaves", "Isocitrate nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+
+p2 <- make_box_dot_plot_cv("m_63", "leaves", "Isocitrate CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_63", "leaves", "Isocitrate", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Isocitrate_leaves_panK4-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p3 <- make_box_dot_plot_levene("m_63", "leaves", "Isocitrate levene transformed", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+ggsave("Isocitrate_leaves_panK4-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+saveRDS(p1, "Isocit_nominal_pank4-1_leaves.RDS")
+saveRDS(p2, "Isocit_cv_pank4-1_leaves.RDS")
+saveRDS(p3, "Isocit_lt_pank4-1_leaves.RDS")
+saveRDS(leg, "Isocit_leg_pank4-1_leaves.RDS")
+
+# Isocitrate fruits and leaves pank4-1 ------------------------------------
+
+p1 <- make_box_dot_plot("m_70", "fruits", "Isocitrate nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+
+p2 <- make_box_dot_plot_cv("m_70", "fruits", "Isocitrate CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+
+p3 <- make_box_dot_plot("m_70", "leaves", "Isocitrate nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p3
+
+p4 <- make_box_dot_plot_cv("m_70", "leaves", "Isocitrate CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p4
+
+leg <- get_legend(make_box_dot_plot("m_70", "leaves", "Isocitrate", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Isocitrate_fruits_leaves_panK4-1_combined.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+
+
+# GABA CV -----------------------------------------------------------------
+
+p4 <- make_box_dot_plot_cv("m_17", "fruits", "GABA", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p4
+
+saveRDS(last_plot(), "GABA_cv_pank4-1_fruits.RDS")
+
+# beta-alanine CV -----------------------------------------------------------------
+
+p4 <- make_box_dot_plot_cv("m_8", "fruits", "beta-alanine", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p4
+
+saveRDS(last_plot(), "beta-alanine_cv_pank4-1_fruits.RDS")
+
+# Valine CV -----------------------------------------------------------------
+
+p4 <- make_box_dot_plot_cv("m_122", "fruits", "Valine", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p4
+
+saveRDS(last_plot(), "Valine_cv_pank4-1_fruits.RDS")
+
+# Succinic acid fruits PanK4-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_102", "fruits", "Succinic acid nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+
+
+p2 <- make_box_dot_plot_cv("m_102", "fruits", "Succinic acid CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_102", "fruits", "Succinic acid", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Succinic acid_fruits_panK4-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+p3 <- make_box_dot_plot_levene("m_102", "fruits", "Succinic acid levene transformed", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+ggsave("Succinic acid_fruits_panK4-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+saveRDS(p1, "Succinic acid_nominal_pank4-1_fruits.RDS")
+saveRDS(p2, "Succinic acid_cv_pank4-1_fruits.RDS")
+saveRDS(p3, "Succinic acid_lt_pank4-1_fruits.RDS")
+saveRDS(leg, "Succinic acid_leg_pank4-1_fruits.RDS")
+
+# Succinic acid leaves PanK4-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_102", "leaves", "Succinic acid nominal", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p1
+
+p2 <- make_box_dot_plot_cv("m_102", "leaves", "Succinic acid CV", "none", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_102", "leaves", "Succinic acid", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Succinic acid_leaves_panK4-1_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p3 <- make_box_dot_plot_levene("m_102", "leaves", "Succinic acid levene transformed", "bottom", c("MoneyMaker","*panK4-1*"), c(1,2))
+p2
+ggsave("Succinic acid_leaves_panK4-1_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+saveRDS(p1, "Succinic acid_nominal_pank4-1_leaves.RDS")
+saveRDS(p2, "Succinic acid_cv_pank4-1_leaves.RDS")
+saveRDS(p3, "Succinic acid_lt_pank4-1_leaves.RDS")
+saveRDS(leg, "Succinic acid_leg_pank4-1_leaves.RDS")
+
+# Glucose fruits transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_44", "fruits", "Glucose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_44", "fruits", "Glucose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_44", "fruits", "Glucose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Glucose_fruits_combined_transp1-1.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+p3 <- make_box_dot_plot_levene("m_44", "fruits", "Glucose levene transformed", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+ggsave("Glucose_fruits_box_dot_plot_levene_transp1-1.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+saveRDS(p1, "Glucose_nominal_transp1-1_fruits.RDS")
+saveRDS(p2, "Glucose_cv_transp1-1_fruits.RDS")
+saveRDS(p3, "Glucose_lt_transp1-1_fruits.RDS")
+saveRDS(leg, "Glucose_leg_transp1-1_fruits.RDS")
+# Glucose leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_44", "leaves", "Glucose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_44", "leaves", "Glucose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_44", "leaves", "Glucose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Glucose_leaves_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p3 <- make_box_dot_plot_levene("m_44", "leaves", "Glucose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+ggsave("Glucose_leaves_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+saveRDS(p1, "Glucose_nominal_transp1-1_leaves.RDS")
+saveRDS(p2, "Glucose_cv_transp1-1_leaves.RDS")
+saveRDS(p3, "Glucose_lt_transp1-1_leaves.RDS")
+saveRDS(leg, "Glucose_leg_transp1-1_leaves.RDS")
+
+
+# Glucose fruits leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_44", "leaves", "Glucose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_44", "leaves", "Glucose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+p3 <- make_box_dot_plot("m_44", "fruits", "Glucose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p3
+
+p4 <- make_box_dot_plot_cv("m_44", "fruits", "Glucose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p4
+
+leg <- get_legend(make_box_dot_plot("m_44", "fruits", "Glucose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Glucose_fruits_leaves_combined_transp1-1.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+# Maltose fruits leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_72", "leaves", "Maltose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_72", "leaves", "Maltose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+p3 <- make_box_dot_plot("m_72", "fruits", "Maltose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p3
+
+p4 <- make_box_dot_plot_cv("m_72", "fruits", "Maltose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p4
+saveRDS(last_plot(), "Maltose_cv_transp1_fruits.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_72", "fruits", "Maltose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Maltose_fruits_leaves_combined_transp1-1.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+
+# Fructose-6-phosphate fruits leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_27", "leaves", "Fructose-6-phosphate nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_27", "leaves", "Fructose-6-phosphate CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+p3 <- make_box_dot_plot("m_27", "fruits", "Fructose-6-phosphate nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p3
+
+p4 <- make_box_dot_plot_cv("m_27", "fruits", "Fructose-6-phosphate CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p4
+
+saveRDS(last_plot(), "Fructose_6_phosphate_cv_transp1_fruits.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_27", "fruits", "Fructose-6-phosphate", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Fructose-6-phosphate_fruits_leaves_combined_transp1-1.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+p3 <- make_box_dot_plot_levene("m_72", "fruits", "Fructose-6-phosphate levene transformed", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+ggsave("Fructose_fruits_box_dot_plot_levene_transp1-1.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+
+# Glucose-6-phosphate fruits leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_41", "leaves", "Glucose-6-phosphate nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_41", "leaves", "Glucose-6-phosphate CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+p3 <- make_box_dot_plot("m_41", "fruits", "Glucose-6-phosphate nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p3
+
+p4 <- make_box_dot_plot_cv("m_41", "fruits", "Glucose-6-phosphate CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p4
+
+saveRDS(last_plot(), "Glucose_6_phosphate_cv_transp1_fruits.RDS")
+
+leg <- get_legend(make_box_dot_plot("m_41", "fruits", "Glucose-6-phosphate", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Glucose-6-phosphate_fruits_leaves_combined_transp1-1.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+
+
+# Fructose fruits transp1-1-------------------------------------------------------
+
+
+p1 <- make_box_dot_plot("m_29", "fruits", "Fructose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+
+p2 <- make_box_dot_plot_cv("m_29", "fruits", "Fructose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_29", "fruits", "Fructose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Fructose_fruits_combined_transp1-1.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+p3 <- make_box_dot_plot_levene("m_29", "fruits", "Fructose levene transformed", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+ggsave("Fructose_fruits_box_dot_plot_levene_transp1-1.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+saveRDS(p1, "Fructose_nominal_pank4-1_fruits.RDS")
+saveRDS(p2, "Fructose_cv_pank4-1_fruits.RDS")
+saveRDS(p3, "Fructose_lt_pank4-1_fruits.RDS")
+saveRDS(leg, "Fructose_leg_pank4-1_fruits.RDS")
+# Fructose leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_29", "leaves", "Fructose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_29", "leaves", "Fructose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_29", "leaves", "Fructose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Fructose_leaves_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p3 <- make_box_dot_plot_levene("m_29", "leaves", "Fructose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+ggsave("Fructose_leaves_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+saveRDS(p1, "Fructose_nominal_pank4-1_leaves.RDS")
+saveRDS(p2, "Fructose_cv_pank4-1_leaves.RDS")
+saveRDS(p3, "Fructose_lt_pank4-1_leaves.RDS")
+saveRDS(leg, "Fructose_leg_pank4-1_leaves.RDS")
+# Fructose fruits leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_29", "leaves", "Fructose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_29", "leaves", "Fructose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+p3 <- make_box_dot_plot("m_29", "fruits", "Fructose nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p3
+
+p4 <- make_box_dot_plot_cv("m_29", "fruits", "Fructose CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p4
+
+leg <- get_legend(make_box_dot_plot("m_29", "fruits", "Fructose", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Fructose_fruits_leaves_combined_transp1-1.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+
+# Tyrosine fruits transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_112", "fruits", "Tyrosine nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+
+p2 <- make_box_dot_plot_cv("m_112", "fruits", "Tyrosine CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_112", "fruits", "Tyrosine", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Tyrosine_fruits_combined_transp1-1.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+p3 <- make_box_dot_plot_levene("m_112", "fruits", "Tyrosine levene transformed", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+ggsave("Tyrosine_fruits_box_dot_plot_levene_transp1-1.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+saveRDS(p1, "Tyrosine_nominal_pank4-1_fruits.RDS")
+saveRDS(p2, "Tyrosine_cv_pank4-1_fruits.RDS")
+saveRDS(p3, "Tyrosine_lt_pank4-1_fruits.RDS")
+saveRDS(leg, "Tyrosine_leg_pank4-1_fruits.RDS")
+# Tyrosine leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_112", "leaves", "Tyrosine nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_112", "leaves", "Tyrosine CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+leg <- get_legend(make_box_dot_plot("m_112", "leaves", "Tyrosine", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Tyrosine_leaves_combined.png",width = 16.5, height = 10, units = "cm", dpi = 300)
+
+
+p3 <- make_box_dot_plot_levene("m_112", "leaves", "Tyrosine", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+ggsave("Tyrosine_leaves_box_dot_plot_levene.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+saveRDS(p1, "Tyrosine_nominal_pank4-1_leaves.RDS")
+saveRDS(p2, "Tyrosine_cv_pank4-1_leaves.RDS")
+saveRDS(p3, "Tyrosine_lt_pank4-1_leaves.RDS")
+saveRDS(leg, "Tyrosine_leg_pank4-1_leaves.RDS")
+
+# Tyrosine fruits leaves transp1-1-------------------------------------------------------
+
+p1 <- make_box_dot_plot("m_112", "leaves", "Tyrosine nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p1
+
+p2 <- make_box_dot_plot_cv("m_112", "leaves", "Tyrosine CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p2
+
+p3 <- make_box_dot_plot("m_112", "fruits", "Tyrosine nominal", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p3
+
+p4 <- make_box_dot_plot_cv("m_112", "fruits", "Tyrosine CV", "none", c("MoneyMaker","*transp1-1*"), c(1,4))
+p4
+
+leg <- get_legend(make_box_dot_plot("m_112", "fruits", "Tyrosine", "bottom", c("MoneyMaker","*transp1-1*"), c(1,4)))
+leg
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Tyrosine_fruits_leaves_combined_transp1-1.png",width = 16.5, height = 20, units = "cm", dpi = 300)
+
+
+# Clustering k-means---------------------------------------------------------
+
+library(pheatmap)
+
+heat_base <- fc_1 %>%
+ filter(#Compound_Class %in% c("Phospholipid", "Phosphatidylcholine"),
+ genotype %in% c("MoneyMaker"),
+ tissue == "fruits")
+
+heat.GC <- heat_base %>%
+ 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)),
+ z_score = (mean_fc - mean(mean_fc))/(max(mean_fc)-min(mean_fc))) %>%
+ ungroup() %>%
+ mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+ pivot_wider(id_cols = c(Compound_Name, Compound_Class, met),
+ names_from = group,
+ values_from = z_score) %>%
+ arrange(Compound_Class, Compound_Name) %>%
+ as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>%
+ select(contains("0.4"), contains("0.6"),
+ contains("0.8"),
+ contains("1"),
+ contains("leaves"), contains("fruits")) %>% as.matrix()
+
+set.seed(484)
+k_means <- kmeans(mat.heat.GC, centers = 6)
+
+clusters <- k_means$cluster
+
+clust_tidy <- tibble("clust" = clusters) %>%
+ mutate(met = names(clusters))
+
+heat_base %>%
+ left_join(clust_tidy) %>%
+ left_join(classes_tidy) %>%
+ 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)),
+ z_score = (mean_fc - mean(mean_fc))/(max(mean_fc)-min(mean_fc))) %>%
+ ungroup() %>%
+ group_by(clust, treatment) %>%
+ mutate(n = n(),
+ ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other"))) %>%
+ ungroup() %>%
+ ggplot(aes(x = treatment, y = z_score, color = ChEBI)) +
+ geom_line(aes(group = met)) +
+ geom_label(aes(x = 0.7, y = -0.6, label = n), inherit.aes = F) +
+ facet_wrap(vars(clust), nrow = 2) +
+ scale_color_manual(values = cb_scale[5:8]) +
+ com_theme +
+ ylab("range-scaled level") +
+ theme(legend.position = "bottom",
+ strip.text = element_blank())
+
+ggsave("clusters_MoneyMaker.pdf",
+ width = 8.25,
+ height = 8.25,
+ units = "cm",
+ dpi = 300)
+
+saveRDS(last_plot(), "clusters_MoneyMaker.RDS")
+
+annotation_row <- heat.GC %>%
+ select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- heat_base %>%
+ 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 <- heat_base %>%
+ 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.value >= 0.05| is.na(p.value), "","X")) %>%
+ pivot_wider(id_cols = c(met),
+ names_from = group,
+ values_from = signif) %>%
+ left_join(met_dat) %>%
+ #left_join(GCid_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 = cb_scale[1], `*panK4-1*` = "brown", `log2-1` = "blue", `*transp1-1*` = cb_scale[4]),
+ 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"))(65),
+ #cellwidth = 16,
+ #cellheight = 4,
+ breaks = seq(-3.25, 3.25, 0.1),
+ clustering_distance_rows = dist((mat.heat.GC), method = "euclidean"),
+ cluster_rows = T,
+ cluster_cols = F,
+ 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,
+ #filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ # "cmQTL_val_1_heatmap_rel_tissue_wt.jpg",
+ # sep = "_")
+)
+
+dev.off()
+
+
+# Correlation network -------------------------------------------------------------
+
+cor_theme <- theme(axis.text = element_blank(),
+ axis.title = element_blank(),
+ panel.background = element_rect(fill = "white"),
+ strip.text = element_text(size = 10, family = "sans"),
+ text = element_text(size = 6, family = "sans"),
+ legend.position = "",
+ legend.title = element_text(size = 8, family = "sans"),
+ legend.text = element_text(size = 6, family = "sans"),
+ plot.margin = unit(c(1,1,1,1), "pt"),
+ legend.box.margin = margin(t = 0, r = 0, b = 0, l = 0 , unit = "mm"),
+ legend.margin = margin(t = 0, r = 0, b = 0, l = 0 , unit = "mm"))
+
+library(viridis)
+library(igraph)
+library(ggnetwork)
+library(jmuOutlier)
+library(ggraph)
+
+#correlation matrices are preformed
+
+## Loop start --------------------------------------------------------------
+
+cor_base <- fc_1_ind %>%
+ filter(#Compound_Class %in% c("Phospholipid", "Phosphatidylcholine"),
+ genotype %in% c("MoneyMaker"),
+ tissue == "fruits")
+
+cor.GC <- cor_base %>%
+ group_by(met) %>%
+ mutate(log_norm = log2(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(genotype, treatment, LIMS_ID),
+ names_from = met,
+ values_from = fc) %>%
+ #arrange(Compound_Class, Compound_Name) %>%
+ as.data.frame()
+
+cor.GC_numeric <- cor.GC %>%
+ select(starts_with("m_"))
+
+g_vec <- genotypes$genotype %>%
+ sort() %>%
+ as.character()
+
+g_vec_lab <- genotypes %>%
+ arrange(genotype) %>%
+ select(genotype_label) %>%
+ as_vector()
+
+setwd(current)
+latest <- str_sort(str_extract(list.files(pattern = "^\\d{6}_correlation$"),
+ pattern = "^\\d{6}_correlation"),
+ decreasing = T)[[1]]
+
+latest_Correlation <- str_c(current, "/", latest)
+
+setwd(latest_Correlation)
+
+## Start -------------------------------------------------------------------
+
+
+for(g in seq_along(genotypes$genotype)) {
+ print(str_c("starting with", g_vec_lab[[g]], sep = " "))
+
+latest_Correlation_date <- str_extract(latest, pattern = "^\\d{6}")
+
+corr_mat <- read_csv(list.files()[str_detect(list.files(),
+ str_c("corr_mat_",
+ g_vec_lab[[g]],
+ ".csv"))])
+
+corr_mat_comb <- corr_mat %>%
+ filter(abs(cor_val) >= 0.75, p_adj <= 0.05) %>%
+ mutate(Correlation = as_factor(if_else(pos_neg == "pos", "+", "-"))) %>%
+ mutate(Correlation = fct_relevel(Correlation, c("+", "-")),
+ weight = abs(cor_val))
+
+#classes_tidy <- met_dat %>%
+# mutate(`ChEBI simplified` = if_else(str_detect(Compound_Class,"carbohydrate"),"carbohydrate or derivative",
+# if_else(str_detect(Compound_Class,"amino_acid"), "amino acid or derivative",
+# if_else(str_detect(Compound_Class,"carboxylic_acid"), "carboxylic acid",
+# "other"))))
+
+annotation_col_corr <- classes_tidy %>%
+ select(ChEBI) %>%
+ as.data.frame()
+
+rownames(annotation_col_corr) <- classes_tidy$met
+
+ann_colors_corr = list(
+ ChEBI= c("AA"= cb_scale[5],
+ "CH" = cb_scale[6],
+ "CA" = cb_scale[7],
+ "other" = cb_scale[8])
+)
+
+
+
+ann_colors_graph = tibble(ChEBI= names(ann_colors_corr$ChEBI),
+ color = ann_colors_corr$ChEBI)
+
+node_cols <- tibble(met = colnames(cor.GC_numeric)) %>%
+ left_join(classes_tidy) %>%
+ left_join(ann_colors_graph) %>%
+ mutate(ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other")))
+
+node_col_vec <- node_cols$color
+
+
+net_raw <- graph_from_data_frame(corr_mat_comb, vertices = node_cols, directed = F)
+
+set.seed(10000)
+
+if (g_vec_lab[[g]] == "MoneyMaker") {
+ layout<- layout_with_graphopt(net_raw , charge = 0.5,
+ spring.length = 1,
+ spring.constant = 100,
+ mass = 100)
+ #layout <- layout_with_fr(net_raw)
+
+ net_metric <- tibble(parameters = c("Number of nodes",
+ "Number of edges",
+ "Average node level centrality degree",
+ "Graph level centrality degree",
+ "Average node level centrality closeness",
+ "Average node level centrality betweenness",
+ "Graph level centrality betweenness",
+ "Graph diameter",
+ "Average distance"
+ ))
+
+}
+
+node_cols <- tibble(met = colnames(cor.GC_numeric)) %>%
+ left_join(classes_tidy) %>%
+ left_join(ann_colors_graph)%>%
+ mutate(ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other")))
+
+node_col_vec <- node_cols$color
+
+
+net_raw <- graph_from_data_frame(corr_mat_comb, vertices = node_cols, directed = F)
+
+ggnetwork(net_raw, layout = layout) %>%
+ mutate(ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other"))) %>%
+ ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
+ geom_edges(aes(color = Correlation),curvature = 0.15) +
+ geom_nodes(aes(fill = ChEBI), color = "black", shape = 21, size = 2) +
+ geom_nodetext(aes(label = Compound_Name), size = 2) +
+ #geom_nodetext(aes(label = if_else(Compound_Name == "malic_acid","Malic acid", ""))) +
+ theme_blank() +
+ scale_fill_manual(values = cb_scale[5:8], aesthetics = "fill") +
+ scale_color_manual(values = c(cb_scale[3], cb_scale[1]), aesthetics = "color") +
+ cor_theme +
+ theme(legend.position = c(0.5, 0.5))
+
+ggsave(str_c(g_vec_lab[[g]],
+ "_Correlation_network_wt_layout.pdf"),
+ width = 8.25,
+ height = 8.25,
+ units = "cm",
+ dpi = 300)
+
+ggnetwork(net_raw, layout = layout) %>%
+ mutate(ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other"))) %>%
+ ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
+ geom_edges(aes(color = Correlation),curvature = 0.15) +
+ geom_nodes(aes(fill = ChEBI), color = "black", shape = 21, size = 2) +
+ theme_blank() +
+ scale_fill_manual(values = cb_scale[5:8], aesthetics = "fill") +
+ scale_color_manual(values = c(cb_scale[3], cb_scale[1]), aesthetics = "color") +
+ cor_theme +
+ theme(legend.position = "right") +
+ scale_x_continuous(expand = expansion(mult = 0.02)) +
+ scale_y_continuous(expand = expansion(mult = 0.02))
+
+ggsave(str_c(g_vec_lab[[g]],
+ "_Correlation_network_wt_layout_unlab.pdf"),
+ width = 8.25,
+ height = 8.25,
+ units = "cm",
+ dpi = 300)
+
+saveRDS(last_plot(),str_c(g_vec_lab[[g]],
+ "_Correlation_network_wt_layout_unlab.RDS"))
+
+
+### Own layout and node filtered --------------------------------------------
+
+ node_cols <- tibble(met = colnames(cor.GC_numeric)) %>%
+ left_join(classes_tidy) %>%
+ left_join(ann_colors_graph) %>%
+ filter(met %in% corr_mat_comb$met1 | met %in% corr_mat_comb$met2) %>%
+ mutate(ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other")))
+
+ net_raw <- graph_from_data_frame(corr_mat_comb, vertices = node_cols, directed = F)
+
+ net_metric <- net_metric %>% add_column(.after = g,
+ "current" =
+ c(net_raw %>% length(),
+ edge.attributes(net_raw)$cor_val %>% length(),
+ degree(net_raw) %>% mean(),
+ centr_degree(net_raw)$centralization,
+ centr_clo(net_raw)$res %>% mean(na.rm = T), #returns NaN for some nodes
+ centr_betw(net_raw)$res %>% mean(),
+ centr_betw(net_raw)$centralization,
+ diameter(net_raw), # does not work with negative weight scores
+ mean_distance(net_raw) # does not work with negative weight scores
+ )
+ )
+
+ colnames(net_metric)[g+1] <- g_vec_lab[[g]]
+ set.seed(10000)
+ layout_own <- layout_with_fr(net_raw)
+
+ ggnetwork(net_raw, layout = layout_own) %>%
+ mutate(ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other"))) %>%
+ ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
+ geom_edges(aes(color = Correlation),curvature = 0.15) +
+ geom_nodes(aes(fill = ChEBI), color = "black", shape = 21, size = 2) +
+ geom_nodetext(aes(label = Compound_Name), size = 2) +
+ #geom_nodetext(aes(label = if_else(Compound_Name == "malic_acid","Malic acid", ""))) +
+ theme_blank() +
+ scale_fill_manual(values = cb_scale[5:8], aesthetics = "fill") +
+ scale_color_manual(values = c(cb_scale[3], cb_scale[1]), aesthetics = "color") +
+ cor_theme +
+ theme(legend.position = c(0.5, 0.5))
+
+ ggsave(str_c(g_vec_lab[[g]],
+ "_Correlation_network_own.pdf"),
+ width = 8.25,
+ height = 8.25,
+ units = "cm",
+ dpi = 300)
+
+ ggnetwork(net_raw, layout = layout_own) %>%
+ mutate(ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other"))) %>%
+ ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
+ geom_edges(aes(color = Correlation),curvature = 0.15) +
+ geom_nodes(aes(fill = ChEBI), color = "black", shape = 21, size = 2) +
+ #geom_nodetext(aes(label = Compound_Name), size = 2) +
+ #geom_nodetext_repel(aes(label = if_else(Compound_Name == "D-maltose","Malt", "")),inherit.aes = T) +
+ theme_blank() +
+ scale_fill_manual(values = cb_scale[5:8], aesthetics = "fill") +
+ scale_color_manual(values = c(cb_scale[3], cb_scale[1]), aesthetics = "color") +
+ cor_theme +
+ theme(legend.position = "right") +
+ scale_x_continuous(expand = expansion(mult = 0.02)) +
+ scale_y_continuous(expand = expansion(mult = 0.02))
+
+ ggsave(str_c(g_vec_lab[[g]],
+ "_Correlation_network_own_unlab.pdf"),
+ width = 8.25,
+ height = 8.25,
+ units = "cm",
+ dpi = 300)
+
+ saveRDS(last_plot(), str_c(g_vec_lab[[g]],
+ "_Correlation_network_own_unlab.RDS"))
+
+}
+
+## Loop end ----------------------------------------------------------------
+
+write_csv(net_metric, "Network_topology_metrics.csv")
+
+
+stop("current_end")
+
+
+# Test edgebundling -------------------------------------------------------
+
+library(edgebundle)
+test <- edge_bundle_force(net_raw, layout_own, compatibility_threshold = 0.2)
+
+ggplot(test) +
+ geom_path(aes(x = x, y = y, group = group)) +
+ geom_point(data = layout_own %>% as.data.frame(), aes(x = V1, y = V2), color = "red", size = 2)
+
+ggnetwork(net_raw, layout = test) %>%
+ ggplot(aes(x = x, y = y, group = group, xend = xend, yend = yend)) +
+ geom_edges(aes(color = Correlation),curvature = 0.15) +
+ geom_nodes(aes(fill = ChEBI), color = "black", shape = 21, size = 2) +
+ #geom_nodetext(aes(label = Compound_Name), size = 2) +
+ #geom_nodetext_repel(aes(label = if_else(Compound_Name == "D-maltose","Malt", "")),inherit.aes = T) +
+ theme_blank() +
+ scale_fill_manual(values = cb_scale[5:8], aesthetics = "fill") +
+ scale_color_manual(values = c(cb_scale[3], cb_scale[1]), aesthetics = "color") +
+ cor_theme +
+ theme(legend.position = "right") +
+ scale_x_continuous(expand = expansion(mult = 0.02)) +
+ scale_y_continuous(expand = expansion(mult = 0.02))
+
+# Test combine ------------------------------------------------------------
+
+#other methods
+
+set.seed(200)
+plot.igraph(net_raw, vertex.color = node_col_vec, vertex.label = "")
+
+ggraph(net_raw, layout = layout) +
+ geom_edge_arc0(aes(color = pos_neg), strength = 0.1) +
+ geom_node_point(aes(color = `ChEBI simplified`), size = 3)
+# geom_node_text(aes(label = name))
+
+corr_mat_plot <- corr_mat_long %>%
+ filter(abs(cor_val) >= 0.75) %>%
+ select(met1, met2) %>%
+ as.matrix()
+
+corr_mat_plot <- corr_mat
+corr_mat_plot[abs(corr_mat_plot)<0.5] <- 0
+
+net_raw <- graph_from_edgelist(corr_mat_plot, directed = F)
+net_raw <- graph_from_adjacency_matrix(corr_mat_plot,
+ weighted=T,
+ mode="undirected",
+ diag=F)
+
+stop <- T
+
+if (stop == T) {
+ stop(print("Automatically ended"))
+p1 <- make_box_dot_plot("m_81", "fruits", "Phenylalanine", "none", c("MoneyMaker","*log2-1*"))
+p1
+
+p2 <- make_box_dot_plot_levene("m_81", "fruits", "Phenylalanine", "none", c("MoneyMaker","*log2-1*"))
+p2
+
+p3 <- make_box_dot_plot_cv("m_81", "fruits", "Phenylalanine", "none", c("MoneyMaker","*log2-1*"))
+p3
+
+leg46 <- get_legend(make_box_dot_plot("m_81", "fruits", "Phenylalanine", "bottom", c("MoneyMaker","*log2-1*")))
+leg46
+
+main_plot <- plot_grid(p1,p2, nrow = 2, labels = "AUTO",rel_heights = c(1,1))
+main_plot
+
+main_plot_2 <- plot_grid(main_plot, p3, ncol = 2, labels = list("", "C"),rel_widths = c(3,1))
+main_plot_2
+
+comp_plot <- plot_grid(main_plot_2, leg46, ncol = 1, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Phenylalanine_combined.png",width = 15.8, height = 16, units = "cm", dpi = 300)
+
+#sugars?
+
+p1 <- make_box_dot_plot("m_72", "fruits", "Maltose", "none", c("MoneyMaker","*transp1-1*"))
+p1
+
+p2 <- make_box_dot_plot_levene("m_72", "fruits", "Maltose", "none", c("MoneyMaker","*transp1-1*"))
+p2
+
+p3 <- make_box_dot_plot_cv("m_72", "fruits", "Maltose", "none", c("MoneyMaker","*transp1-1*"))
+p3
+
+leg46 <- get_legend(make_box_dot_plot("m_72", "fruits", "Maltose", "bottom", c("MoneyMaker","*transp1-1*")))
+leg46
+
+main_plot <- plot_grid(p1,p2, nrow = 2, labels = "AUTO",rel_heights = c(1,1))
+main_plot
+
+main_plot_2 <- plot_grid(main_plot, p3, ncol = 2, labels = list("", "C"),rel_widths = c(3,1))
+main_plot_2
+
+
+ggsave("Maltose_combined.png",width = 15.8, height = 16, units = "cm", dpi = 300)
+
+
+
+# Recycling ---------------------------------------------------------------
+# Correlation -------------------------------------------------------------
+
+library(viridis)
+library(igraph)
+library(ggnetwork)
+library(jmuOutlier)
+
+cor_base <- fc_1_ind %>%
+ filter(#Compound_Class %in% c("Phospholipid", "Phosphatidylcholine"),
+ genotype %in% c("MoneyMaker"),
+ tissue == "fruits")
+
+cor.GC <- cor_base %>%
+ group_by(met) %>%
+ mutate(log_norm = log2(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(genotype, treatment, LIMS_ID),
+ names_from = met,
+ values_from = fc) %>%
+ #arrange(Compound_Class, Compound_Name) %>%
+ as.data.frame()
+
+cor.GC_numeric <- cor.GC %>%
+ select(starts_with("m_"))
+
+corr_mat <- cor(cor.GC_numeric, method = "spearman")
+
+for (i in seq_along(colnames(cor.GC_numeric))) {
+
+}
+corr_p <- perm.cor.test(x = cor.GC_numeric, method = "spearman")
+
+
+classes_tidy <- met_dat %>%
+ mutate(`ChEBI simplified` = if_else(str_detect(Compound_Class,"carbohydrate"),"carbohydrate or derivative",
+ if_else(str_detect(Compound_Class,"amino_acid"), "amino acid or derivative",
+ if_else(str_detect(Compound_Class,"carboxylic_acid"), "carboxylic acid",
+ "other"))))
+
+annotation_row_corr <- classes_tidy %>%
+ select(`ChEBI simplified`) %>%
+ as.data.frame()
+
+rownames(annotation_row_corr) <- classes_tidy$met
+
+annotation_col_corr <- classes_tidy %>%
+ select(`ChEBI simplified`) %>%
+ as.data.frame()
+
+rownames(annotation_col_corr) <- classes_tidy$met
+
+
+ann_colors_corr = list(
+ `ChEBI simplified`= c("amino acid or derivative"= vir_scale[1],
+ "carbohydrate or derivative" = vir_scale[2],
+ "carboxylic acid" = vir_scale[3],
+ "other" = vir_scale[4])
+)
+
+pheatmap <- pheatmap(corr_mat,
+ color = colorRampPalette(c("blue","white", "red"))(7),
+ cellwidth = 4,
+ cellheight = 4,
+ breaks = seq(-1.05, 1.05, 0.3),
+ cluster_rows = T,
+ cluster_cols = T,
+ annotation_names_row = F,
+ annotation_names_col = F,
+ show_rownames = F,
+ show_colnames = F,
+ annotation_row = annotation_row_corr,
+ annotation_col = annotation_col_corr,
+ annotation_legend =T,
+ legend = F,
+ display_numbers = F,
+ number_color = "black",
+ fontsize_number = 6,
+ annotation_colors = ann_colors_corr,
+ #filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ # "metabolite_correlation_heatmap_simple.jpg",
+ # sep = "_")
+)
+
+library(igraph)
+
+ann_colors_graph = tibble(`ChEBI simplified`= names(ann_colors_corr$`ChEBI simplified`),
+ color = ann_colors_corr$`ChEBI simplified`)
+
+node_cols <- tibble(met = colnames(cor.GC_numeric)) %>%
+ left_join(classes_tidy) %>%
+ left_join(ann_colors_graph)
+
+node_col_vec <- node_cols$color
+
+
+
+corr_mat_long <- corr_mat %>%
+ as_tibble(rownames = "met1") %>%
+ pivot_longer(cols = starts_with("m_"),
+ values_to = "cor_val",
+ names_to = "met2") %>%
+ mutate(forward = str_c(met1, met2, sep = "_"),
+ dup_check = str_c(sort(str_extract_all(forward,
+ "(?<=m_)\\d+",
+ simplify = T)),
+ collapse = "_")) %>%
+ filter(abs(cor_val) >= 0.5, !met1 == met2)
+
+corr_mat_plot <- corr_mat_long %>%
+ select(met1, met2) %>%
+ as.matrix()
+
+corr_mat_plot[abs(corr_mat_plot)<0.5] <- 0
+
+net_raw <- graph_from_edgelist(corr_mat_plot, directed = F)
+
+net_raw <- graph_from_adjacency_matrix(corr_mat_plot,
+ weighted=T,
+ mode="undirected",
+ diag=F)
+set.seed(200)
+plot.igraph(net_raw, vertex.color = node_col_vec, vertex.label = "")
+
+edges <- E(net_raw) #%>%
+names() %>%
+ tibble()
+
+V(net_raw)$Compound_Class <- node_cols$`ChEBI simplified`
+E(net_raw)$cor <- corr_mat_plot
+
+test <- layout.fruchterman.reingold(net_raw)
+
+ggraph(net_raw, layout = test) +
+ geom_edge_link0(aes(color = cor)) +
+ geom_node_point(aes(color = Compound_Class), size = 3)
+
+
+
+# Correlation network -------------------------------------------------------------
+
+cor_theme <- theme(axis.text = element_blank(),
+ axis.title = element_blank(),
+ panel.background = element_rect(fill = "white"),
+ strip.text = element_text(size = 10, family = "sans"),
+ text = element_text(size = 6, family = "sans"),
+ legend.position = "",
+ legend.title = element_text(size = 8, family = "sans"),
+ legend.text = element_text(size = 6, family = "sans"),
+ plot.margin = unit(c(1,1,1,1), "pt"),
+ legend.box.margin = margin(t = 0, r = 0, b = 0, l = 0 , unit = "mm"),
+ legend.margin = margin(t = 0, r = 0, b = 0, l = 0 , unit = "mm"))
+
+
+library(viridis)
+library(igraph)
+library(ggnetwork)
+library(jmuOutlier)
+library(ggraph)
+
+
+
+cor_base <- fc_1_ind %>%
+ filter(#Compound_Class %in% c("Phospholipid", "Phosphatidylcholine"),
+ genotype %in% c("MoneyMaker"),
+ tissue == "fruits")
+
+cor.GC <- cor_base %>%
+ group_by(met) %>%
+ mutate(log_norm = log2(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(genotype, treatment, LIMS_ID),
+ names_from = met,
+ values_from = fc) %>%
+ #arrange(Compound_Class, Compound_Name) %>%
+ as.data.frame()
+
+cor.GC_numeric <- cor.GC %>%
+ select(starts_with("m_"))
+
+corr_mat <- cor(cor.GC_numeric, method = "spearman")
+
+x <- vector(mode = "list", length = length(colnames(cor.GC_numeric)))
+names(x) <- colnames(cor.GC_numeric)
+
+start <- Sys.time()
+for (i in seq_along(colnames(cor.GC_numeric))) {
+ #names(x[[i]]) <- colnames(cor.GC_numeric)
+ for (n in seq_along(colnames(cor.GC_numeric))) {
+ if(i > n) {
+ x[[i]][[n]] <- NA
+ } else {
+ x[[i]][[n]] <- perm.cor.test(cor.GC_numeric[[i]], cor.GC_numeric[[n]],
+ method = "spearman",
+ num.sim = 100
+ )$p.value
+ }
+ }
+}
+
+stop <- Sys.time()
+stop-start
+
+for (i in seq_along(colnames(cor.GC_numeric))) {
+ names(x[[i]]) <- colnames(cor.GC_numeric)
+}
+
+corr_p <- x %>%
+ reduce(bind_rows)
+
+corr_p_long <- corr_p %>%
+ mutate(met1 = colnames(cor.GC_numeric)) %>%
+ pivot_longer(cols = starts_with("m_"),
+ values_to = "cor_p_val",
+ names_to = "met2") %>%
+ mutate(#cor_num = row_number(),
+ met1_num = as.numeric(str_extract(met1,"(?<=m_)\\d+")),
+ met2_num = as.numeric(str_extract(met2,"(?<=m_)\\d+"))) %>%
+ mutate(small = if_else(met1_num < met2_num, met1_num, met2_num),
+ big = if_else(met1_num > met2_num, met1_num, met2_num),
+ sorted = str_c(small, big, sep = "_")) %>%
+ filter(!is.na(cor_p_val)) %>%
+ #distinct(sorted, .keep_all = T) %>%
+ mutate(p_adj = p.adjust(cor_p_val, method = "fdr")) %>%
+ filter(!met1 == met2) %>%
+ # group_by(cor_num) %>%
+ mutate(cor_p_val_plot = if_else(p_adj > 0.05, 1, p_adj))
+
+corr_mat_long <- corr_mat %>%
+ as_tibble(rownames = "met1") %>%
+ pivot_longer(cols = starts_with("m_"),
+ values_to = "cor_val",
+ names_to = "met2") %>%
+ mutate(#cor_num = row_number(),
+ met1_num = as.numeric(str_extract(met1,"(?<=m_)\\d+")),
+ met2_num = as.numeric(str_extract(met2,"(?<=m_)\\d+"))) %>%
+ mutate(small = if_else(met1_num < met2_num, met1_num, met2_num),
+ big = if_else(met1_num > met2_num, met1_num, met2_num),
+ sorted = str_c(small, big, sep = "_")) %>%
+ distinct(sorted, .keep_all = T) %>%
+ filter(!met1 == met2) %>%
+ mutate(pos_neg = as_factor(if_else(cor_val > 0, "pos", "neg")),
+ cor_val_plot = if_else(abs(cor_val)< 0.4, 0, cor_val),
+ weight = cor_val)
+
+corr_mat_comb <- corr_mat_long %>%
+ left_join(corr_p_long) %>%
+ filter(abs(cor_val_plot) >= 0.4 , cor_p_val_plot <= 0.05) %>%
+ mutate(ChEBI = fct_relevel(ChEBI, c("AA", "CH", "CA", "other")))
+
+
+
+#
+# Clustering difference---------------------------------------------------------
+
+library(pheatmap)
+
+heat_base <- fc_1 %>%
+ filter(#Compound_Class %in% c("Phospholipid", "Phosphatidylcholine"),
+ genotype %in% c("MoneyMaker", "*panK4-1*"),
+ tissue == "fruits")
+
+heat.GC <- heat_base %>%
+ 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)),
+ log_norm_level = log_norm_level + 1) %>%
+ ungroup() %>%
+ mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")))
+
+heat_GC_WT <- heat.GC %>%
+ filter(genotype == "MoneyMaker") %>%
+ select(met, treatment, log_norm_level_wt = log_norm_level)
+
+heat_GC_diff <- heat.GC %>%
+ filter(genotype == "*panK4-1*") %>%
+ left_join(heat_GC_WT) %>%
+ mutate(log_norm_diff = log_norm_level_wt - log_norm_level)
+
+heat.GC <- heat_GC_diff %>%
+ pivot_wider(id_cols = c(Compound_Name, Compound_Class, met),
+ names_from = group,
+ values_from = log_norm_diff) %>%
+ arrange(Compound_Class, Compound_Name) %>%
+ as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>%
+ select(contains("0.4"), contains("0.6"),
+ contains("0.8"),
+ contains("1"),
+ contains("leaves"), contains("fruits")) %>% as.matrix()
+
+plot(hclust(dist(mat.heat.GC)))
+
+all_clust <- hclust(dist(mat.heat.GC, method = "euclidean"))
+
+clusters <- cutree(all_clust, k = 4)
+
+clust_tidy <- tibble("clust" = clusters) %>%
+ mutate(met = names(clusters))
+
+heat_GC_diff %>%
+ left_join(clust_tidy) %>%
+ 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() %>%
+ ggplot(aes(x = treatment, y = log_norm_diff, color = Compound_Class)) +
+ geom_line(aes(group = met)) +
+ facet_wrap(vars(clust))
+
+annotation_row <- heat.GC %>%
+ select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- heat_base %>%
+ 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 <- heat_base %>%
+ 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.value >= 0.05| is.na(p.value), "","X")) %>%
+ pivot_wider(id_cols = c(met),
+ names_from = group,
+ values_from = signif) %>%
+ left_join(met_dat) %>%
+ #left_join(GCid_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 = cb_scale[1], `*panK4-1*` = "brown", `log2-1` = "blue", `*transp1-1*` = cb_scale[4]),
+ 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"))(65),
+ #cellwidth = 16,
+ #cellheight = 4,
+ breaks = seq(-3.25, 3.25, 0.1),
+ clustering_distance_rows = dist((mat.heat.GC), method = "euclidean"),
+ cluster_rows = T,
+ cluster_cols = F,
+ 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,
+ #filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ # "cmQTL_val_1_heatmap_rel_tissue_wt.jpg",
+ # sep = "_")
+)
+
+dev.off()
+
+
+
+
+
+# Clustering ---------------------------------------------------------
+
+library(pheatmap)
+
+heat_base <- fc_1 %>%
+ filter(#Compound_Class %in% c("Phospholipid", "Phosphatidylcholine"),
+ genotype %in% c("MoneyMaker"),
+ tissue == "fruits")
+
+heat.GC <- heat_base %>%
+ 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)),
+ z_score = (mean_fc - mean(mean_fc))/(max(mean_fc)-min(mean_fc))) %>%
+ ungroup() %>%
+ mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+ pivot_wider(id_cols = c(Compound_Name, Compound_Class, met),
+ names_from = group,
+ values_from = z_score) %>%
+ arrange(Compound_Class, Compound_Name) %>%
+ as.data.frame()
+
+rownames(heat.GC) <- heat.GC$met
+
+mat.heat.GC <- heat.GC %>%
+ select(contains("0.4"), contains("0.6"),
+ contains("0.8"),
+ contains("1"),
+ contains("leaves"), contains("fruits")) %>% as.matrix()
+
+plot(hclust(dist(mat.heat.GC, method = "maximum"), method = "ward.D2"))
+
+all_clust <- hclust(dist(mat.heat.GC, method = "maximum"))
+
+clusters <- cutree(all_clust, k = 4)
+
+clust_tidy <- tibble("clust" = clusters) %>%
+ mutate(met = names(clusters))
+
+heat_base %>%
+ left_join(clust_tidy) %>%
+ 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)),
+ z_score = (mean_fc - mean(mean_fc))/(sd(mean_fc))) %>%
+ ungroup() %>%
+ ggplot(aes(x = treatment, y = z_score, color = Compound_Class)) +
+ geom_line(aes(group = met)) +
+ facet_wrap(vars(clust))
+
+annotation_row <- heat.GC %>%
+ select(Compound_Class)
+
+rownames(annotation_row) <- heat.GC$met
+
+annotation_col <- heat_base %>%
+ 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 <- heat_base %>%
+ 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.value >= 0.05| is.na(p.value), "","X")) %>%
+ pivot_wider(id_cols = c(met),
+ names_from = group,
+ values_from = signif) %>%
+ left_join(met_dat) %>%
+ #left_join(GCid_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 = cb_scale[1], `*panK4-1*` = "brown", `log2-1` = "blue", `*transp1-1*` = cb_scale[4]),
+ 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"))(65),
+ #cellwidth = 16,
+ #cellheight = 4,
+ breaks = seq(-3.25, 3.25, 0.1),
+ clustering_distance_rows = dist((mat.heat.GC), method = "euclidean"),
+ cluster_rows = T,
+ cluster_cols = F,
+ 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,
+ #filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ # "cmQTL_val_1_heatmap_rel_tissue_wt.jpg",
+ # sep = "_")
+)
+
+dev.off()
+
+
+
+
+testtheme <- theme(axis.text.x = element_markdown(angle = 45, hjust = 1, margin = unit(c(1,0,0,0), "mm")),
+ axis.title.x = element_blank(),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ text = element_text(size = 8, family = "Arial Narrow"),
+ legend.title = element_blank(),
+ legend.text = element_markdown(),
+ legend.position = "right",
+ legend.direction = "vertical",
+ plot.margin = unit(c(1,0,1,2), "mm"),
+ legend.margin = margin(t = 0, r = 2, b = 0, l = 2 , unit = "mm"))
+
+
+make_box_point_plot <- function(plot_met, plot_tissue, plot_label, plot_legend){
+
+ sig_bar <- sig_GC %>%
+ left_join(genotypes, by = c("group1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("group2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y) %>%
+ mutate(y.position = 1.1 * tot_val) %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ fc_1_ind_plot <- fc_1_ind %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ plot_out <- fc_1_ind %>%
+ filter(met == plot_met, tissue == plot_tissue) %>%
+ ggplot(aes(x = genotype, y = fc)) +
+ geom_boxplot(position = "dodge", aes(fill = genotype), color = "black") +
+ geom_point(aes(fill = genotype, color = genotype),size = 2, shape = 21) +
+ 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 = 10),
+ legend.text = element_markdown(),
+ legend.position = plot_legend) +
+ ylab("Mean fold-change")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ scale_color_manual(values = bw_scale, aesthetics = "color")
+
+ plot_out
+
+}
+
+make_box_beeswarm_plot <- function(plot_met, plot_tissue, plot_label, plot_legend){
+
+ sig_bar <- sig_GC %>%
+ left_join(genotypes, by = c("group1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("group2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y) %>%
+ mutate(y.position = 1.1 * tot_val) %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ fc_1_ind_plot <- fc_1_ind %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ plot_out <- fc_1_ind %>%
+ filter(met == plot_met, tissue == plot_tissue) %>%
+ ggplot(aes(x = genotype, y = fc)) +
+ geom_boxplot(position = "dodge", aes(fill = genotype), color = "black") +
+ geom_beeswarm(aes(fill = genotype, color = genotype),size = 2, shape = 21, cex = 2, priority = "density") +
+ 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 = 10),
+ legend.text = element_markdown(),
+ legend.position = plot_legend) +
+ ylab("Mean fold-change")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ scale_color_manual(values = bw_scale, aesthetics = "color")
+
+ plot_out
+
+}
+
+make_box_dot_line_plot <- function(plot_met, plot_tissue, plot_label, plot_legend, plot_genotypes){
+
+ sig_bar <- sig_GC %>%
+ left_join(genotypes, by = c("group1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("group2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y) %>%
+ mutate(y.position = 1.1 * tot_val,
+ treatment = as_factor(treatment)) %>%
+ filter(met == plot_met, tissue == plot_tissue,
+ group1 %in% plot_genotypes, group2 %in% plot_genotypes)
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ fc_1_plot <- fc_1%>%
+ filter(met == plot_met, tissue == plot_tissue, genotype %in% plot_genotypes) %>%
+ mutate(treatment = as_factor(treatment))
+
+ plot_out <- fc_1_ind %>%
+ filter(met == plot_met, tissue == plot_tissue, genotype %in% plot_genotypes) %>%
+ mutate(treatment = as_factor(treatment)) %>%
+ ggplot(position = "dodge", aes(x = treatment, y = fc)) +
+ geom_boxplot(aes(fill = genotype), color = "black") +
+ geom_dotplot(aes(fill = genotype, color = genotype), position = position_dodge(0.75),
+ binaxis = "y", stackdir = "center", binwidth = 0.1, dotsize = 1) +
+ geom_smooth(aes(x = treatment, y = fc, color = genotype, group = genotype),method = "lm", position = position_dodge(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 = 10),
+ legend.text = element_markdown(),
+ legend.position = plot_legend) +
+ ylab("Mean fold-change")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ scale_color_manual(values = bw_scale, aesthetics = "color")
+
+ plot_out
+
+}
+
+p1 <- make_box_dot_line_plot("m_70", "fruits", "Malic acid", "none", c("MoneyMaker","*panK4-1*"))
+p1
+ggsave("Malic_acid_box_dot_line_plot_nominal.png",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+make_col_plot_cv <- function(plot_met, plot_tissue, plot_label, plot_legend){
+
+ sig_bar <- sig_GC_cv %>%
+ left_join(genotypes, by = c("alias1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("alias2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y_cv) %>%
+ mutate(y.position = 1.1 * tot_val) %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ plot_out <- fc_cv %>%
+ filter(met == plot_met, tissue == plot_tissue) %>%
+ 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 = 10),
+ legend.text = element_markdown(),
+ legend.position = plot_legend) +
+ ylab("CV fold-change")+
+ ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+ plot_out
+
+}
+
+make_col_plot <- function(plot_met, plot_tissue, plot_label, plot_legend){
+
+ sig_bar <- sig_GC %>%
+ left_join(genotypes, by = c("group1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("group2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y) %>%
+ mutate(y.position = 1.1 * tot_val) %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ fc_1_ind_plot <- fc_1_ind %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ plot_out <- fc_1 %>%
+ filter(met == plot_met, tissue == plot_tissue) %>%
+ ggplot(aes(x = genotype, y = mean_fc)) +
+ geom_col(position = "dodge", aes(fill = genotype), color = "black") +
+ geom_jitter(aes(fill = genotype, color = genotype, y = fc), size = 2, shape = 21, data = fc_1_ind_plot) +
+ 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 = 10),
+ legend.text = element_markdown(),
+ legend.position = plot_legend) +
+ ylab("Mean fold-change")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ scale_color_manual(values = bw_scale, aesthetics = "color")
+
+ plot_out
+
+}
+
+# Figure 46 ---------------------------------------------------------------
+
+p1 <- make_col_plot("m_70", "fruits", "Malic acid", "none")
+p1
+
+
+p2 <- make_col_plot_cv("m_119", "fruits", "Malic acid CV", "none")
+p2
+
+
+leg46 <- get_legend(make_col_plot("m_123", "fruits", "sucrose", "bottom"))
+leg46
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg46, ncol = 1, rel_heights = c(10,1))
+comp_plot
+
+ggsave("figure_46.wmf",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+# Figure 47 ---------------------------------------------------------------
+
+p1 <- make_col_plot("m_74", "fruits", "Phenylalanine", "none")
+p1
+
+p2 <- make_col_plot_cv("m_74", "fruits", "Phe CV", "none")
+p2
+
+
+leg46 <- get_legend(make_col_plot("m_74", "fruits", "sucrose", "bottom"))
+leg46
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO",rel_widths = c(3,1))
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg46, ncol = 1, rel_heights = c(10,1))
+comp_plot
+
+ggsave("figure_47.wmf",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+# Figure 48 ---------------------------------------------------------------
+
+p1 <- make_col_plot_cv("m_25", "fruits", "F6P CV", "none")
+p1
+
+p2 <- make_col_plot_cv("m_38", "fruits", "G6P CV", "none")
+p2
+
+p3 <- make_col_plot_cv("m_65", "fruits", "Maltose CV", "none")
+p3
+
+leg48 <- get_legend(make_col_plot("m_65", "fruits", "sucrose", "bottom"))
+leg48
+
+main_plot <- plot_grid(p1,p2,p3, ncol = 3, labels = "AUTO")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg48, ncol = 1, rel_heights = c(10,1))
+comp_plot
+
+ggsave("figure_48.wmf",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+}
+
+
\ No newline at end of file
diff --git a/workflows/Phenotype_analysis/210114_cmQTL_val_phenotype_analysis.R b/workflows/Phenotype_analysis/210114_cmQTL_val_phenotype_analysis.R
new file mode 100644
index 0000000000000000000000000000000000000000..d6f45eba2d021d9cd41d76a942257d68683584c2
--- /dev/null
+++ b/workflows/Phenotype_analysis/210114_cmQTL_val_phenotype_analysis.R
@@ -0,0 +1,826 @@
+rm(list = ls())
+
+
+# Libraries ---------------------------------------------------------------
+
+library(tidyverse)
+library(openxlsx)
+library(lubridate)
+library(broom)
+library(ggtext)
+library(glue)
+library(ggpubr)
+library(markdown)
+library(cowplot)
+library(RVAideMemoire)
+
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+getwd()
+
+pt <- readxl::read_xlsx("../200923_samplelist_WIJESI-030820-13_cmQTL_validation.xlsx")
+
+line_info <- readxl::read_xlsx("../WIJESI-030820-13- Culturereport-051220-7624.xlsx")
+
+genotypes <- readxl::read_xlsx("../Genotype_names.xlsx")
+
+pt_tidy <- pt %>%
+ mutate(col_pot = as.integer(col_pot),
+ comment = if_else(is.na(comment), "NA", comment),
+ irrigation = as_factor(irrigation),
+ #LIMS_ID = as.double(str_replace(LIMS_ID, "(?<=244)1", "")),
+ fruit_num_dropped = if_else(is.na(fruit_num_dropped), 0, fruit_num_dropped),
+ weight_sartorius_dropped = if_else(is.na(weight_sartorius_dropped), 0, weight_sartorius_dropped),
+ total_fruit_num = fruit_num_dropped + fruit_num_harvest,
+ total_fruit_weight = weight_sartorius_dropped + weight_sartorius_harvest,
+ #total_fruit_weight = if_else(str_detect(comment, "beefsteak"), NaN, total_fruit_weight),
+ av_fruit_weight = total_fruit_weight/total_fruit_num,
+ biomass = total_fruit_weight + shoot_fresh_weight,
+ harvest_index = total_fruit_weight/biomass,
+ days_until_flower = sow_date %--% flower_date %/% days(1),
+ days_until_redripe = sow_date %--% redripe_date %/% days(1)) %>%
+ left_join(line_info, by = c("LIMS_ID" = "aliquot")) %>%
+ left_join(genotypes) %>%
+ mutate(plantline = as_factor(plantline),
+ genotype = fct_relevel(genotype, c("MoneyMaker", "panK4-1", "log2-1", "transp1-1")),
+ irrigation = fct_relevel(irrigation, c("0.4", "0.6", "0.8", "1")),
+ group = str_c(genotype, irrigation, sep = "_"))
+
+pt_tidy %>%
+ mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+ genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*"))) %>%
+ ggplot(aes(x = col_pot, y = irrigation, color = genotype)) +
+ geom_point( size = 8) +
+ xlab("Column") +
+ theme(text = element_text(size = 18),
+ legend.text = element_markdown())
+
+#ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+# "cmQTL_val1_layout.jpg",
+# sep = "_"),
+# width = 27,
+# height = 9,
+# units = "cm",
+# dpi = 300)
+
+pt_vars <- c(#"fruit_num_dropped", "weight_sartorius_dropped",
+ "total_fruit_num", "total_fruit_weight",
+ "av_fruit_weight",# "days_until_flower",
+ "days_until_redripe","shoot_fresh_weight", "harvest_index",
+ "biomass", "shoot_dry_weight")
+
+pt_long <- pt_tidy %>%
+ pivot_longer(cols = all_of(pt_vars),
+ names_to = "pt_var",
+ values_to = "value") %>%
+ group_by(plantline, alias, genotype, irrigation, pt_var) %>%
+ mutate(rep = row_number()) %>%
+ filter(!is.na(value)) %>%
+ arrange(pt_var) %>%
+ ungroup()
+
+pt_long %>%
+ group_by(plantline, alias, irrigation, pt_var) %>%
+ ggplot(aes(irrigation, value, color = genotype)) +
+ geom_boxplot() +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+pt_long %>%
+ group_by(plantline, alias, irrigation, pt_var) %>%
+ ggplot(aes(irrigation, value, fill = alias)) +
+ #stat_summary(aes(color = alias)) +
+ #geom_boxplot(position = "identity")+
+ #geom_dotplot(binaxis = "y", stackdir = "center") +
+ geom_smooth(aes(group = alias, color = alias), method = "lm", formula = y ~ poly(x,2), se = T) +
+ facet_wrap(facets = vars(pt_var),scales = "free", nrow = 2)
+
+WT <- pt_long %>%
+ filter(plantline == "967514") %>%
+ rename(WT_value = value) %>%
+ select(irrigation, pt_var, WT_value, rep)
+
+# Summary statistics ------------------------------------------------------
+
+pt_vars_label <- c("Total fruit number per plant", "Total fruit weight per plant/ g",
+ "Average fruit weight/ g", "Days until ripe fruit", "Shoot fresh weight/ g",
+ "Harvest index", "Biomass", "Shoot dry weight/ g")
+
+pt_vars_plot <- cbind(pt_vars, pt_vars_label) %>% as_tibble()
+
+sum_stat <- pt_long %>%
+ left_join(WT, by = c("irrigation", "pt_var", "rep")) %>%
+ group_by(plantline, genotype, alias, irrigation, pt_var) %>%
+ summarise(mean = mean(value),
+ median = median(value),
+ sd = sd(value),
+ n = n(),
+ se = sd/sqrt(n),
+ ttest = t.test(value, WT_value)$p.value,
+ iqr = IQR(value),
+ Q1 = quantile(value)["25%"],
+ Q3 = quantile(value)["75%"]) %>%
+ ungroup() %>%
+ group_by(irrigation, pt_var) %>%
+ mutate(p_adj = p.adjust(ttest, method = "fdr"),
+ pt_var = as_factor(pt_var),
+ alias = as_factor(alias),
+ irrigation = as_factor(irrigation),
+ genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+ genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*"))) %>%
+ ungroup()
+
+
+
+#dry weight biomass still missing
+#AOV TukeyHSD?
+#test manova
+
+pt_tidy_numeric <- pt_tidy %>%
+ select(all_of(pt_vars))
+
+pt_norm <- pt_long %>%
+ group_by(irrigation, genotype, pt_var) %>%
+ summarise(norm = shapiro.test(value)$p.value)
+
+pt_norm_wt <- pt_norm %>%
+ filter(genotype == "MoneyMaker") %>%
+ select(irrigation, comp = genotype, pt_var, wt_norm = norm)
+
+pt_norm_mut <- pt_norm %>%
+ filter(genotype != "MoneyMaker") %>%
+ select(irrigation, genotype, pt_var, norm)
+
+pt_norm_comb <- pt_norm_wt %>%
+ left_join(pt_norm_mut)
+
+pt_variance <- map(.x = pt_tidy_numeric, .f = ~pairwise.var.test(.x, pt_tidy$group, p.method = "none"))
+
+pt_variance_tidy <- pt_variance %>%
+ map(tidy) %>%
+ map2(.y = names(pt_variance), .f = ~.x %>% mutate(pt_vars = .y)) %>%
+ purrr::reduce(bind_rows) %>%
+ separate(group1,
+ into = c("genotype1", "irrigation1"),
+ sep = "_") %>%
+ separate(group2,
+ into = c("genotype2", "irrigation2"),
+ sep = "_") %>%
+ filter(irrigation1 == irrigation2, genotype1 == "MoneyMaker"& genotype2 != "MoneyMaker"|
+ genotype1 != "MoneyMaker" & genotype2 == "MoneyMaker") %>%
+ mutate(genotype = if_else(genotype1 == "MoneyMaker", genotype2, genotype1),
+ comp = "MoneyMaker") %>%
+ select(genotype, comp, irrigation = irrigation1, pt_var = pt_vars, p.value)
+
+prereq <- pt_norm_comb %>%
+ left_join(pt_variance_tidy) %>%
+ mutate(fulfill = if_else(wt_norm >= 0.05, norm >= 0.05, p.value >= 0.05))
+
+prereq_sum <- prereq %>%
+ group_by(fulfill) %>%
+ summarise(n = n()) %>%
+ ungroup() %>%
+ mutate(percent = 100*n/sum(n))
+
+pt_aov <- map(.x = pt_tidy_numeric, .f = ~aov(.x ~ alias* irrigation, data = pt_tidy))
+
+pt_aov_tidy <- pt_aov %>%
+ map(tidy) %>%
+ map2(.y = names(pt_aov), .f = ~.x %>% mutate(pt_vars = .y)) %>%
+ purrr::reduce(bind_rows)
+
+pt_aov_table <- pt_aov_tidy %>%
+ left_join(pt_vars_plot) %>%
+ filter(term != "Residuals") %>%
+ select(term, pt_vars_label, p.value) %>%
+ pivot_wider(names_from = term, values_from = p.value)
+
+write.xlsx(pt_aov_table, str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "anova_table.xlsx",
+ sep = "_"))
+
+pt_tuk <- map(.x = pt_aov, .f = ~TukeyHSD(.x)) %>%
+ map(.f = tidy) %>%
+ map2(.y = names(pt_aov), .f = ~.x %>% mutate(var = .y)) %>%
+ purrr::reduce(bind_rows)
+
+# General differences -----------------------------------------------------
+
+sig_pt_overall <- pt_tuk %>%
+ filter(term == "alias") %>%
+ separate(col = contrast, into = c("group1", "group2"), sep = "-")
+
+lines <- sum_stat %>% distinct(alias,plantline) %>%
+ mutate(plantline = as.double(as.character(plantline))) %>%
+ left_join(genotypes)
+
+sig_pt_overall_table <- sig_pt_overall %>%
+ filter(group1 == "967514 MM WT" | group2 == "967514 MM WT") %>%
+ ungroup() %>%
+ left_join(lines, by = c("group1" = "alias")) %>%
+ select(var, group1, group2, adj.p.value, genotype1 = genotype) %>%
+ left_join(lines, by = c("group2" = "alias")) %>%
+ select(var,adj.p.value, genotype1, genotype2 = genotype) %>%
+ left_join(pt_vars_plot, by = c("var" = "pt_vars")) %>%
+ pivot_wider(names_from = genotype2,
+ values_from = adj.p.value)
+
+write.xlsx(sig_pt_overall_table, str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "overall_significance.xlsx",
+ sep = "_"))
+
+# Pairwise significant changes --------------------------------------------
+
+
+sig_pt_groups <- pt_tuk %>%
+ filter(term == "alias:irrigation") %>%
+ separate(col = contrast, into = c("group1", "group2"), sep = "-")
+
+sig_pt <- sig_pt_groups %>%
+ separate(group1, into = c("alias1", "irrigation1"), sep = ":") %>%
+ separate(group2, into = c("alias2", "irrigation2"), sep = ":") %>%
+ filter(irrigation1 == irrigation2) %>%
+ filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT") %>%
+ ungroup() %>%
+ mutate(p.signif = if_else(adj.p.value <= 0.05, "*", "ns")) %>%
+ left_join(sum_stat, by = c("var" = "pt_var", "irrigation1" = "irrigation", "alias1" = "alias")) %>%
+ select(var, irrigation1, irrigation2, alias1, alias2, p.signif, mean1 = mean, se1 = se, genotype1 = genotype) %>%
+ left_join(sum_stat, by = c("var" = "pt_var", "irrigation2" = "irrigation", "alias2" = "alias")) %>%
+ select(var, irrigation1, irrigation2, alias1, alias2, p.signif, mean1, mean2 = mean, se1, se2 = se, genotype1, genotype2 = genotype) %>%
+ group_by(irrigation1, var) %>%
+ mutate(tot_val1 = mean1 + se1,
+ tot_val2 = mean2 + se2,
+ y.position = 1.1*(max(tot_val1, tot_val2))) %>%
+ rename(pt_var = var,
+ group1 = genotype1,
+ group2 = genotype2,
+ irrigation = irrigation1) %>%
+ mutate(pt_var = as_factor(pt_var),
+ group1 = as_factor(group1),
+ group2 = as_factor(group2),
+ irrigation = as_factor(irrigation)) %>%
+ ungroup()
+
+pt_vars
+
+plot_out <- vector("list", length = length(pt_vars))
+
+for(pt in seq_along(pt_vars)) {
+ sig_pt_temp <- sig_pt %>%
+ filter(pt_var == pt_vars[[pt]])
+
+plot_out[[pt]] <- sum_stat %>%
+ filter(pt_var == pt_vars[[pt]]) %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.75) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "right") +
+ ggtitle(pt_vars_label[[pt]]) +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+}
+
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_phenotypes.pdf"),
+ width = 9,
+ height = 4)
+
+for (pt in seq_along(pt_vars)) {
+ print(plot_out[[pt]])
+}
+
+dev.off()
+
+
+
+sum_stat %>%
+ ggplot(aes(irrigation, mean, color = alias)) +
+ geom_point() +
+ geom_line(aes(group = plantline)) +
+ geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd)) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+
+
+sum_stat %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ ggplot(aes(irrigation, mean, color = genotype)) +
+ geom_line(aes(group = genotype)) +
+ geom_pointrange(aes(ymin = mean - se, ymax = mean + se, shape = genotype, fill = genotype)) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ 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 = 8),
+ legend.text = element_markdown(),
+ legend.position = "right") +
+ scale_shape_manual(values = c(1,2,22,23))
+
+sum_stat %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ ggplot(aes(irrigation, mean, color = genotype)) +
+ geom_line(aes(group = genotype)) +
+ geom_point(aes(shape = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se, color = genotype, width = 0.25)) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ 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 = 8),
+ legend.text = element_markdown(),
+ legend.position = "right") +
+ scale_shape_manual(values = c(1,2,22,23))
+
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "cmQTL_val1_phenotypic.jpg",
+ sep = "_"),
+ width = 183,
+ height = 100,
+ units = "mm",
+ dpi = 300)
+
+
+# Overview for supplemental -----------------------------------------------
+
+pt_long %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ group_by(plantline, genotype, irrigation, pt_var) %>%
+ ggplot(aes(irrigation, value, fill = genotype)) +
+ #stat_summary(aes(color = genotype)) +
+ #geom_violin(aes(color = genotype),fill = NA, position = "identity", trim = T, scale = "width", alpha = 1)+
+ geom_dotplot(binaxis = "y", stackdir = "center") +
+ geom_smooth(aes(group = genotype, color = genotype), method = "lm", formula = y ~ poly(x,2), se = F, size = 1.5) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ 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 = 8),
+ legend.text = element_markdown(),
+ legend.position = "right") +
+ scale_shape_manual(values = c(1,2,22,23))
+
+
+
+sum_stat %>%
+ ggplot(aes(irrigation, mean, color = alias)) +
+ geom_point() +
+ geom_line(aes(group = plantline)) +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "cmQTL_val1_phenotypic.jpg",
+ sep = "_"),
+ width = 183,
+ height = 100,
+ units = "mm",
+ dpi = 300)
+
+library(bootstrap)
+
+theta <- function(x){
+ sd(x)/mean(x)
+}
+
+CV <- sum_stat %>%
+ ungroup() %>%
+ group_by(plantline, genotype, alias, pt_var) %>%
+ summarise(grand_mean = mean(mean),
+ grand_sd = sd(mean)) %>%
+ mutate(cv = grand_sd/grand_mean)
+
+cv_jack <- sum_stat %>%
+ group_by(plantline, genotype, alias, pt_var) %>%
+ summarise(cv = jackknife(mean, theta)$jack.values) %>%
+ mutate(jack_rep = row_number()) %>%
+ ungroup()
+
+cv_jack_wide <- cv_jack %>%
+ pivot_wider(id_cols = c(plantline, genotype, alias, pt_var, jack_rep),
+ names_from = pt_var,
+ values_from = cv)
+
+cv_jack_numeric <- cv_jack_wide %>%
+ select(-c(alias, genotype, plantline, jack_rep))
+
+cv_jack_sum_stat <- cv_jack %>%
+ group_by(plantline, genotype, alias, pt_var) %>%
+ summarise(mean_cv = mean(cv),
+ sd = sd(cv),
+ n = n(),
+ se = sd/sqrt(n))
+
+
+pt_cv_jack <- map(.x = cv_jack_numeric, .f = ~aov(.x ~ alias, data = cv_jack_wide))
+
+pt_jack_t <- map(.x = cv_jack_numeric, .f = ~pairwise.t.test(x = .x , g = cv_jack_wide$alias)) %>%
+ map(.f = tidy) #%>%
+
+pt_jack_t_tidy <- pt_jack_t %>%
+ map2(.y = names(pt_jack_t), .f = ~.x %>% mutate(var = .y)) %>%
+ purrr::reduce(bind_rows)
+
+sig_pt_cv_groups <- pt_jack_t_tidy %>%
+ filter(group1 == "967514 MM WT" | group2 == "967514 MM WT") %>%
+ 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_sum_stat, by = c("var" = "pt_var", "group1" = "alias")) %>%
+ select(var, group1, group2, p.signif, mean_cv1 = mean_cv, se1 = se, genotype1 = genotype) %>%
+ left_join(cv_jack_sum_stat, by = c("var" = "pt_var", "group2" = "alias")) %>%
+ select(var, group1, group2, p.signif, mean_cv1, mean_cv2 = mean_cv, se1, se2 = se, genotype1, genotype2 = genotype) %>%
+ group_by(var) %>%
+ mutate(tot_val1 = mean_cv1 + se1,
+ tot_val2 = mean_cv2 + se2,
+ y.position = 1.1*(max(tot_val1, tot_val2))) %>%
+ rename(pt_var = var) %>%
+ mutate(pt_var = as_factor(pt_var),
+ group1 = genotype1,
+ group2 = genotype2) %>%
+ ungroup()
+
+cv_jack %>%
+ ggplot(aes(alias, cv, fill = alias)) +
+ geom_boxplot() +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+
+
+
+CV %>%
+ ggplot(aes(alias, cv, fill = alias)) +
+ geom_col() +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+x = rnorm(10)
+
+
+
+# Jackknife plots---------------------------------------------------------------
+
+
+
+
+pt_vars
+pt_vars_label <- c("Total fruit number per plant", "Total fruit weight per plant/ g",
+ "Average fruit weight/ g", "Days until ripe fruit", "Shoot fresh weight/ g",
+ "Harvest index", "Biomass", "Shoot dry weight/ g")
+
+plot_out <- vector("list", length = length(pt_vars))
+
+for(pt in seq_along(pt_vars)) {
+ sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == pt_vars[[pt]])
+
+ plot_out[[pt]] <- cv_jack_sum_stat %>%
+ filter(pt_var == pt_vars[[pt]]) %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.75) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "right") +
+ ggtitle(pt_vars_label[[pt]]) +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV")
+}
+
+plot_out[[1]]
+
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_phenotypes_cv.pdf"),
+ width = 9,
+ height = 4)
+
+
+
+
+
+
+for (pt in seq_along(pt_vars)) {
+ print(plot_out[[pt]])
+}
+
+dev.off()
+
+
+
+plot(bootstrap::jackknife(x = x, theta = theta)$jack.values)
+
+# Figure 40 ---------------------------------------------------------------
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "total_fruit_weight")
+
+p1<- sum_stat %>%
+ filter(pt_var == "total_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "",
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ xlab("") +
+ ylim(0,1150) +
+ ggtitle("Total fruit weight per plant/ g") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+p1
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "av_fruit_weight")
+
+p2<- sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "",
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ xlab("") +
+ ylim(c(0,65))+
+ ggtitle("Average fruit weight/ g") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+p2
+
+leg40 <- get_legend(sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "right") +
+ ggtitle("Average fruit weight/ g") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill"))
+leg40
+
+main_plot <- plot_grid(p1,p2, ncol = 1, labels = "AUTO")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg40, ncol = 2, rel_widths = c(3,1))
+comp_plot
+
+ggsave("figure_40.wmf",width = 15.8, height = 12, units = "cm", dpi = 300)
+
+# Figure 41 ---------------------------------------------------------------
+
+
+sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == "av_fruit_weight")
+
+p1 <- cv_jack_sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "",
+ plot.title = element_text(size = 14),
+ plot.margin = unit(c(0,0,1,0), "mm")) +
+ ggtitle("Average fruit weight/ g") +
+ xlab("") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV")
+p1
+
+sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == "harvest_index")
+
+p2 <- cv_jack_sum_stat %>%
+ filter(pt_var == "harvest_index") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "",
+ plot.title = element_text(size = 14),
+ plot.margin = unit(c(0,0,1,0), "mm")) +
+ ggtitle("Harvest index") +
+ xlab("") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV")
+p2
+
+leg41 <- get_legend(cv_jack_sum_stat %>%
+ filter(pt_var == "harvest_index") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "bottom") +
+ ggtitle("Harvest index") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV"))
+
+
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg41, ncol = 1, rel_heights = c(10,1)) + theme(plot.margin = unit(c(0,0,5,0), "mm"))
+comp_plot
+
+ggsave("figure_41.wmf",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+
+# Figure means presentation -----------------------------------------------
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "total_fruit_weight")
+
+p1<- sum_stat %>%
+ filter(pt_var == "total_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "",
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ xlab("") +
+ ylim(0,1150) +
+ ggtitle("Total fruit weight per plant/ g") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+p1
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "av_fruit_weight")
+
+p2<- sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "",
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ xlab("") +
+ ylim(c(0,65))+
+ ggtitle("Average fruit weight/ g") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+p2
+
+leg40 <- get_legend(sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "right") +
+ ggtitle("Average fruit weight/ g") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill"))
+leg40
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg40, ncol = 2, rel_widths = c(6,1))
+comp_plot
+
+ggsave("Av_tot_fruit_weight_mean.wmf",width = 33.5, height = 12, units = "cm", dpi = 300)
+
+# Figure CV average fruit weight ---------------------------------------------------------------
+
+
+sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == "av_fruit_weight")
+
+p1 <- cv_jack_sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "",
+ plot.title = element_text(size = 14),
+ plot.margin = unit(c(0,0,0,0), "mm")) +
+ ggtitle("Average fruit weight/ g") +
+ xlab("") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV")
+p1
+
+
+
+leg41 <- get_legend(cv_jack_sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "right") +
+ ggtitle("Harvest index") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV"))
+
+
+
+main_plot <- plot_grid(p1, ncol = 1, labels = "")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg41, ncol = 2, rel_widths = c(3,2)) + theme(plot.margin = unit(c(0,0,0,0), "mm"))
+comp_plot
+
+ggsave("av_fruit_weight_cv_presentation.wmf",width = 10, height = 10, units = "cm", dpi = 300)
+
+
diff --git a/workflows/Phenotype_figures/Figure maker_phenotypes.R b/workflows/Phenotype_figures/Figure maker_phenotypes.R
new file mode 100644
index 0000000000000000000000000000000000000000..ca17cf9dc3867f77104d0a2a87386cb22a1ff0ca
--- /dev/null
+++ b/workflows/Phenotype_figures/Figure maker_phenotypes.R
@@ -0,0 +1,1448 @@
+
+# Library and Directory ---------------------------------------------------
+
+
+rm(list = ls())
+library(tidyverse)
+library(openxlsx)
+library(lubridate)
+library(broom)
+library(ggtext)
+library(glue)
+library(ggpubr)
+library(markdown)
+library(cowplot)
+library(ggbeeswarm)
+library(extrafont)
+library(devEMF)
+library(RVAideMemoire)
+
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+getwd()
+
+# Directory setting -------------------------------------------------------
+
+
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+getwd()
+
+current <- getwd()
+source <- str_c(current,"/..")
+
+cur_date <- str_c(str_replace_all(Sys.Date(),"^.{2}|-",""))
+
+out <- str_c(cur_date, "Figures", sep = "_")
+
+if (file.exists(out)) {
+ cat("The folder already exists")
+} else {
+ dir.create(out)
+}
+
+out_dir <- str_c(current, out, sep = "/")
+
+setwd(source)
+
+# Data loading ------------------------------------------------------------
+
+
+
+pt <- readxl::read_xlsx("200923_samplelist_WIJESI-030820-13_cmQTL_validation.xlsx")
+
+line_info <- readxl::read_xlsx("WIJESI-030820-13- Culturereport-051220-7624.xlsx")
+
+genotypes <- readxl::read_xlsx("Genotype_names.xlsx")
+
+setwd(out_dir)
+# Data combination --------------------------------------------------------
+
+com_theme <- theme(axis.text.x = element_markdown(angle = 45, hjust = 1, size = 8, family = "sans"),
+ axis.text.y = element_text(size = 8, family = "sans"),
+ axis.title.x = element_blank(),
+ axis.title.y = element_text(size = 8, family = "sans"),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ strip.text = element_text(size = 8, family = "sans"),
+ text = element_text(size = 8, family = "sans"),
+ legend.title = element_blank(),
+ legend.text = element_markdown(size = 8),
+ plot.margin = unit(c(1,1,1,1), "mm"),
+ legend.margin = margin(t = 0, r = 2, b = 0, l = 2 , unit = "mm"))
+
+
+cb_scale <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442",
+ "#0072B2", "#D55E00", "#CC79A7","#000000")
+bw_scale <- c("black", "black", "black", "black", "black", "black", "black")
+
+
+pt_tidy <- pt %>%
+ mutate(col_pot = as.integer(col_pot),
+ comment = if_else(is.na(comment), "NA", comment),
+ irrigation = as_factor(irrigation),
+ #LIMS_ID = as.double(str_replace(LIMS_ID, "(?<=244)1", "")),
+ fruit_num_dropped = if_else(is.na(fruit_num_dropped), 0, fruit_num_dropped),
+ weight_sartorius_dropped = if_else(is.na(weight_sartorius_dropped), 0, weight_sartorius_dropped),
+ total_fruit_num = fruit_num_harvest, #fruit_num_dropped,
+ total_fruit_weight = weight_sartorius_harvest, #weight_sartorius_dropped,
+ #total_fruit_weight = if_else(str_detect(comment, "beefsteak"), NaN, total_fruit_weight),
+ av_fruit_weight = total_fruit_weight/total_fruit_num,
+ biomass = total_fruit_weight + shoot_fresh_weight,
+ harvest_index = total_fruit_weight/biomass,
+ days_until_flower = sow_date %--% flower_date %/% days(1),
+ days_until_redripe = sow_date %--% redripe_date %/% days(1)) %>%
+ left_join(line_info, by = c("LIMS_ID" = "aliquot")) %>%
+ left_join(genotypes) %>%
+ mutate(plantline = as_factor(plantline),
+ genotype = fct_relevel(genotype, c("MoneyMaker", "panK4-1", "log2-1", "transp1-1")),
+ irrigation = fct_relevel(irrigation, c("0.4", "0.6", "0.8", "1")),
+ group = str_c(genotype, irrigation, sep = "_"))
+
+pt_tidy %>%
+ mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+ genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*"))) %>%
+ ggplot(aes(x = col_pot, y = irrigation, color = genotype)) +
+ geom_point( size = 8) +
+ xlab("Column") +
+ theme(text = element_text(size = 18),
+ legend.text = element_markdown())
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "cmQTL_val1_layout.jpg",
+ sep = "_"),
+ width = 27,
+ height = 9,
+ units = "cm",
+ dpi = 300)
+
+pt_vars <- c(#"fruit_num_dropped", "weight_sartorius_dropped",
+ "total_fruit_num", "total_fruit_weight",
+ "av_fruit_weight",# "days_until_flower",
+ "days_until_redripe","shoot_fresh_weight", "harvest_index",
+ "biomass", "shoot_dry_weight")
+
+pt_long <- pt_tidy %>%
+ pivot_longer(cols = all_of(pt_vars),
+ names_to = "pt_var",
+ values_to = "value") %>%
+ group_by(plantline, alias, genotype, irrigation, pt_var) %>%
+ mutate(rep = row_number()) %>%
+ filter(!is.na(value)) %>%
+ arrange(pt_var) %>%
+ ungroup()
+
+pt_long %>%
+ group_by(plantline, alias, irrigation, pt_var) %>%
+ ggplot(aes(irrigation, value, color = genotype)) +
+ geom_boxplot() +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+pt_long %>%
+ group_by(plantline, alias, irrigation, pt_var) %>%
+ ggplot(aes(irrigation, value, fill = alias)) +
+ #stat_summary(aes(color = alias)) +
+ #geom_boxplot(position = "identity")+
+ #geom_dotplot(binaxis = "y", stackdir = "center") +
+ geom_smooth(aes(group = alias, color = alias), method = "lm", formula = y ~ poly(x,2), se = T) +
+ facet_wrap(facets = vars(pt_var),scales = "free", nrow = 2)
+
+WT <- pt_long %>%
+ filter(plantline == "967514") %>%
+ rename(WT_value = value) %>%
+ select(irrigation, pt_var, WT_value, rep)
+
+# Summary statistics ------------------------------------------------------
+
+pt_vars_label <- c("Fruit number per plant", "Total fruit weight",
+ "Average fruit weight", "Days until ripe fruit", "Shoot fresh weight",
+ "Harvest index", "Biomass", "Shoot dry weight")
+
+pt_vars_plot <- cbind(pt_vars, pt_vars_label) %>% as_tibble()
+
+sum_stat <- pt_long %>%
+ left_join(WT, by = c("irrigation", "pt_var", "rep")) %>%
+ group_by(plantline, genotype, alias, irrigation, pt_var) %>%
+ summarise(mean = mean(value),
+ median = median(value),
+ sd = sd(value),
+ n = n(),
+ se = sd/sqrt(n),
+ ttest = t.test(value, WT_value)$p.value,
+ iqr = IQR(value),
+ Q1 = quantile(value)["25%"],
+ Q3 = quantile(value)["75%"]) %>%
+ ungroup() %>%
+ group_by(irrigation, pt_var) %>%
+ mutate(p_adj = p.adjust(ttest, method = "fdr"),
+ pt_var = as_factor(pt_var),
+ alias = as_factor(alias),
+ irrigation = as_factor(irrigation),
+ genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+ genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*"))) %>%
+ ungroup()
+
+
+# Estimation of normality and equal variances -----------------------------
+
+#dry weight biomass still missing
+#AOV TukeyHSD?
+#test manova
+
+pt_tidy_numeric <- pt_tidy %>%
+ select(all_of(pt_vars))
+
+pt_norm <- pt_long %>%
+ group_by(irrigation, genotype, pt_var) %>%
+ summarise(norm = shapiro.test(value)$p.value)
+
+pt_norm_wt <- pt_norm %>%
+ filter(genotype == "MoneyMaker") %>%
+ select(irrigation, comp = genotype, pt_var, wt_norm = norm)
+
+pt_norm_mut <- pt_norm %>%
+ filter(genotype != "MoneyMaker") %>%
+ select(irrigation, genotype, pt_var, norm)
+
+pt_norm_comb <- pt_norm_wt %>%
+ left_join(pt_norm_mut)
+
+pt_variance <- map(.x = pt_tidy_numeric, .f = ~pairwise.var.test(.x, pt_tidy$group, p.method = "none"))
+
+pt_variance_tidy <- pt_variance %>%
+ map(tidy) %>%
+ map2(.y = names(pt_variance), .f = ~.x %>% mutate(pt_vars = .y)) %>%
+ purrr::reduce(bind_rows) %>%
+ separate(group1,
+ into = c("genotype1", "irrigation1"),
+ sep = "_") %>%
+ separate(group2,
+ into = c("genotype2", "irrigation2"),
+ sep = "_") %>%
+ filter(irrigation1 == irrigation2, genotype1 == "MoneyMaker"& genotype2 != "MoneyMaker"|
+ genotype1 != "MoneyMaker" & genotype2 == "MoneyMaker") %>%
+ mutate(genotype = if_else(genotype1 == "MoneyMaker", genotype2, genotype1),
+ comp = "MoneyMaker") %>%
+ select(genotype, comp, irrigation = irrigation1, pt_var = pt_vars, p.value)
+
+prereq <- pt_norm_comb %>%
+ left_join(pt_variance_tidy) %>%
+ mutate(fulfill = if_else(wt_norm >= 0.05, norm >= 0.05, p.value >= 0.05))
+
+prereq_sum <- prereq %>%
+ group_by(fulfill) %>%
+ summarise(n = n()) %>%
+ ungroup() %>%
+ mutate(percent = 100*n/sum(n))
+
+pt_aov <- map(.x = pt_tidy_numeric, .f = ~aov(.x ~ alias* irrigation, data = pt_tidy))
+
+pt_aov_tidy <- pt_aov %>%
+ map(tidy) %>%
+ map2(.y = names(pt_aov), .f = ~.x %>% mutate(pt_vars = .y)) %>%
+ purrr::reduce(bind_rows)
+
+pt_aov_table <- pt_aov_tidy %>%
+ left_join(pt_vars_plot) %>%
+ filter(term != "Residuals") %>%
+ select(term, pt_vars_label, p.value) %>%
+ pivot_wider(names_from = term, values_from = p.value)
+
+write.xlsx(pt_aov_table, str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "anova_table.xlsx",
+ sep = "_"))
+
+pt_tuk <- map(.x = pt_aov, .f = ~TukeyHSD(.x)) %>%
+ map(.f = tidy) %>%
+ map2(.y = names(pt_aov), .f = ~.x %>% mutate(var = .y)) %>%
+ purrr::reduce(bind_rows)
+
+# General differences -----------------------------------------------------
+
+sig_pt_overall <- pt_tuk %>%
+ filter(term == "alias") %>%
+ separate(col = contrast, into = c("group1", "group2"), sep = "-")
+
+lines <- sum_stat %>% distinct(alias,plantline) %>%
+ mutate(plantline = as.double(as.character(plantline))) %>%
+ left_join(genotypes)
+
+sig_pt_overall_table <- sig_pt_overall %>%
+ filter(group1 == "967514 MM WT" | group2 == "967514 MM WT") %>%
+ ungroup() %>%
+ left_join(lines, by = c("group1" = "alias")) %>%
+ select(var, group1, group2, adj.p.value, genotype1 = genotype) %>%
+ left_join(lines, by = c("group2" = "alias")) %>%
+ select(var,adj.p.value, genotype1, genotype2 = genotype) %>%
+ left_join(pt_vars_plot, by = c("var" = "pt_vars")) %>%
+ pivot_wider(names_from = genotype2,
+ values_from = adj.p.value)
+
+write.xlsx(sig_pt_overall_table, str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "overall_significance.xlsx",
+ sep = "_"))
+
+# Pairwise significant changes --------------------------------------------
+
+per_comp_y <- pt_long %>%
+ group_by(pt_var, irrigation) %>%
+ summarise(max_y = max(value),
+ min_y = min(value))
+
+sig_pt_groups <- pt_tuk %>%
+ filter(term == "alias:irrigation") %>%
+ separate(col = contrast, into = c("group1", "group2"), sep = "-")
+
+sig_pt <- sig_pt_groups %>%
+ separate(group1, into = c("alias1", "irrigation1"), sep = ":") %>%
+ separate(group2, into = c("alias2", "irrigation2"), sep = ":") %>%
+ filter(irrigation1 == irrigation2) %>%
+ filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT") %>%
+ ungroup() %>%
+ mutate(p.signif = if_else(adj.p.value <= 0.05, "*", "ns")) %>%
+ left_join(sum_stat, by = c("var" = "pt_var", "irrigation1" = "irrigation", "alias1" = "alias")) %>%
+ select(var, irrigation1, irrigation2, alias1, alias2, p.signif, mean1 = mean, se1 = se, genotype1 = genotype, adj.p.value) %>%
+ left_join(sum_stat, by = c("var" = "pt_var", "irrigation2" = "irrigation", "alias2" = "alias")) %>%
+ select(var, irrigation1, irrigation2, alias1, alias2, p.signif, mean1, mean2 = mean, se1, se2 = se, genotype1, genotype2 = genotype, adj.p.value) %>%
+ rename(pt_var = var,
+ group1 = genotype1,
+ group2 = genotype2,
+ irrigation = irrigation1) %>%
+ left_join(per_comp_y) %>%
+ mutate(y.position = 1.1*(max_y)) %>%
+ mutate(pt_var = as_factor(pt_var),
+ group1 = as_factor(group1),
+ group2 = as_factor(group2),
+ irrigation = as_factor(irrigation)) %>%
+ ungroup()
+
+genotypes_label <- sum_stat %>%
+ distinct(plantline, genotype)
+
+pt_vars
+
+plot_out <- vector("list", length = length(pt_vars))
+
+for(pt in seq_along(pt_vars)) {
+
+ binwidth <- per_comp_y %>%
+ filter(pt_var == pt_vars[[pt]]) %>%
+ summarise(min_y_comp = min(min_y),
+ max_y_comp = max(max_y)) %>%
+ ungroup() %>%
+ mutate(binwidth = (max_y_comp - min_y_comp)/50) %>%
+ select(binwidth) %>%
+ as_vector() %>%
+ as.numeric()
+
+ sig_pt_temp <- sig_pt %>%
+ filter(pt_var == pt_vars[[pt]])
+
+ plot_out[[pt]] <- pt_long %>%
+ select(-genotype) %>%
+ left_join(genotypes_label) %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ filter(pt_var == pt_vars[[pt]]) %>%
+ ggplot(aes(genotype, value)) +
+ geom_boxplot(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_dotplot(aes(fill = genotype), color = "black", binaxis = "y", stackdir = "center", binwidth = binwidth) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ hide.ns = T) +
+ theme(axis.text.x = element_markdown(angle = 45, hjust = 1, margin = unit(c(1,0,0,0), "mm")),
+ axis.title.x = element_blank(),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ text = element_text(size = 8, family = "Arial Narrow"),
+ legend.title = element_blank(),
+ legend.text = element_markdown(),
+ legend.position = "right",
+ legend.direction = "vertical",
+ plot.margin = unit(c(1,0,1,2), "mm"),
+ legend.margin = margin(t = 0, r = 2, b = 0, l = 2 , unit = "mm")) +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_color_manual(values = cb_scale, aesthetics = "color") +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ xlab(label = "") +
+ ggtitle(pt_vars_label[[pt]])
+}
+
+plot_out[[pt]]
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_phenotypes.pdf"),
+ width = 6.5,
+ height = 3)
+
+for (pt in seq_along(pt_vars)) {
+ print(plot_out[[pt]])
+}
+
+dev.off()
+
+if (!is.null(dev.list())) {
+ dev.off()
+}
+
+sum_stat %>%
+ ggplot(aes(irrigation, mean, color = alias)) +
+ geom_point() +
+ geom_line(aes(group = plantline)) +
+ geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd)) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+
+
+sum_stat %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ ggplot(aes(irrigation, mean, color = genotype)) +
+ geom_line(aes(group = genotype)) +
+ geom_pointrange(aes(ymin = mean - se, ymax = mean + se, shape = genotype, fill = genotype)) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ 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 = 8),
+ legend.text = element_markdown(),
+ legend.position = "right") +
+ scale_shape_manual(values = c(1,2,22,23))
+
+sum_stat %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ ggplot(aes(irrigation, mean, color = genotype)) +
+ geom_line(aes(group = genotype)) +
+ geom_point(aes(shape = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se, color = genotype, width = 0.25)) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ 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 = 8),
+ legend.text = element_markdown(),
+ legend.position = "right") +
+ scale_shape_manual(values = c(1,2,22,23))
+
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "cmQTL_val1_phenotypic.jpg",
+ sep = "_"),
+ width = 183,
+ height = 100,
+ units = "mm",
+ dpi = 300)
+
+
+# Overview for supplemental -----------------------------------------------
+
+bins <- pt_long %>%
+ group_by(pt_var) %>%
+ summarise(binwidth = max(value)/60)
+
+cb_scale <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442",
+ "#0072B2", "#D55E00","#000000")
+
+pt_long %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ group_by(plantline, genotype, irrigation, pt_var) %>%
+ ggplot(aes(irrigation, value, fill = genotype)) +
+ geom_dotplot(binaxis = "y", stackdir = "center") +
+ geom_smooth(aes(group = genotype, color = genotype), method = "lm", formula = y ~ poly(x,2), se = F, size = 1) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ 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 = 8, family = "Arial Narrow"),
+ legend.text = element_markdown(),
+ legend.position = "right") +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ scale_color_manual(values = cb_scale, aesthetics = "color")
+
+sum_stat %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ ggplot(aes(irrigation, mean, color = genotype)) +
+ geom_line(aes(group = genotype), size = 0.3) +
+ geom_point(shape = 1, size = 2) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se, color = genotype, width = 0.25), size = 0.3) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ 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 = 8, family = "Arial Narrow"),
+ legend.text = element_markdown(),
+ legend.position = "right") +
+ #scale_shape_manual(values = c(1,2,22,23)) +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ scale_color_manual(values = cb_scale, aesthetics = "color")
+
+ggsave("Overview_phenotype_dot_line.svg",width = 17.0, height = 8, units = "cm", dpi = 300)
+
+overview <- pt_long %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+ genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*"))) %>%
+ group_by(plantline, genotype, irrigation, pt_var, pt_vars_label) %>%
+ summarise(mean = mean(value)) %>%
+ ggplot(aes(irrigation, mean, fill = genotype)) +
+ #geom_dotplot(binaxis = "y", stackdir = "center") +
+ geom_line(aes(group = genotype, color = genotype), size = 0.5) +
+ geom_point(aes(group = genotype, color = genotype), size = 0.5) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ theme(axis.text.x = element_markdown(angle = 45, hjust = 1, size = 8, family = "sans"),
+ axis.title.x = element_text(size = 8, family = "sans"),
+ axis.title.y = element_text(size = 8, family = "sans"),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ strip.text = element_text(size = 6, family = "sans"),
+ text = element_text(size = 8, family = "sans"),
+ legend.title = element_text(size = 10, family = "sans"),
+ legend.text = element_markdown(size = 8),
+ legend.position = "right") +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ scale_color_manual(values = cb_scale, aesthetics = "color"); overview
+
+ggsave("Overview_phenotype_line_mean.png", width = 17.0, height = 8, units = "cm", dpi = 300, scale = 1)
+
+pt_long %>%
+ #filter(genotype == "MoneyMaker") %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ group_by(plantline, genotype, irrigation, pt_var) %>%
+ ggplot(aes(irrigation, value)) +
+ geom_boxplot(fill = "orange", size = 0.25) +
+ geom_smooth(aes(group = pt_var), method = "lm", formula = y ~ poly(x,2), se = F, size = 1) +
+ facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ com_theme +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ scale_color_manual(values = cb_scale, aesthetics = "color")
+
+ggsave("Overview_phenotype_box_line_all_g_aggregated.png",width = 17.0, height = 8, units = "cm", dpi = 300)
+saveRDS(last_plot(), "Overview_phenotype_box_line_all_g_aggregated.RDS")
+
+sum_stat_out <- sum_stat %>%
+ mutate(genotype = str_remove_all(genotype, "\\*")) %>%
+ left_join(pt_norm) %>%
+ left_join(pt_variance_tidy) %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ left_join(sig_pt, by = c("irrigation", "pt_var", "alias" = "alias2")) %>%
+ mutate(genotype = as_factor(genotype),
+ genotype = fct_relevel(genotype, c("MoneyMaker", "panK4-1", "log2-1", "transp1-1"))) %>%
+ arrange(pt_vars_label, irrigation, genotype) %>%
+ select(Genotype = genotype, Irrigation = irrigation, Trait = pt_vars_label,
+ Mean = mean, SD = sd, SE = se, N = n, Normality = norm, `Equal variance` = p.value,
+ `P-value` = adj.p.value)
+
+write.xlsx(sum_stat_out, "Summary_statistics.xlsx")
+
+pt_raw_out <- pt_long %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ arrange(pt_vars_label, irrigation, genotype) %>%
+ select(Genotype = genotype, Irrigation = irrigation, Replicate = rep, Trait = pt_vars_label,
+ Value = value)
+
+write.xlsx(pt_raw_out, "Raw_data_phenotypes.xlsx")
+
+
+# Figure part panK4-1 total fruit number ----------------------------------
+
+A <- sum_stat %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ filter(pt_var == "total_fruit_num", genotype == "MoneyMaker" | genotype == "*panK4-1*") %>%
+ ggplot(aes(irrigation, mean, color = genotype)) +
+ geom_line(aes(group = genotype)) +
+ geom_point(aes(color = genotype),
+ size = 1) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se, color = genotype),
+ width = 0.1) +
+ geom_text(aes(x = "0.8", y = 12, label = "TukeyHSD: 1.30*10^-3"), inherit.aes = F, size = 2) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ #facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ com_theme +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_color_manual(values = cb_scale, aesthetics = "color")
+
+A
+
+ggsave("Total_fruit_number_panK4.png",
+ width = 8.25,
+ height = 3,
+ units = "cm",
+ dpi = 300)
+
+saveRDS(A, "Total_fruit_number_panK4.RDS")
+
+# Figure part log2-1 average fruit weight ----------------------------------
+
+
+B <- sum_stat %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ filter(pt_var == "av_fruit_weight", genotype == "MoneyMaker" | genotype == "*log2-1*") %>%
+ ggplot(aes(irrigation, mean, color = genotype)) +
+ geom_line(aes(group = genotype)) +
+ geom_point(aes(color = genotype),
+ size = 1) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se, color = genotype),
+ width = 0.1) +
+ geom_text(aes(x = "0.6", y = 17, label = "TukeyHSD: 6.69*10^-6"), inherit.aes = F, size = 2) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ #facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ com_theme +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_color_manual(values = cb_scale[c(1,3)], aesthetics = "color")
+
+B
+
+ggsave("Average_fruit_weight_log2-1.png",
+ width = 8.25,
+ height = 3,
+ units = "cm",
+ dpi = 300)
+
+saveRDS(B, "Average_fruit_weight_log2-1.RDS")
+
+# Figure part transp1-1 days until red ripe ----------------------------------
+
+
+C <- sum_stat %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ filter(pt_var == "days_until_redripe", genotype == "MoneyMaker" | genotype == "*transp1-1*") %>%
+ ggplot(aes(irrigation, mean, color = genotype)) +
+ geom_line(aes(group = genotype)) +
+ geom_point(aes(color = genotype),
+ size = 1) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se, color = genotype),
+ width = 0.1) +
+ geom_text(aes(x = "0.6", y = 110, label = "TukeyHSD: 3.80*10^-3"), inherit.aes = F, size = 2) +
+ #geom_text(aes(label = is_sig), nudge_y = 10) +
+ #facet_wrap(nrow = 2, ncol = 4, facets = vars(pt_vars_label),scales = "free")+
+ com_theme +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_color_manual(values = cb_scale[c(1,4)], aesthetics = "color")
+
+C
+
+ggsave("Days_until_red_ripe_transp1-1.png",
+ width = 8.25,
+ height = 3,
+ units = "cm",
+ dpi = 300)
+
+saveRDS(C, "Days_until_red_ripe_transp1-1.RDS")
+# Combining plots ---------------------------------------------------------
+
+
+main_plot <- plot_grid(A,B,C, nrow = 1, labels = "AUTO", label_size = 12)
+
+main_plot
+
+ggsave("Selected_phenotypes.png",
+ width = 17.0,
+ height = 6,
+ units = "cm",
+ dpi = 300)
+
+
+sum_stat %>%
+ ggplot(aes(irrigation, mean, color = alias)) +
+ geom_point() +
+ geom_line(aes(group = plantline)) +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+ggsave(str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+ "cmQTL_val1_phenotypic.jpg",
+ sep = "_"),
+ width = 183,
+ height = 100,
+ units = "mm",
+ dpi = 300)
+
+# CV analysis -------------------------------------------------------------
+
+
+
+library(bootstrap)
+
+theta <- function(x){
+ sd(x)/mean(x)
+}
+
+CV <- sum_stat %>%
+ ungroup() %>%
+ group_by(plantline, genotype, alias, pt_var) %>%
+ summarise(grand_mean = mean(mean),
+ grand_sd = sd(mean)) %>%
+ mutate(cv = grand_sd/grand_mean)
+
+cv_jack <- sum_stat %>%
+ group_by(plantline, genotype, alias, pt_var) %>%
+ summarise(cv = jackknife(mean, theta)$jack.values) %>%
+ mutate(jack_rep = row_number()) %>%
+ ungroup()
+
+cv_jack_wide <- cv_jack %>%
+ pivot_wider(id_cols = c(plantline, genotype, alias, jack_rep),
+ names_from = pt_var,
+ values_from = cv)
+
+cv_jack_numeric <- cv_jack_wide %>%
+ select(-c(alias, genotype, plantline, jack_rep))
+
+cv_jack_sum_stat <- cv_jack %>%
+ group_by(plantline, genotype, alias, pt_var) %>%
+ summarise(mean_cv = mean(cv),
+ sd = sd(cv),
+ n = n(),
+ se = sd/sqrt(n))
+
+
+pt_cv_jack <- map(.x = cv_jack_numeric, .f = ~aov(.x ~ alias, data = cv_jack_wide))
+
+pt_jack_t <- map(.x = cv_jack_numeric, .f = ~pairwise.t.test(x = .x , g = cv_jack_wide$alias)) %>%
+ map(.f = tidy) #%>%
+
+pt_jack_t_tidy <- pt_jack_t %>%
+ map2(.y = names(pt_jack_t), .f = ~.x %>% mutate(var = .y)) %>%
+ purrr::reduce(bind_rows)
+
+sig_pt_cv_groups <- pt_jack_t_tidy %>%
+ filter(group1 == "967514 MM WT" | group2 == "967514 MM WT") %>%
+ 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_sum_stat, by = c("var" = "pt_var", "group1" = "alias")) %>%
+ select(var, group1, group2, p.signif, mean_cv1 = mean_cv, se1 = se, genotype1 = genotype) %>%
+ left_join(cv_jack_sum_stat, by = c("var" = "pt_var", "group2" = "alias")) %>%
+ select(var, group1, group2, p.signif, mean_cv1, mean_cv2 = mean_cv, se1, se2 = se, genotype1, genotype2 = genotype) %>%
+ group_by(var) %>%
+ mutate(tot_val1 = mean_cv1 + se1,
+ tot_val2 = mean_cv2 + se2,
+ y.position = 1.1*(max(tot_val1, tot_val2))) %>%
+ rename(pt_var = var) %>%
+ mutate(pt_var = as_factor(pt_var),
+ group1 = genotype1,
+ group2 = genotype2) %>%
+ ungroup()
+
+cv_jack %>%
+ ggplot(aes(alias, cv, fill = alias)) +
+ geom_boxplot() +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+
+
+
+CV %>%
+ ggplot(aes(alias, cv, fill = alias)) +
+ geom_col() +
+ facet_wrap(facets = vars(pt_var),scales = "free")
+
+x = rnorm(10)
+
+# Figure Average fruit weight ---------------------------------------------
+
+
+## Nominal all lines -----------------------------------------------------
+
+binwidth <- pt_long %>%
+ filter(pt_var == "av_fruit_weight",
+ genotype %in% c("MoneyMaker", "*panK4-1*")) %>%
+ group_by(pt_var) %>%
+ summarise(min_y_comp = min(value),
+ max_y_comp = max(value)) %>%
+ ungroup() %>%
+ mutate(binwidth = (max_y_comp - min_y_comp)/50) %>%
+ select(binwidth) %>%
+ as_vector() %>%
+ as.numeric()
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "av_fruit_weight",
+ group1 %in% c("MoneyMaker", "*panK4-1*"),
+ group2 %in% c("MoneyMaker", "*panK4-1*"))
+
+A <- pt_long %>%
+ select(-genotype) %>%
+ left_join(genotypes_label) %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ filter(pt_var == "av_fruit_weight", genotype %in% c("MoneyMaker", "*panK4-1*")) %>%
+ ggplot(aes(genotype, value)) +
+ geom_boxplot(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_dotplot(aes(fill = genotype), color = "black", binaxis = "y", stackdir = "center", binwidth = binwidth, dotsize = 1.5) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ hide.ns = T) +
+ theme(axis.text.x = element_markdown(angle = 45, hjust = 1, size = 8, family = "sans"),
+ axis.title.x = element_text(size = 8, family = "sans"),
+ axis.title.y = element_text(size = 8, family = "sans"),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ strip.text = element_text(size = 6, family = "sans"),
+ text = element_text(size = 8, family = "sans"),
+ legend.title = element_text(size = 10, family = "sans"),
+ legend.text = element_markdown(size = 8),
+ legend.position = "bottom",
+ legend.direction = "horizontal",
+ plot.margin = unit(c(1,0,1,2), "mm"),
+ legend.margin = margin(t = 0, r = 2, b = 0, l = 2 , unit = "mm")) +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_color_manual(values = cb_scale, aesthetics = "color") +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ xlab(label = "")
+
+A
+
+leg <- get_legend(A)
+
+A <- A + theme(legend.position = "")
+A
+
+## CV all lines ------------------------------------------------------------
+
+per_comp_y_cv <- cv_jack %>%
+ group_by(pt_var) %>%
+ summarise(max_y = max(cv),
+ min_y = min(cv))
+
+binwidth <- per_comp_y_cv %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ summarise(min_y_comp = min(min_y),
+ max_y_comp = max(max_y)) %>%
+ ungroup() %>%
+ mutate(binwidth = (max_y_comp - min_y_comp)/50) %>%
+ select(binwidth) %>%
+ as_vector() %>%
+ as.numeric()
+
+sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == "av_fruit_weight")
+
+B <- cv_jack %>%
+ select(-genotype) %>%
+ left_join(genotypes_label) %>%
+ left_join(pt_vars_plot, by = c("pt_var" = "pt_vars")) %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, cv)) +
+ geom_boxplot(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_dotplot(aes(fill = genotype), color = "black", binaxis = "y", stackdir = "center", binwidth = binwidth, dotsize = 1.5) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ hide.ns = T) +
+ theme(axis.text.x = element_markdown(angle = 45, hjust = 1, size = 8, family = "sans"),
+ axis.title.x = element_text(size = 8, family = "sans"),
+ axis.title.y = element_text(size = 8, family = "sans"),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ strip.text = element_text(size = 6, family = "sans"),
+ text = element_text(size = 8, family = "sans"),
+ legend.title = element_text(size = 10, family = "sans"),
+ legend.text = element_markdown(size = 8),
+ legend.position = "",
+ legend.direction = "vertical",
+ plot.margin = unit(c(1,1,1,2), "mm"),
+ legend.margin = margin(t = 0, r = 2, b = 0, l = 2 , unit = "mm")) +
+ scale_shape_manual(values = c(1,2,22,23)) +
+ scale_color_manual(values = cb_scale, aesthetics = "color") +
+ scale_fill_manual(values = cb_scale, aesthetics = "fill") +
+ xlab(label = "")
+
+B
+
+
+# Combining plots ---------------------------------------------------------
+
+
+main_plot <- plot_grid(A,B, nrow = 1, labels = "AUTO", rel_widths = c(2,1), label_size = 12)
+main_plot
+
+comp_plot <- plot_grid(main_plot, ncol = 1, leg, rel_heights = c(10,1))
+comp_plot
+
+ggsave("Average_fruit_weight_all_lines.png",
+ width = 17.0,
+ height = 8,
+ units = "cm",
+ dpi = 300)
+
+# Jackknife plots---------------------------------------------------------------
+
+
+
+
+pt_vars
+pt_vars_label <- c("Total fruit number per plant", "Total fruit weight per plant",
+ "Average fruit weight", "Days until ripe fruit", "Shoot fresh weight",
+ "Harvest index", "Biomass", "Shoot dry weight")
+
+plot_out <- vector("list", length = length(pt_vars))
+
+for(pt in seq_along(pt_vars)) {
+ sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == pt_vars[[pt]])
+
+ plot_out[[pt]] <- cv_jack_sum_stat %>%
+ filter(pt_var == pt_vars[[pt]]) %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.75) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "right") +
+ ggtitle(pt_vars_label[[pt]]) +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV")
+}
+
+plot_out[[1]]
+
+pdf(file = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),"cmQTL_val1_phenotypes_cv.pdf"),
+ width = 9,
+ height = 4)
+
+
+
+
+
+
+for (pt in seq_along(pt_vars)) {
+ print(plot_out[[pt]])
+}
+
+dev.off()
+
+
+
+plot(bootstrap::jackknife(x = x, theta = theta)$jack.values)
+
+stop(print("current stop"))
+
+# Figure 40 ---------------------------------------------------------------
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "total_fruit_weight")
+
+p1<- sum_stat %>%
+ filter(pt_var == "total_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "",
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ xlab("") +
+ ylim(0,1150) +
+ ggtitle("Total fruit weight per plant") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+p1
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "av_fruit_weight")
+
+p2<- sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "",
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ xlab("") +
+ ylim(c(0,65))+
+ ggtitle("Average fruit weight") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+p2
+
+leg40 <- get_legend(sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "right") +
+ ggtitle("Average fruit weight") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill"))
+leg40
+
+main_plot <- plot_grid(p1,p2, ncol = 1, labels = "AUTO")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg40, ncol = 2, rel_widths = c(3,1))
+comp_plot
+
+ggsave("figure_40.wmf",width = 15.8, height = 12, units = "cm", dpi = 300)
+
+# Figure 41 ---------------------------------------------------------------
+
+
+sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == "av_fruit_weight")
+
+p1 <- cv_jack_sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "",
+ plot.title = element_text(size = 14),
+ plot.margin = unit(c(0,0,1,0), "mm")) +
+ ggtitle("Average fruit weight") +
+ xlab("") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV")
+p1
+
+sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == "harvest_index")
+
+p2 <- cv_jack_sum_stat %>%
+ filter(pt_var == "harvest_index") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "",
+ plot.title = element_text(size = 14),
+ plot.margin = unit(c(0,0,1,0), "mm")) +
+ ggtitle("Harvest index") +
+ xlab("") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV")
+p2
+
+leg41 <- get_legend(cv_jack_sum_stat %>%
+ filter(pt_var == "harvest_index") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "bottom") +
+ ggtitle("Harvest index") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV"))
+
+
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "AUTO")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg41, ncol = 1, rel_heights = c(10,1)) + theme(plot.margin = unit(c(0,0,5,0), "mm"))
+comp_plot
+
+ggsave("figure_41.wmf",width = 15.8, height = 8, units = "cm", dpi = 300)
+
+
+# Figure means presentation -----------------------------------------------
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "total_fruit_weight")
+
+p1<- sum_stat %>%
+ filter(pt_var == "total_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "",
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ xlab("") +
+ ylim(0,1150) +
+ ggtitle("Total fruit weight per plant") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+p1
+
+sig_pt_temp <- sig_pt %>%
+ filter(pt_var == "av_fruit_weight")
+
+p2<- sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "",
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ xlab("") +
+ ylim(c(0,65))+
+ ggtitle("Average fruit weight") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+p2
+
+leg40 <- get_legend(sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean - se, ymax = mean + se), position = position_dodge(width = 0.9), width = 0.25) +
+ facet_grid(cols = vars(irrigation), scales = "fixed") +
+ stat_pvalue_manual(sig_pt_temp, 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(),
+ legend.position = "right") +
+ ggtitle("Average fruit weight") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill"))
+leg40
+
+main_plot <- plot_grid(p1,p2, ncol = 2, labels = "")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg40, ncol = 2, rel_widths = c(6,1))
+comp_plot
+
+ggsave("Av_tot_fruit_weight_mean.wmf",width = 33.5, height = 12, units = "cm", dpi = 300)
+
+# Figure CV average fruit weight ---------------------------------------------------------------
+
+
+sig_pt_temp <- sig_pt_cv_groups %>%
+ filter(pt_var == "av_fruit_weight")
+
+p1 <- cv_jack_sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "",
+ plot.title = element_text(size = 14),
+ plot.margin = unit(c(0,0,0,0), "mm")) +
+ ggtitle("Average fruit weight") +
+ xlab("") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV")
+p1
+
+
+
+leg41 <- get_legend(cv_jack_sum_stat %>%
+ filter(pt_var == "av_fruit_weight") %>%
+ ggplot(aes(genotype, mean_cv)) +
+ geom_col(position = "dodge", color = "black", aes(fill = genotype)) +
+ geom_errorbar(aes(ymin = mean_cv - se, ymax = mean_cv + se), position = position_dodge(width = 0.9), width = 0.25) +
+ stat_pvalue_manual(sig_pt_temp, label = "p.signif", y.position = "y.position",
+ step.increase = 0.1,
+ 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(),
+ legend.position = "right") +
+ ggtitle("Harvest index") +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill") +
+ ylab("mean CV"))
+
+
+
+main_plot <- plot_grid(p1, ncol = 1, labels = "")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg41, ncol = 2, rel_widths = c(3,2)) + theme(plot.margin = unit(c(0,0,0,0), "mm"))
+comp_plot
+
+ggsave("av_fruit_weight_cv_presentation.wmf",width = 10, height = 10, units = "cm", dpi = 300)
+
+
+
+
+setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
+getwd()
+
+library(tidyverse)
+library(ggpubr)
+library(glue)
+library(ggtext)
+library(cowplot)
+
+rm(list = ls())
+fc_1 <- read_csv("211012_fc_1.csv") %>%
+ mutate(genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+sig_lip <- read_csv("211012_sig_lip.csv")
+
+fc_cv <- read.csv("211012_fc_cv.csv")%>%
+ mutate(genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+sig_lip_cv <- read_csv("211012_sig_lip_cv_groups.csv")
+
+met_dat <- fc_1 %>%
+ distinct(met, Compound_Name, Compound_Class)
+
+genotypes <- fc_1 %>%
+ distinct(alias, genotype)
+
+per_comp_y <- fc_1 %>%
+ group_by(tissue, treatment, met) %>%
+ summarise(tot_val = max(mean_fc + se))
+
+make_col_plot <- function(plot_met, plot_tissue, plot_label, plot_legend){
+
+ sig_bar <- sig_lip %>%
+ left_join(genotypes, by = c("group1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("group2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y) %>%
+ mutate(y.position = 1.1 * tot_val) %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ plot_out <- fc_1 %>%
+ filter(met == plot_met, tissue == plot_tissue) %>%
+ 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 = 10),
+ legend.text = element_markdown(),
+ legend.position = plot_legend,
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ ylab("Mean fold-change")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+ plot_out
+
+}
+
+per_comp_y_cv <- fc_cv %>%
+ group_by(tissue, met) %>%
+ summarise(tot_val = max(mean_fc + se))
+
+make_col_plot_cv <- function(plot_met, plot_tissue, plot_label, plot_legend){
+
+ sig_bar <- sig_lip_cv %>%
+ left_join(genotypes, by = c("alias1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("alias2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y_cv) %>%
+ mutate(y.position = 1.1 * tot_val) %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ #ylim_top <- 1.1*max(sig_bar$y.position)
+
+ plot_out <- fc_cv %>%
+ filter(met == plot_met, tissue == plot_tissue) %>%
+ 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 = 10),
+ legend.text = element_markdown(),
+ legend.position = plot_legend) +
+ ylab("CV fold-change")+
+ #ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+ plot_out
+
+}
+
+#make_col_plot("m_1", "fruits", "test")
+
+make_col_plot_no_x <- function(plot_met, plot_tissue, plot_label, plot_legend){
+
+ sig_bar <- sig_lip %>%
+ left_join(genotypes, by = c("group1" = "alias")) %>%
+ rename(genotype1 = genotype) %>%
+ left_join(genotypes, by = c("group2" = "alias")) %>%
+ rename(genotype2 = genotype) %>%
+ select(-group1, -group2) %>%
+ #mutate(genotype1 = as_factor(if_else(genotype1 == "MoneyMaker", glue("{genotype1}"),glue("*{genotype1}*"))),
+ # genotype2 = as_factor(if_else(genotype2 == "MoneyMaker", glue("{genotype2}"),glue("*{genotype2}*")))) %>%
+ rename(group1 = genotype1, group2 = genotype2) %>%
+ left_join(per_comp_y) %>%
+ mutate(y.position = 1.1 * tot_val) %>%
+ filter(met == plot_met, tissue == plot_tissue)
+
+ ylim_top <- 1.1*max(sig_bar$y.position)
+
+ plot_out <- fc_1 %>%
+ filter(met == plot_met, tissue == plot_tissue) %>%
+ 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_blank(),
+ panel.background = element_rect(fill = "white"),
+ panel.border = element_rect(color = "black",fill = NA),
+ text = element_text(size = 10),
+ legend.text = element_markdown(),
+ legend.position = plot_legend,
+ plot.margin = unit(c(0,0,0,0), "cm")) +
+ ylab("Mean fold-change")+
+ ylim(c(0, ylim_top))+
+ xlab("") +
+ ggtitle(label = plot_label) +
+ scale_fill_grey(start = 1, end = 0, aesthetics = "fill")
+
+ plot_out
+
+}
+
+# Figure 58 ---------------------------------------------------------------
+
+p1 <- make_col_plot_no_x("Cluster_08623", "leaves", "TAG 48:1", "none")
+p1
+
+p2 <- make_col_plot_no_x("Cluster_08498", "leaves", "TAG 48:3 (1)", "none")
+p2
+
+p3 <- make_col_plot_no_x("Cluster_08438", "leaves", "TAG 48:4", "none")
+p3
+
+p4 <- make_col_plot_no_x("Cluster_09298", "leaves", "TAG 50:3 (1)", "none")
+p4
+
+p5 <- make_col_plot("Cluster_09251", "leaves", "TAG 50:4 (1)", "none")
+p5
+
+p6 <- make_col_plot("Cluster_10916", "leaves", "TAG 54:5", "none")
+p6
+
+leg58 <- get_legend(make_col_plot("Cluster_08623", "leaves", "sucrose", "bottom"))
+leg58
+
+main_plot <- plot_grid(p1,p2,p3,p4,p5,p6, ncol = 2, rel_heights = c(1,1,1.35), labels = "AUTO")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg58, ncol = 1, rel_heights = c(10,1))
+comp_plot
+
+ggsave("figure_58.wmf",width = 15.8, height = 16, units = "cm", dpi = 300)
+
+# Figure 61 ---------------------------------------------------------------
+
+p1 <- make_col_plot_cv("Cluster_09079", "leaves", "TAG 50:7 (2)", "none")
+p1
+
+p2 <- make_col_plot_cv("Cluster_10388", "leaves", "TAG 52:0", "none")
+p2
+
+p3 <- make_col_plot_cv("Cluster_09920", "leaves", "TAG 52:7 (2)", "none")
+p3
+
+p4 <- make_col_plot_cv("Cluster_11255", "leaves", "TAG 54:0", "none")
+p4
+
+leg61 <- get_legend(make_col_plot("Cluster_09079", "leaves", "sucrose", "bottom"))
+leg61
+
+main_plot <- plot_grid(p1,p2,p3,p4, ncol = 4, labels = "AUTO")
+main_plot
+
+comp_plot <- plot_grid(main_plot, leg61, ncol = 1, rel_heights = c(10,1))
+comp_plot
+
+ggsave("figure_61.wmf",width = 15.8, height = 8, units = "cm", dpi = 300)