diff --git a/.Rhistory b/.Rhistory
index 97ff8ae14324425df3fa0294343461689c85ea38..23867cfd3e2ace9be782c58cb62d9c5057e795e6 100644
--- a/.Rhistory
+++ b/.Rhistory
@@ -1,512 +1,512 @@
-group = vector(mode = "list", length = length(x))
-for(i in seq_along(x)){
-if(str_detect(x[[i]], "Characteristic|Parameter|Factor|Component")){
-group[[i]] <- tibble(id = x[[i]])
-if(str_detect(x[[i+1]], "Unit")) {
-group[[i]] <- group[[i]] %>% mutate(isa_unit = x[[i+1]])
-if(str_detect(x[[i+2]], "Term Source REF")) {
-group[[i]] <- group[[i]] %>% mutate(term_source = x[[i+2]])
-if(str_detect(x[[i+3]], "Term Accession Number")){
-group[[i]] <- group[[i]] %>% mutate(term_acc = x[[i+3]])
}
}
-} else if(str_detect(x[[i+1]], "Term Source REF")) {
-group[[i]] <- group[[i]] %>% mutate(term_source = x[[i+1]])
-if(str_detect(x[[i+2]], "Term Accession Number")){
-group[[i]] <- group[[i]] %>% mutate(term_acc = x[[i+2]])
-}
-}
-} else if(str_detect(x[[i]], "Source REF|Accession Number|Unit",negate = T)){
-group[[i]] <- tibble(id = x[[i]])
-}
-}
-group_tib <- group %>%
-discard(.p = is.null) %>%
-reduce(.f = bind_rows)
-assign(output_name,group_tib, pos = 1)
-}
-break_lines(sam_dat1, "output_samdat")
-break_lines(isa_ext, "output_ext")
-break_lines(isa_gc, "output_gc")
-View(output_gc)
-output <- bind_rows(output_samdat,
-output_ext,
-output_gc,
-output_ms)
-beak_lines(isa_ms, "output_ms")
-break_lines(isa_ms, "output_ms")
-output <- bind_rows(output_samdat,
-output_ext,
-output_gc,
-output_ms)
-here()
-getwd()
-setwd(here())# Not recommended but convenient in Rstudio to start from root
-write.xlsx(output,
-here("assays/cmQTL_val1_GH_2020_GC_MS/dataset/isa.dataset_temp.xlsx"))
-here(out,
-"cmQTL_val_1_2_met_dat_GC.csv"
+dev.off()
+heat.GC <- fc %>%
+group_by(Compound_Name, met) %>%
+mutate(log_norm = log2(mean_fc),
+log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>%
+ungroup() %>%
+#mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+pivot_wider(id_cols = c(Compound_Name, met),
+names_from = group,
+values_from = log_norm) %>%
+left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+as.data.frame()
+norm_MM <- means %>%
+filter(genotype == "MoneyMaker") %>%
+select(tissue, treatment, met, MM_mean = mean)
+fc <- means_tec_rep %>%
+left_join(norm_MM) %>%
+mutate(fc = mean_tec_rep/MM_mean) %>%
+group_by(tissue, treatment, alias, genotype, met) %>%
+summarise(mean_fc = mean(fc),
+sd = sd(fc),
+n = n()) %>%
+mutate(se = sd/sqrt(n),
+group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+left_join(sig_GC, by = c("tissue", "treatment", "met", "alias" = "group2")) %>%
+group_by(met) %>%
+filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>%
+left_join(met_dat) %>%
+arrange(Compound_Class, Compound_Name) %>%
+ungroup() %>%
+mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+fc <- means_tec_rep %>%
+left_join(norm_MM) %>%
+mutate(fc = mean_tec_rep/MM_mean) %>%
+group_by(tissue, treatment, alias, genotype, met) %>%
+summarise(mean_fc = mean(fc),
+sd = sd(fc),
+n = n()) %>%
+mutate(se = sd/sqrt(n),
+group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+left_join(sig_GC, by = c("tissue", "treatment", "met", "alias" = "group2")) %>%
+group_by(met) %>%
+filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>%
+left_join(met_dat) %>%
+left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+ungroup() %>%
+mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+heat.GC <- fc %>%
+group_by(Compound_Name, met) %>%
+mutate(log_norm = log2(mean_fc),
+log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>%
+ungroup() %>%
+#mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+pivot_wider(id_cols = c(Compound_Name, met),
+names_from = group,
+values_from = log_norm) %>%
+left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+as.data.frame()
+heat.GC <- fc %>%
+group_by(Compound_Name, met) %>%
+mutate(log_norm = log2(mean_fc),
+log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>%
+ungroup() %>%
+#mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+pivot_wider(id_cols = c(Compound_Name, met),
+names_from = group,
+values_from = log_norm) %>%
+#left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+as.data.frame()
+heat.GC <- fc %>%
+group_by(Compound_Name, met) %>%
+mutate(log_norm = log2(mean_fc),
+log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>%
+ungroup() %>%
+#mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+pivot_wider(id_cols = c(Compound_Name, Compound_Class, met),
+names_from = group,
+values_from = log_norm) %>%
+#left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+as.data.frame()
+rownames(heat.GC) <- heat.GC$met
+mat.heat.GC <- heat.GC %>%
+select(contains("leaf"), contains("fruit")) %>% as.matrix()
+annotation_row <- heat.GC %>%
+select(Compound_Class)
+rownames(annotation_row) <- heat.GC$met
+annotation_col <- sam_dat %>%
+filter(exp == 1, class == "sample") %>%
+distinct(tissue, treatment, genotype) %>%
+mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+as.data.frame()
+rownames(annotation_col) <- annotation_col$group
+annotation_col <- annotation_col %>%
+select(-group, treatment, tissue, genotype)
+heat.GC_signif <- means %>%
+distinct(tissue, treatment, alias, genotype, met) %>%
+#filter(genotype != "MoneyMaker") %>%
+left_join(sig_GC, by = c("met", "tissue", "treatment", "alias" = "group2")) %>%
+left_join(sig_mets) %>%
+filter(sig == T) %>%
+mutate(group = as_factor(str_c(tissue, treatment, genotype, sep = "_")),
+signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+pivot_wider(id_cols = c(met),
+names_from = group,
+values_from = signif) %>%
+left_join(met_dat) %>%
+left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+as.data.frame()
+heat_cols <- colnames(mat.heat.GC)
+rownames(heat.GC_signif) <- heat.GC_signif$met
+mat.heat.GC_signif <- heat.GC_signif %>%
+select(all_of(heat_cols)) %>% as.matrix()
+ann_colors = list(
+tissue = c(fruit = "red",leaf = "darkgreen"),
+genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+pheatmap.GC <- pheatmap(mat.heat.GC,
+#color = RColorBrewer::brewer.pal(7,"BuRd"),
+#cellwidth = 16,
+#cellheight = 4,
+#breaks = c(-2.5,-1.5, -0.5 ,0.5,1.5,2.5),
+cluster_rows = T,
+cluster_cols = T,
+annotation_names_row = F,
+show_rownames = F,
+annotation_row = annotation_row,
+annotation_col = annotation_col,
+display_numbers = mat.heat.GC_signif,
+number_color = "black",
+fontsize_number = 6,
+angle_col = 45,
+fontsize_col = 6,
+annotation_colors = ann_colors,
+filename = str_c(str_replace_all(Sys.Date(),"^.{2}|-",""),
+"cmQTL_val_1_heatmap_rel_tissue_wt.jpg",
+sep = "_")
)
-here(out,
-"cmQTL_val_1_2_sam_dat_GC.csv"
+pheatmap.GC <- pheatmap(mat.heat.GC,
+#color = RColorBrewer::brewer.pal(7,"BuRd"),
+#cellwidth = 16,
+#cellheight = 4,
+#breaks = c(-2.5,-1.5, -0.5 ,0.5,1.5,2.5),
+cluster_rows = T,
+cluster_cols = T,
+annotation_names_row = F,
+show_rownames = F,
+annotation_row = annotation_row,
+annotation_col = annotation_col,
+display_numbers = mat.heat.GC_signif,
+number_color = "black",
+fontsize_number = 6,
+angle_col = 45,
+fontsize_col = 6,
+annotation_colors = ann_colors,
+filename = here(out,
+"cmQTL_val_1_heatmap_rel_tissue_wt.jpg")
)
-View(isa_ms)
-View(isa_ms_tidy)
-1176-1145
-source("C:/Users/Micha/Nextcloud/Data_Plant/cmqtl_val1_arc/workflows/GC_MS_normalization/210927_primary_normalization_with_split.R", echo=TRUE)
-#PCA testing chunk####
-##https://www.intechopen.com/books/metabolomics-fundamentals-and-applications/processing-and-visualization-of-metabolomics-data-using-r
-#exp1
-pca_base <- features_all %>%
-ungroup() %>%
-filter(exp == 1) %>%
-arrange(batch_pca, daily_num) %>%
-mutate(log_imp = log2(imp)) %>%
-pivot_wider(id_cols = c(tissue, treatment, batch_GC, daily_num, class, machine_num_GC),
+GC_tidy <- means_tec_rep_lt %>%
+pivot_wider(id_cols = c( genotype, alias, treatment, tissue, LIMS_ID),
names_from = met,
-values_from = linnorm_fw_log) %>%
-filter(!is.na(m_1))
-pca <- summary(prcomp(
-pca_base%>%
-select(starts_with("m_"))))
-sam_vars <- pca_base %>% select(-starts_with("m_")) %>% colnames()
-pca_plot <- as_tibble(pca$x) %>%
-mutate(join_num = 1:nrow(pca$x)) %>%
-left_join(pca_base %>%
-mutate(join_num=1:nrow(pca$x))) %>%
-#left_join(sam_dat) %>%
-mutate(batch_GC= as_factor(batch_GC),
-class= as_factor(class)) %>%
-select(all_of(sam_vars), everything())
-exp_var <- as_tibble(pca[["importance"]])
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = tissue)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = tissue)) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1_tissue_lin_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = treatment)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = treatment)) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1_treatment_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
-dev.off()
-pca_base <- features_all %>%
+values_from = lev_t) %>%
+mutate(group = str_c(tissue, treatment, alias, sep = "_"))
+miss_per_clust <- GC_long %>%
+group_by(met, tissue) %>%
+summarise(na = sum(is.na(area)),
+n = n()) %>%
ungroup() %>%
-filter(exp == 1) %>%
-arrange(batch_pca, daily_num) %>%
-mutate(log_imp = log2(imp)) %>%
-pivot_wider(id_cols = c(tissue, treatment, batch_GC, daily_num, class, machine_num_GC),
-names_from = met,
-values_from = loess_norm_fw_log) %>%
-filter(!is.na(m_1))
-pca <- summary(prcomp(
-pca_base%>%
-select(starts_with("m_"))))
-sam_vars <- pca_base %>% select(-starts_with("m_")) %>% colnames()
-pca_plot <- as_tibble(pca$x) %>%
-mutate(join_num = 1:nrow(pca$x)) %>%
-full_join(pca_base %>%
-mutate(join_num=1:nrow(pca$x))) %>%
-#left_join(sam_dat) %>%
-mutate(batch_pca= as_factor(batch_GC),
-class= as_factor(class)) %>%
-select(all_of(sam_vars), everything())
-exp_var <- as_tibble(pca[["importance"]])
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = tissue)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = tissue)) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1_tissue_loess_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = treatment)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = treatment)) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1_treatment_loess_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
-pca_base <- features_all %>%
+mutate(percent_na_clust = na/n*100) %>%
+select(met, tissue, percent_na_clust)
+miss_any_treat <- miss_per_treat %>%
+group_by(met, tissue) %>%
+summarise(miss_comp_treat = if_else(any(percent_na == 100), T,F))
+means_tec_rep_lt <- means_tec_rep %>%
+group_by(met, tissue, treatment, genotype, alias) %>%
+mutate(lev_t = abs(log10(mean_tec_rep) - log10(median(mean_tec_rep)))) %>%
+ungroup()
+means_lt <- means_tec_rep_lt %>%
ungroup() %>%
-filter(exp == 2) %>%
-arrange(batch_pca, daily_num) %>%
-mutate(log_imp = log2(isnorm)) %>%
-pivot_wider(id_cols = c(tissue, treatment, batch_GC, daily_num, class, machine_num_GC),
-names_from = met,
-values_from = log_imp)
-pca <- summary(prcomp(
-pca_base%>%
-select(starts_with("m_"))))
-sam_vars <- pca_base %>% select(-starts_with("m_")) %>% colnames()
-pca_plot <- as_tibble(pca$x) %>%
-mutate(join_num = 1:nrow(pca$x)) %>%
-full_join(pca_base %>%
-mutate(join_num=1:nrow(pca$x))) %>%
-#left_join(sam_dat) %>%
-mutate(batch_pca= as_factor(batch_GC),
-class= as_factor(class)) %>%
-select(all_of(sam_vars), everything())
-exp_var <- as_tibble(pca[["importance"]])
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = class)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = class)) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = batch_GC)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = batch_GC)) +
-geom_text(aes(x = PC1, y = PC2,
-label = if_else(machine_num_GC %in% exclude_samples, machine_num_GC, ""))) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1_treatment_loess_outlier_exp2_loess_norm_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
-for(tiss in seq_along(tissues)){
-pca_base <- features_all %>%
+group_by(met, tissue, treatment, genotype, alias) %>%
+summarise(mean = mean(lev_t),
+sd = sd(lev_t),
+n = n()) %>%
ungroup() %>%
-filter(exp == 1, tissue == tissues[[tiss]]) %>%
-arrange(batch_pca, daily_num) %>%
-mutate(log_imp = log2(imp)) %>%
-pivot_wider(id_cols = c(tissue, treatment,batch_GC, daily_num, class, machine_num_GC),
+mutate(se = sd/sqrt(n))
+GC_tidy <- means_tec_rep_lt %>%
+pivot_wider(id_cols = c(genotype, alias, treatment, tissue, LIMS_ID),
names_from = met,
-values_from = loess_norm_fw_log) %>%
-filter(!is.na(m_1))
-pca <- summary(prcomp(
-pca_base%>%
-select(starts_with("m_"))))
-sam_vars <- pca_base %>% select(-starts_with("m_")) %>% colnames()
-pca_plot <- as_tibble(pca$x) %>%
-mutate(join_num = 1:nrow(pca$x)) %>%
-full_join(pca_base %>%
-mutate(join_num=1:nrow(pca$x))) %>%
-left_join(sam_dat) %>%
-mutate(batch_pca= as_factor(batch_GC),
-class= as_factor(class)) %>%
-select(all_of(sam_vars), everything())
-exp_var <- as_tibble(pca[["importance"]])
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = treatment)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = treatment)) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1", tissues[[tiss]], "treatment_loess_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
+values_from = lev_t) %>%
+mutate(group = str_c(tissue, treatment, alias, sep = "_"))
+GC_tidy_numeric <- GC_tidy %>%
+select(all_of(mets))
+GC_t <- map(.x = GC_tidy_numeric,
+.f = ~pairwise.t.test(x = .x ,
+g = GC_tidy$group,
+p.adjust.method = "none",
+pool.sd = F)) %>%
+map(.f = tidy)
+GC_t_tidy <- GC_t %>%
+map2(.y = names(GC_t), .f = ~.x %>% mutate(var = .y)) %>%
+purrr::reduce(bind_rows)
+sig_GC_lt<- GC_t_tidy %>%
+separate(group1, into = c("tissue1", "treatment1", "alias1"), sep = "_") %>%
+separate(group2, into = c("tissue2", "treatment2", "alias2"), sep = "_") %>%
+filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT",
+tissue1 == tissue2, treatment1 == treatment2) %>%
+group_by(var) %>%
+mutate(adj.p.value = p.adjust(p.value)) %>%
+#mutate(adj.p.value = p.value * 121) %>%
+ungroup() %>%
+mutate(p.signif = if_else(p.value <= 0.05, "*", "ns")) %>%
+left_join(means_lt, by = c("var" = "met", "treatment1" = "treatment", "alias1" = "alias", "tissue1" = "tissue")) %>%
+select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1 = mean, se1 = se, p.value, adj.p.value) %>%
+left_join(means_lt, by = c("var" = "met", "treatment2" = "treatment", "alias2" = "alias", "tissue2" = "tissue")) %>%
+select(var, treatment1, treatment2, alias1, alias2, tissue1, tissue2, p.signif, mean1, mean2 = mean, se1, se2 = se, p.value, adj.p.value) %>%
+group_by(treatment1, tissue1, var) %>%
+mutate(tot_val1 = mean1 + se1,
+tot_val2 = mean2 + se2,
+y.position = 1.1*(max(tot_val1, tot_val2))) %>%
+ungroup() %>%
+rename(met = var,
+group1 = alias1,
+group2 = alias2,
+treatment = treatment1,
+tissue = tissue1) %>%
+left_join(miss_per_treat) %>%
+mutate(met = as_factor(met),
+group1 = as_factor(group1),
+group2 = as_factor(group2),
+treatment = as_factor(treatment),
+tissue = as_factor(tissue),
+p.signif = if_else(percent_na >= 60, "ns", p.signif))
+sig_mets <- sig_GC_lt %>%
+filter(p.signif == "*") %>%
+distinct(met) %>%
+mutate(sig = T)
+norm_MM_1_lt <- means_lt %>%
+filter(genotype == "MoneyMaker", treatment == 1) %>%
+select(tissue, met, MM_mean = mean)
+fc_1_lt <- means_tec_rep_lt %>%
+#left_join(norm_MM_1) %>%
+#mutate(fc = mean_tec_rep/MM_mean) %>%
+group_by(tissue, treatment, alias, genotype, met) %>%
+summarise(mean_fc = mean(lev_t),
+sd = sd(lev_t),
+n = n()) %>%
+mutate(se = sd/sqrt(n),
+group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+left_join(sig_GC_lt, by = c("tissue", "treatment", "met", "alias" = "group2")) %>%
+group_by(met) %>%
+#filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>%
+left_join(met_dat) %>%
+arrange(Compound_Class, Compound_Name) %>%
+ungroup() %>%
+mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+fc_1_lt <- means_tec_rep_lt %>%
+#left_join(norm_MM_1) %>%
+#mutate(fc = mean_tec_rep/MM_mean) %>%
+group_by(tissue, treatment, alias, genotype, met) %>%
+summarise(mean_fc = mean(lev_t),
+sd = sd(lev_t),
+n = n()) %>%
+mutate(se = sd/sqrt(n),
+group = as_factor(str_c(tissue, treatment, genotype, sep = "_"))) %>%
+left_join(sig_GC_lt, by = c("tissue", "treatment", "met", "alias" = "group2")) %>%
+group_by(met) %>%
+#filter(any(p.signif == "*"), all(is.finite(mean_fc)), tissue != "flowers") %>%
+left_join(met_dat) %>%
+left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+ungroup() %>%
+mutate(genotype = as_factor(if_else(genotype == "MoneyMaker", glue("{genotype}"),glue("*{genotype}*"))),
+genotype = fct_relevel(genotype, c("MoneyMaker", "*panK4-1*", "*log2-1*", "*transp1-1*")))
+# CV analysis -------------------------------------------------------------
+library(bootstrap)
+theta <- function(x){
+sd(x)/mean(x)
}
-sam_vars <- pca_base %>% select(-starts_with("m_")) %>% colnames()
-pca_plot <- as_tibble(pca$x) %>%
-mutate(join_num = 1:nrow(pca$x)) %>%
-full_join(pca_base %>%
-mutate(join_num=1:nrow(pca$x))) %>%
-left_join(sam_dat) %>%
-mutate(batch_pca= as_factor(batch_pca),
-class= as_factor(class)) %>%
-select(all_of(sam_vars), everything())
-exp_var <- as_tibble(pca[["importance"]])
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = tissue)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = tissue)) +
-geom_text(aes(x = PC1, y = PC2,
-label = if_else(machine_num_GC == "21106rA_31", machine_num_GC, "")),
-nudge_y = -2,
-nudge_x = 4) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1_tissue_imp_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = batch_GC)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = batch_GC)) +
-geom_text(aes(x = PC1, y = PC2,
-label = if_else(machine_num_GC == "21106rA_31", machine_num_GC, "")),
-nudge_y = -2,
-nudge_x = 4) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1_batch_imp_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
-pca_plot %>%
-# filter(treatment == "HL") %>%
-ggplot()+
-geom_jitter(aes(x=PC1, y=PC2, color = class)) +
-stat_ellipse(aes(x=PC1, y=PC2, color = class)) +
-geom_text(aes(x = PC1, y = PC2,
-label = if_else(machine_num_GC == "21106rA_31", machine_num_GC, "")),
-nudge_y = -2,
-nudge_x = 4) +
-ylab(str_c("PC2 ", "(", exp_var$PC2[[2]]*100, "%)")) +
-xlab(str_c("PC1 ", "(", exp_var$PC1[[2]]*100, "%)"))
-ggsave(here(out,
-"cmQTL_val1_class_imp_PCA.jpg"
-),
-width = 183,
-height = 100,
-units = "mm",
-dpi = 300)
-#cmQTLval2
-pca_base <- features_all %>%
+cv <- means %>%
ungroup() %>%
-filter(exp == 2, batch_GC == "3_split") %>%
-arrange(batch_GC, daily_num) %>%
-mutate(log_imp = log2(isnorm)) %>%
-pivot_wider(id_cols = c(tissue, treatment, batch_GC, daily_num, class, machine_num_GC),
+group_by(genotype, tissue, alias, met) %>%
+summarise(grand_mean = mean(mean),
+grand_sd = sd(mean),
+n = n()) %>%
+mutate(cv = grand_sd/grand_mean)
+cv_jack <- means %>%
+group_by(genotype, tissue, alias, met) %>%
+summarise(cv = jackknife(mean, theta)$jack.values) %>%
+mutate(jack_rep = row_number()) %>%
+ungroup()
+cv_jack <- means %>%
+group_by(genotype, tissue, alias, met) %>%
+reframe(cv = jackknife(mean, theta)$jack.values) %>%
+mutate(jack_rep = row_number()) %>%
+ungroup()
+cv_jack_mean <- cv_jack %>%
+group_by(genotype,tissue, alias, met) %>%
+summarise(mean_cv = mean(cv),
+sd_cv = sd(cv),
+n = n()) %>%
+ungroup() %>%
+mutate(se_cv = sd_cv/sqrt(n))
+cv_jack_wide <- cv_jack %>%
+pivot_wider(id_cols = c(genotype, tissue, alias, jack_rep),
names_from = met,
-values_from = log_imp)%>%
-filter(!is.na(m_1))
-pca <- summary(prcomp(
-pca_base%>%
-select(starts_with("m_"))))
-sam_vars <- pca_base %>% select(-starts_with("m_")) %>% colnames()
-skip <- T
-if(skip == T) {
-} else{
-# Plot all metabolites ----------------------------------------------------
-plotmets <- features_all %>% distinct(met) %>% as_vector
-plotmet_labs <- plotmets %>% as_tibble() %>%
-left_join(met_dat, by= c("value" = "met")) %>%
-#mutate(peak_num = base::rank(HMDB_clear_name, ties.method = "first"),
-# dup = sum(peak_num),
-# HMDB_clear_name_unique = if_else(dup>1, str_c(HMDB_clear_name, peak_num), HMDB_clear_name)) %>%
-#ungroup() %>%
-select(Compound_Name) %>% as_vector()
-plot_out <- vector("list", length = length(plotmets))
-for (meta in seq_along(plotmets)) {
-plot_out [[meta]] <- features_all %>%
-mutate(xint = if_else(daily_num == 4, run_num_GC-3.5, max(run_num_GC)),
-is_miss = as_factor(if_else(is.na(area), T, F))) %>%
-filter(met == plotmets[[meta]]) %>%
-filter(class!="blank") %>%
-ggplot(aes(x=run_num_GC, y=isnorm)) +
-geom_point(aes(color = class, shape = is_miss)) +
-geom_point(aes(y=predict), color="black", size=0.1) +
-geom_vline(aes(xintercept = xint))+
-# facet_grid(rows = vars(treatment), cols = vars(rep), scales = "free") +
-ggtitle(label = plotmet_labs[[meta]])
-}
-pdf(file = here(out,"_loess_fit.pdf"))
-for (meta in seq_along(plotmets)) {
-print(plot_out[[meta]])
-}
-dev.off()
-plot_out <- vector("list", length = length(plotmets))
-for (meta in seq_along(plotmets)) {
-plot_out [[meta]] <- features_all %>%
-mutate(xint = if_else(daily_num == 4, run_num_GC-3.5, max(run_num_GC))) %>%
-filter(met == plotmets[[meta]]) %>%
-filter(class!="blank") %>%
-ggplot(aes(x=run_num_GC, y=loess_norm_fw_log)) +
-geom_point(aes(color=class)) +
-#geom_point(aes(y=predp), color="black", size=0.1) +
-geom_vline(aes(xintercept = xint))+
-# facet_grid(rows = vars(treatment), cols = vars(rep), scales = "free") +
-ggtitle(label = plotmet_labs[[meta]])
-}
-pdf(file = here(out,"_loess_norm_fw_log.pdf"))
-for (meta in seq_along(plotmets)) {
-print(plot_out[[meta]])
-}
-dev.off()
-for (meta in seq_along(plotmets)) {
-plot_out [[meta]] <- features_all %>%
-mutate(is_miss = as_factor(if_else(is.na(area), T, F))) %>%
-filter(met == plotmets[[meta]]) %>%
-filter(class!="blank") %>%
-ggplot(aes(x=daily_num, y=isnorm)) +
-geom_point(aes(color = class, shape = is_miss)) +
-geom_point(aes(y=predict), color="black", size=0.1) +
-facet_wrap(facets = vars(batch_GC), nrow = 1, scales = "free") +
-ggtitle(label = plotmet_labs[[meta]]) +
-theme(axis.ticks.y = element_blank(),
-axis.text.y = element_blank(),
-panel.spacing = unit(-0.1,"lines"))
-}
-pdf(file = here(out,"_loess_norm_no_scale.pdf"),
-height = 3.5)
-for (meta in seq_along(plotmets)) {
-print(plot_out[[meta]])
-}
+values_from = cv) %>%
+mutate(group = str_c(tissue, alias, sep = "_"))
+cv_jack_numeric <- cv_jack_wide %>%
+select(-c(genotype, tissue, alias, jack_rep, group))
+GC_cv_jack <- map(.x = cv_jack_numeric, .f = ~aov(.x ~ alias, data = cv_jack_wide))
+GC_jack_t <- map(.x = cv_jack_numeric, .f = ~pairwise.t.test(x = .x , g = cv_jack_wide$group, p.adjust.method = "none")) %>%
+map(.f = tidy)
+GC_jack_t_tidy <- GC_jack_t %>%
+map2(.y = names(GC_jack_t), .f = ~.x %>% mutate(var = .y)) %>%
+purrr::reduce(bind_rows)
+sig_GC_cv_groups <- GC_jack_t_tidy %>%
+separate(group1, into = c("tissue1", "alias1"), sep = "_") %>%
+separate(group2, into = c("tissue2", "alias2"), sep = "_") %>%
+filter(alias1 == "967514 MM WT" | alias2 == "967514 MM WT",
+tissue1 == tissue2) %>%
+group_by(var) %>%
+mutate(adj.p.value = p.adjust(p.value)) %>%
+ungroup() %>%
+mutate(p.signif = if_else(adj.p.value <= 0.05, "*", "ns")) %>%
+left_join(cv_jack_mean, by = c("var" = "met", "alias1" = "alias", "tissue1" = "tissue")) %>%
+select(var, alias1, alias2, p.signif, mean_cv1 = mean_cv, se_cv1 = se_cv, tissue1, tissue2, p.value, adj.p.value) %>%
+left_join(cv_jack_mean, by = c("var" = "met", "alias2" = "alias", "tissue2" = "tissue")) %>%
+select(var, alias1, alias2, p.signif, mean_cv1, mean_cv2 = mean_cv, se_cv1, se_cv2 = se_cv, tissue1, tissue2, p.value, adj.p.value) %>%
+group_by(var, tissue1) %>%
+mutate(tot_val1 = mean_cv1 + se_cv1,
+tot_val2 = mean_cv2 + se_cv2,
+y.position = 1.1*(max(tot_val1, tot_val2))) %>%
+ungroup() %>%
+left_join(genotypes, by = c("alias1" = "alias")) %>%
+rename(genotype1 = genotype) %>%
+left_join(genotypes, by = c("alias2" = "alias")) %>%
+rename(genotype2 = genotype) %>%
+rename(met = var,
+group1 = genotype1,
+group2 = genotype2,
+tissue = tissue1) %>%
+left_join(miss_any_treat) %>%
+left_join(miss_per_clust) %>%
+group_by(met, tissue) %>%
+mutate(p.signif = if_else(miss_comp_treat == T| percent_na_clust >= 50, "ns", p.signif),
+met = as_factor(met),
+group1 = as_factor(group1),
+group2 = as_factor(group2)) %>%
+ungroup()
+sig_mets_cv <- sig_GC_cv_groups %>%
+filter(p.signif == "*") %>%
+distinct(met, .keep_all = T) %>%
+mutate(sig = T) %>%
+select(met, sig)
+heat.GC <- cv_jack_mean%>%
+group_by(tissue, met) %>%
+mutate(total_norm = mean_cv/mean(mean_cv),
+log_norm = log2(total_norm),
+log_norm = if_else(is.infinite(log_norm), 0, log_norm),
+log_norm_level = (log_norm - mean(log_norm))/(max(log_norm)-min(log_norm))) %>%
+ungroup() %>%
+left_join(sig_mets_cv) %>%
+filter(sig == T) %>%
+mutate(group = as_factor(str_c(tissue, genotype, sep = "_"))) %>%
+pivot_wider(id_cols = c(met),
+names_from = group,
+values_from = log_norm_level) %>%
+left_join(met_dat) %>%
+left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+as.data.frame()
+rownames(heat.GC) <- heat.GC$met
+mat.heat.GC <- heat.GC %>%
+select(contains("leaf"), contains("fruit")) %>% as.matrix()
+annotation_row <- heat.GC %>%
+select(Compound_Class)
+rownames(annotation_row) <- heat.GC$met
+annotation_col <- sam_dat %>%
+filter(exp == 1, class == "Sample") %>%
+distinct(tissue, genotype) %>%
+mutate(group = as_factor(str_c(tissue, genotype, sep = "_"))) %>%
+as.data.frame()
+rownames(annotation_col) <- annotation_col$group
+annotation_col <- annotation_col %>%
+select(-group, tissue, genotype)
+heat.GC_signif <- cv_jack_mean %>%
+distinct(tissue, alias, genotype, met) %>%
+#filter(genotype != "MoneyMaker") %>%
+left_join(sig_GC_cv_groups, by = c("met", "tissue" = "tissue2", "alias" = "alias2")) %>%
+left_join(sig_mets_cv) %>%
+filter(sig == T) %>%
+mutate(group = as_factor(str_c(tissue, genotype, sep = "_")),
+signif = if_else(p.signif == "ns"| is.na(p.signif), "","X")) %>%
+pivot_wider(id_cols = c(met),
+names_from = group,
+values_from = signif) %>%
+left_join(met_dat) %>%
+left_join(GC_classes) %>%
+arrange(Compound_Class, Compound_Name) %>%
+as.data.frame()
+heat_cols <- colnames(mat.heat.GC)
+rownames(heat.GC_signif) <- heat.GC_signif$met
+mat.heat.GC_signif <- heat.GC_signif %>%
+select(all_of(heat_cols)) %>% as.matrix()
+ann_colors = list(
+tissue = c(fruit = "red",leaf = "darkgreen"),
+genotype = c(MoneyMaker = "yellow", `panK4-1` = "brown", `log2-1` = "blue", `transp1-1` = "grey"),
+treatment = c("0.4" = "red", "0.6" = "orange", "0.8" = "yellow", "1" = "green"))
+pheatmap.GC <- pheatmap(mat.heat.GC,
+#color = plasma(14),
+#cellwidth = 16,
+#cellheight = 4,
+#breaks = c(-6.5,-5.5,-4.5,-3.5,-2.5,-1.5, -0.5 ,0.5,1.5,2.5,3.5,4.5,5.5,6.5),
+cluster_rows = T,
+cluster_cols = T,
+annotation_names_row = F,
+show_rownames = F,
+annotation_row = annotation_row,
+annotation_col = annotation_col,
+display_numbers = mat.heat.GC_signif,
+number_color = "black",
+fontsize_number = 6,
+angle_col = 45,
+fontsize_col = 6,
+annotation_colors = ann_colors,
+filename = here(out,
+"cmQTL_val1_cv_heatmap_rel_tissue.jpg",
+)
+)
+pheatmap.GC <- pheatmap(mat.heat.GC,
+#color = plasma(14),
+#cellwidth = 16,
+#cellheight = 4,
+#breaks = c(-6.5,-5.5,-4.5,-3.5,-2.5,-1.5, -0.5 ,0.5,1.5,2.5,3.5,4.5,5.5,6.5),
+cluster_rows = T,
+cluster_cols = T,
+annotation_names_row = F,
+show_rownames = F,
+annotation_row = annotation_row,
+annotation_col = annotation_col,
+display_numbers = mat.heat.GC_signif,
+number_color = "black",
+fontsize_number = 6,
+angle_col = 45,
+fontsize_col = 6,
+annotation_colors = ann_colors,
+filename = here(out,
+"cmQTL_val1_cv_heatmap_rel_tissue.jpg")
+)
+cv %>%
+filter(genotype == "MoneyMaker") %>%
+mutate(`cv > 1` = if_else(cv > 1, "yes", "no")) %>%
+ggplot(aes(x = tissue, y = cv)) +
+geom_dotplot(aes(fill = `cv > 1`), stackdir = "center", binaxis = "y", binwidth = 0.1) +
+geom_hline(aes(yintercept = 1), color = "red") +
+theme(panel.background = element_rect(fill = "white"),
+panel.border = element_rect(color = "black",fill = NA),
+legend.position = "bottom",
+legend.text = element_text(size = 8),
+text = element_text(size = 10)) +
+ylim(c(-0.15, 2))
dev.off()
-for (meta in seq_along(plotmets)) {
-plot_out [[meta]] <- features_all %>%
-mutate(xint = if_else(daily_num == 4, run_num_GC-3.5, max(run_num_GC)),
-is_miss = as_factor(if_else(is.na(area), T, F))) %>%
-filter(met == plotmets[[meta]]) %>%
-filter(class!="blank") %>%
-ggplot(aes(x=run_num_GC, y=isnorm)) +
-geom_point(aes(color = class, shape = is_miss)) +
-geom_point(aes(y=predp), color="black", size=0.1) +
-geom_vline(aes(xintercept = xint))+
-# facet_grid(rows = vars(treatment), cols = vars(rep), scales = "free") +
-ggtitle(label = plotmet_labs[[meta]])
-}
-pdf(file = here(out,"_linnorm.pdf"))
-for (meta in seq_along(plotmets)) {
-print(plot_out[[meta]])
-}
+cv %>%
+filter(genotype == "MoneyMaker") %>%
+mutate(`cv > 1` = if_else(cv > 1, "yes", "no")) %>%
+ggplot(aes(x = tissue, y = cv)) +
+geom_dotplot(aes(fill = `cv > 1`), stackdir = "center", binaxis = "y", binwidth = 0.1) +
+geom_hline(aes(yintercept = 1), color = "red") +
+theme(panel.background = element_rect(fill = "white"),
+panel.border = element_rect(color = "black",fill = NA),
+legend.position = "bottom",
+legend.text = element_text(size = 8),
+text = element_text(size = 10)) +
+ylim(c(-0.15, 2))
dev.off()
-rescale <- features_all %>%
-group_by(met, tissue, exp) %>%
-summarise(rescale = median(imp))
-sam_vars <- c(#"plantline","extraction_num", "sample_num", #in case I add them back
-"alias", "LIMS_ID",
-"treatment", "tissue", "batch_GC", "run_date_GC",
-"machine_num_GC",
-"class", "run_num_GC", "sample_weight", "exp", "genotype")
-features_out <- features_all %>%
-left_join(rescale) %>%
-filter(class != "blank",
-(component %in% take_split & str_detect(batch_GC, "split"))|
-(!component %in% take_split & str_detect(batch_GC, "split", negate = T))|
-tissue == "leaves") %>%
-group_by(met, tissue, exp) %>%
-mutate(loess_norm_med = loess_norm_fw/median(loess_norm_fw),
-rescaled = loess_norm_med*rescale) %>%
-select(all_of(sam_vars), met, Compound_Name, loess_norm_fw, loess_norm_med, area, rescaled)
-#features_out %>%
-# filter(met == "m_02177", tissue == "fruits", exp == 1) %>%
-# ggplot(aes(x = run_num_GC, y = loess_norm_fw)) +
-# geom_point()
-}
-write_csv(features_out,
-here(out,
-"cmQTL_val_1_2_feat_dat_GC.csv"
-))
-here()
-write_csv(features_out,
-here(out,
-"cmQTL_val_1_2_feat_dat_GC.csv"
-))
-features_out
-write_csv(features_out,
-here(out,
-"cmQTL_val_1_2_feat_dat_GC.csv"
-))
-getwd()
-here()
-write_csv(features_out,
-here(out,
-"cmQTL_val_1_2_feat_dat_GC.csv"
-))
-here::dr_here()
-here::i_am("workflows/GC_MS_normalization/210927_primary_normalization_with_split.R")
-write_csv(features_out,
-here(out,
-"cmQTL_val_1_2_feat_dat_GC.csv"
-))
-here(out,
-"cmQTL_val_1_2_feat_dat_GC.csv"
-)
-out
-out <- here("runs/GC-MS normalization")
-out
-write_csv(features_out,
-here(out,
-"cmQTL_val_1_2_feat_dat_GC.csv"
-))
-write_csv(met_dat,
-here(out,
-"cmQTL_val_1_2_met_dat_GC.csv"
-))
-write_csv(sam_dat,
-here(out,
-"cmQTL_val_1_2_sam_dat_GC.csv"
-))
-write_csv(missingness,
-here(out,
-"cmQTL_val_1_2_miss_GC.csv"
-))
-View(features_out)
-source("C:/Users/Micha/Nextcloud/Data_Plant/cmqtl_val1_arc/workflows/GC_MS_normalization/210927_primary_normalization_with_split.R", echo=TRUE)
-View(features_out)
-View(features_all)
-features_out <- features_all %>%
-left_join(rescale) %>%
-filter(class != "blank",
-(component %in% take_split & str_detect(batch_GC, "split"))|
-(!component %in% take_split & str_detect(batch_GC, "split", negate = T))|
-tissue == "leaves") %>%
-group_by(met, tissue, exp) %>%
-mutate(loess_norm_med = loess_norm_fw/median(loess_norm_fw),
-rescaled = loess_norm_med*rescale) %>%
-select(all_of(sam_vars), met, Compound_Name, loess_norm_fw, loess_norm_med, area, rescaled)
-}
-write_csv(features_out,
-here(out,
-"cmQTL_val_1_2_feat_dat_GC.csv"
-))
-write_csv(met_dat,
-here(out,
-"cmQTL_val_1_2_met_dat_GC.csv"
-))
-write_csv(sam_dat,
-here(out,
-"cmQTL_val_1_2_sam_dat_GC.csv"
-))
-write_csv(missingness,
-here(out,
-"cmQTL_val_1_2_miss_GC.csv"
-))
-source("C:/Users/Micha/Nextcloud/Data_Plant/cmqtl_val1_arc/workflows/GC_MS_normalization/210927_primary_normalization_with_split.R", echo=TRUE)
-View(features_out)
-View(met_dat)
+cv %>%
+filter(genotype == "MoneyMaker") %>%
+mutate(`cv > 1` = if_else(cv > 1, "yes", "no")) %>%
+ggplot(aes(x = tissue, y = cv)) +
+geom_dotplot(aes(fill = `cv >