From 540b6c5bf7e4ca7689ae1eeb5333b2428e4839bd Mon Sep 17 00:00:00 2001
From: Tanushree Tunstall <tanu@tunstall.in>
Date: Thu, 29 Oct 2020 18:50:25 +0000
Subject: [PATCH] generating plots for all sample types withoout stats

---
 Header_TT.R |   2 +
 boxplot.R   | 351 ++++++++++++++++++++++++++++++++--------------------
 2 files changed, 217 insertions(+), 136 deletions(-)

diff --git a/Header_TT.R b/Header_TT.R
index 8e3f031..787b6c6 100755
--- a/Header_TT.R
+++ b/Header_TT.R
@@ -7,3 +7,5 @@ library(rstatix)
 library(Hmisc)
 library(qwraps2)
 source("legend_adjustment.R")
+
+# https://www.datanovia.com/en/blog/how-to-add-p-values-onto-a-grouped-ggplot-using-the-ggpubr-r-package/#perform-all-pairwise-comparisons
\ No newline at end of file
diff --git a/boxplot.R b/boxplot.R
index 5d576ea..035ea4d 100644
--- a/boxplot.R
+++ b/boxplot.R
@@ -22,150 +22,229 @@ outfile_bp = paste0("boxplots_", my_sample_type, ".pdf")
 output_boxplot = paste0(outdir_plots, outfile_bp); output_boxplot
 
 #===============================
-# data assignment for stats
+# data assignment for plots
 #================================
-wf = serum_wf[serum_wf$flustat == 1,]
-lf = serum_lf[serum_lf$flustat == 1,]
-lf$timepoint = paste0("t", lf$timepoint)
+#wf_fp_npa = npa_wf[npa_wf$flustat == 1,]
+lf_fp_npa = npa_lf[npa_lf$flustat == 1,]
+lf_fp_npa$timepoint = paste0("t", lf_fp_npa$timepoint)
+lf_fp_npa$timepoint = as.factor(lf_fp_npa$timepoint)
+lf_fp_npa$obesity = as.factor(lf_fp_npa$obesity)
 
-#==================
-# Data for plotting
-#==================
-# convert these to factor
-lf$obesity = as.factor(lf$obesity)
-lf$timepoint = as.factor(lf$timepoint)
+#wf_fp_sam = samm_wf[samm_wf$flustat == 1,]
+lf_fp_sam = sam_lf[sam_lf$flustat == 1,]
+lf_fp_sam$timepoint = paste0("t", lf_fp_sam$timepoint)
+lf_fp_sam$timepoint = as.factor(lf_fp_sam$timepoint)
+lf_fp_sam$obesity = as.factor(lf_fp_sam$obesity)
+
+#wf_fp_serum = serum_wf[serum_wf$flustat == 1,]
+lf_fp_serum = serum_lf[serum_lf$flustat == 1,]
+lf_fp_serum$timepoint = paste0("t", lf_fp_serum$timepoint)
+lf_fp_serum$timepoint = as.factor(lf_fp_serum$timepoint)
+lf_fp_serum$obesity = as.factor(lf_fp_serum$obesity)
 
 ########################################################################
-my_comparisons <- list( c("0", "1") )
-########################################################################
-lf_test = lf[lf$mediator == "eotaxin",]
-
-str(lf_test)
-p = ggplot(lf_test, aes(x = timepoint, y = value, fill = obesity  ))  + 
-  geom_boxplot(width = 0.5)
-p
-# see default ggplot palette
-ggplot_build(p)$data
-
-p + stat_compare_means(comparisons = my_comparisons
-                       , method = "wilcox.test"
-                       , paired = F
-                       , label = "p.format") 
-
-########################################################################
-library(ggpubr)
-library(rstatix)
-
-# stats
-stat.test <- lf_test %>%
-  group_by(timepoint) %>%
-  wilcox_test(value ~ obesity) %>%
-  adjust_pvalue(method = "bonferroni") %>%
-  add_significance("p.adj")
-stat.test
-
-
-# add stats to grouped boxplots
-str(lf_test)
-lf_test$obesity = as.factor(lf_test$obesity) # ensure factor
-
-bxp <- ggboxplot(lf_test, x = "timepoint", y = "value", 
-  color = "obesity", palette = c("#00BFC4", "#F8766D"))
-bxp
-
-stat.test <- stat.test %>%
-  add_xy_position(x = "timepoint", dodge = 0.8)
-bxp + stat_pvalue_manual(
-  stat.test,  label = "p", tip.length = 0
-)
-################################
-lf_test2 = lf[order(lf$mediator),] # 2 meds
-lf_test2 = lf_test2[1:798,]
-table(lf_test2$mediator)
-str(lf_test2)
-
-stat.test <- lf_test2 %>%
-  group_by(timepoint, mediator) %>%
-  wilcox_test(value ~ obesity, paired = F) %>%
-  #adjust_pvalue(method = "bonferroni") %>%
-  #add_significance("p.adj")
-  add_significance("p")
-
-stat.test
-
-bxp <- ggboxplot(lf_test2, x = "timepoint", y = "value", 
-                 color = "obesity", palette = c("#00BFC4", "#F8766D")) +
-  facet_wrap(~mediator, scales = "free_y")
-
-bxp
-
-stat.test <- stat.test %>%
-  add_xy_position(x = "timepoint", dodge = 0.8)
-bxp + stat_pvalue_manual(stat.test,  label = "p.signif", tip.length = 0)
-
-########################################################################
-lf_comp = lf[-which(is.na(lf$value)),]
-lf_test_comp = lf_comp[lf_comp$mediator == "eotaxin",]
-str(lf_comp)
-
-p2 = ggplot(lf_test_comp, aes(x = timepoint, y = value, fill = obesity ))  + 
-  geom_boxplot(width = 0.5)
-p2 
-
-########################################################################
-# Output plots as one pdf
-cat("Output plots will be in:", output_boxplot)
-pdf(output_boxplot, width = 25, height = 15)
-
-#stats data
-str(lf)
-stat.test <- lf %>%
-  group_by(timepoint, mediator) %>%
-  wilcox_test(value ~ obesity, paired = F) %>%
-  #adjust_pvalue(method = "bonferroni") %>%
-  #add_significance("p.adj")
-  add_significance("p")
-stat.test
-
-bxp <- ggboxplot(lf, x = "timepoint", y = "value", 
-                 color = "obesity", palette = c("#00BFC4", "#F8766D")) +
-  facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+
-  scale_y_log10()
-bxp
-
-bxp + stat.test <- stat.test %>%
-  add_xy_position(x = "timepoint", dodge = 0.8)
-
-bxp + stat_pvalue_manual(stat.test,  label = "p.signif", tip.length = 0)
-
-dev.off()
-
-######output 
-
-# Output plots as one pdf
 cat("Output plots will be in:", output_boxplot)
 pdf(output_boxplot, width = 20, height = 15)
-my_title1 = "Boxplots: serum"
 
-bxp <- ggboxplot(lf, x = "timepoint", y = "value", 
+#=======================================================================
+#                         SAM
+#=======================================================================
+if (is.factor(lf_fp_sam$timepoint) && is.factor(lf_fp_sam$timepoint)){
+  cat ("PASS: required groups are factors")
+}
+table(lf_fp_sam$mediator)
+lf_fp_sam = lf_fp_sam[!lf_fp_sam$mediator == "vitd",]
+
+#------------------------------------------
+title_sam_linear = "Boxplot: sam (Linear)"
+#-----------------------------------------
+bxp_sam_linear <- ggboxplot(lf_fp_sam, x = "timepoint", y = "value", 
                  color = "obesity", palette = c("#00BFC4", "#F8766D")) +
   facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+
-  scale_y_log10()
-
-bxp + theme(axis.text.x = element_text(size = 15)
-            , axis.text.y = element_text(size = 15
-                                         , angle = 0
-                                         , hjust = 1
-                                         , vjust = 0)
-            , axis.title.x = element_text(size = 15)
-            , axis.title.y = element_text(size = 15)
-            , plot.title = element_text(size = 20, hjust = 0.5)
-            , strip.text.x = element_text(size = 15, colour = "black")
-            , legend.title = element_text(color = "black", size = 20)
-            , legend.text = element_text(size = 15)
-            , legend.direction = "horizontal") +
-  labs(title = my_title1
+  #scale_y_log10() + 
+  theme(axis.text.x = element_text(size = 15)
+                          , axis.text.y = element_text(size = 15
+                                                       , angle = 0
+                                                       , hjust = 1
+                                                       , vjust = 0)
+                          , axis.title.x = element_text(size = 15)
+                          , axis.title.y = element_text(size = 15)
+                          , plot.title = element_text(size = 20, hjust = 0.5)
+                          , strip.text.x = element_text(size = 15, colour = "black")
+                          , legend.title = element_text(color = "black", size = 20)
+                          , legend.text = element_text(size = 15)
+                          , legend.direction = "horizontal") +
+  labs(title = title_sam_linear
        , x = ""
-       , y = "Levels (Log10)") 
-#shift_legend2(bxp)
+       , y = "Levels")
+
+bxp_sam_linear
+
+#------------------------------------
+title_sam_log = "Boxplot: sam (Log)"
+#-----------------------------------
+
+bxp_sam_log <- ggboxplot(lf_fp_sam, x = "timepoint", y = "value", 
+                            color = "obesity", palette = c("#00BFC4", "#F8766D")) +
+  facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+
+  scale_y_log10() + 
+  theme(axis.text.x = element_text(size = 15)
+        , axis.text.y = element_text(size = 15
+                                     , angle = 0
+                                     , hjust = 1
+                                     , vjust = 0)
+        , axis.title.x = element_text(size = 15)
+        , axis.title.y = element_text(size = 15)
+        , plot.title = element_text(size = 20, hjust = 0.5)
+        , strip.text.x = element_text(size = 15, colour = "black")
+        , legend.title = element_text(color = "black", size = 20)
+        , legend.text = element_text(size = 15)
+        , legend.direction = "horizontal") +
+  labs(title = title_sam_log
+       , x = ""
+       , y = "Levels (Log)")
+
+bxp_sam_log
+
+#=======================================================================
+#                        SERUM
+#=======================================================================
+if (is.factor(lf_fp_serum$timepoint) && is.factor(lf_fp_serum$timepoint)){
+  cat ("PASS: required groups are factors")
+}
+
+table(lf_fp_serum$mediator)
+
+#------------------------------------------
+title_serum_linear = "Boxplot: serum (Linear)"
+#-----------------------------------------
+bxp_serum_linear <- ggboxplot(lf_fp_serum, x = "timepoint", y = "value", 
+                              color = "obesity", palette = c("#00BFC4", "#F8766D")) +
+  facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+
+  #scale_y_log10() + 
+  theme(axis.text.x = element_text(size = 15)
+        , axis.text.y = element_text(size = 15
+                                     , angle = 0
+                                     , hjust = 1
+                                     , vjust = 0)
+        , axis.title.x = element_text(size = 15)
+        , axis.title.y = element_text(size = 15)
+        , plot.title = element_text(size = 20, hjust = 0.5)
+        , strip.text.x = element_text(size = 15, colour = "black")
+        , legend.title = element_text(color = "black", size = 20)
+        , legend.text = element_text(size = 15)
+        , legend.direction = "horizontal") +
+  labs(title = title_serum_linear
+       , x = ""
+       , y = "Levels")
+
+bxp_serum_linear
+
+#------------------------------------
+title_serum_log = "Boxplot: serum (Log)"
+#-----------------------------------
+
+bxp_serum_log <- ggboxplot(lf_fp_serum, x = "timepoint", y = "value", 
+                           color = "obesity", palette = c("#00BFC4", "#F8766D")) +
+  facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+
+  scale_y_log10() + 
+  theme(axis.text.x = element_text(size = 15)
+        , axis.text.y = element_text(size = 15
+                                     , angle = 0
+                                     , hjust = 1
+                                     , vjust = 0)
+        , axis.title.x = element_text(size = 15)
+        , axis.title.y = element_text(size = 15)
+        , plot.title = element_text(size = 20, hjust = 0.5)
+        , strip.text.x = element_text(size = 15, colour = "black")
+        , legend.title = element_text(color = "black", size = 20)
+        , legend.text = element_text(size = 15)
+        , legend.direction = "horizontal") +
+  labs(title = title_serum_log
+       , x = ""
+       , y = "Levels (Log)")
+
+bxp_serum_log
+#=======================================================================
+#                         NPA
+#=======================================================================
+if (is.factor(lf_fp_npa$timepoint) && is.factor(lf_fp_npa$timepoint)){
+  cat ("PASS: required groups are factors")
+}
+
+table(lf_fp_npa$mediator)
+head(lf_fp_npa$value[lf_fp_npa$mediator == "vitd"])
+lf_fp_npa = lf_fp_npa[!lf_fp_npa$mediator == "vitd",]
+
+#------------------------------------------
+title_npa_linear = "Boxplot: NPA (Linear)"
+#-----------------------------------------
+bxp_npa_linear <- ggboxplot(lf_fp_npa, x = "timepoint", y = "value", 
+                            color = "obesity", palette = c("#00BFC4", "#F8766D")) +
+  facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = T)+
+  #scale_y_log10() + 
+  theme(axis.text.x = element_text(size = 15)
+        , axis.text.y = element_text(size = 15
+                                     , angle = 0
+                                     , hjust = 1
+                                     , vjust = 0)
+        , axis.title.x = element_text(size = 15)
+        , axis.title.y = element_text(size = 15)
+        , plot.title = element_text(size = 20, hjust = 0.5)
+        , strip.text.x = element_text(size = 15, colour = "black")
+        , legend.title = element_text(color = "black", size = 20)
+        , legend.text = element_text(size = 15)
+        , legend.direction = "horizontal") +
+  labs(title = title_npa_linear
+       , x = ""
+       , y = "Levels")
+
+bxp_npa_linear
+
+#------------------------------------
+title_npa_log = "Boxplot: NPA (Log)"
+#-----------------------------------
+
+bxp_npa_log <- ggboxplot(lf_fp_npa, x = "timepoint", y = "value", 
+                         color = "obesity", palette = c("#00BFC4", "#F8766D")) +
+  facet_wrap(~mediator, nrow = 7, ncol = 5, scales = "free_y", shrink = F)+
+  scale_y_log10() + 
+  theme(axis.text.x = element_text(size = 15)
+        , axis.text.y = element_text(size = 15
+                                     , angle = 0
+                                     , hjust = 1
+                                     , vjust = 0)
+        , axis.title.x = element_text(size = 15)
+        , axis.title.y = element_text(size = 15)
+        , plot.title = element_text(size = 20, hjust = 0.5)
+        , strip.text.x = element_text(size = 15, colour = "black")
+        , legend.title = element_text(color = "black", size = 20)
+        , legend.text = element_text(size = 15)
+        , legend.direction = "horizontal") +
+  labs(title = title_npa_log
+       , x = ""
+       , y = "Levels (Log)")
+
+bxp_npa_log
 dev.off()
+#==========================================================================
+#------------------------------------
+title_npa_log_stats = "Boxplot: NPA (Log) + stats"
+#-----------------------------------
+stat.test <- lf_fp_npa %>%
+  group_by(timepoint, mediator) %>%
+  wilcox_test(value ~ obesity, paired = F) %>%
+  add_significance("p")
+stat.test
+
+stat.test <- stat.test %>%
+  add_xy_position(x = "timepoint", dodge = 0.8)
+
+
+bxp_npa_linear + stat_pvalue_manual(stat.test,  label = "p.signif", tip.length = 0)
+
+bxp_npa_log + stat_pvalue_manual(stat.test,  label = "p.signif", tip.length = 0)
+
+
+
+dev.off()
\ No newline at end of file