diff --git a/boxplot_stat_function.R b/boxplot_stat_function.R
index dc209c6..8233c44 100644
--- a/boxplot_stat_function.R
+++ b/boxplot_stat_function.R
@@ -5,28 +5,27 @@ getwd()
 ############################################################
 # TASK: boxplots with stats
 ############################################################
-doMyPlots <- function(x, title  = "TEST") {
-    mediators = levels(as.factor(x$mediator))
+doMyPlots <- function(df) {
+    mediators = levels(as.factor(df$mediator))
     
     plots <- list()
     
     for (i in mediators) {
       cat("Plotting:", i, "\n")
-      single=x[x$mediator==i,]
+      single=df[df$mediator==i,]
       
       max_y = max(single$value, na.rm = T)
       cat("Plotting:", i, "max_y:", max_y, "\n")
-      
+    
       #----------
       # boxplot
       #----------
-      p2 = ggboxplot(single
-                     , x = "timepoint"
-                     , y = "value"
-                     , color = "obesity"
-                     #, palette = c("blue", "red")
-                     , palette = c("#00BFC4", "#F8766D")) +
-        
+      p = ggplot(single)+ geom_boxplot(aes(x = timepoint
+                                           , y = value
+                                           , color = obesity
+                                           #, palette = c("#00BFC4", "#F8766D")
+                                           ))+
+        scale_colour_manual(values=c("blue", "red")) + 
       theme(axis.text.x = element_text(size = 15)
             , axis.text.y = element_text(size = 15
                                          , angle = 0
@@ -38,31 +37,29 @@ doMyPlots <- function(x, title  = "TEST") {
             , plot.subtitle = element_text(size = 20, hjust = 0.5)
             , plot.title = element_text(size = 20, hjust = 0.5)) +
         labs(title = i)
+      
       #--------
       # stats
       #---------
-      stat_npa2 <- single %>%
+      stat_df <- single %>%
         group_by(timepoint, mediator) %>%
         wilcox_test(value ~ obesity, paired = F) %>%
         add_significance("p")
-      stat_npa2
+      stat_df 
       
-      stat_npa2 <- stat_npa2 %>%
+      stat_df <- stat_df %>%
         add_xy_position(x = "timepoint", dodge = 0.8)
-      
-      p2 = p2 + stat_pvalue_manual(stat_npa2
+      cat("printing p")
+      p2 = p + stat_pvalue_manual(stat_df
                                    #, y.position = max_y
                                    , label = "{p} {p.signif}"
-                                   , hide.ns=T
+                                   , hide.ns = T
                                    , tip.length = 0)+
-        scale_y_continuous(expand = expansion(mult = c(0.05, 0.25))) +
-        scale_x_discrete()
-      
+        scale_y_continuous(expand = expansion(mult = c(0.05, 0.25)))
       plots[[i]] <- p2
       
     }
     return(plots)
-
   }
 
 ###############################################################
@@ -75,7 +72,7 @@ source("boxplot_linear.R")
 #=============
 outfile_bp = paste0("boxplots_TEST", ".pdf")
 output_boxplot_stats = paste0(outdir_plots, outfile_bp); output_boxplot_stats
-pdf(output_boxplot_stats, width=22, height=16)
+#pdf(output_boxplot_stats, width=22, height=16)
 
 med_names = c("eotaxin3", "il12p70", "itac", "il13")
 lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,]
@@ -83,8 +80,8 @@ lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,]
 foo=doMyPlots(lf_test)
 #baz=cowplot::plot_grid(plotlist=foo, align = 'hv', ncol=2, nrow=2)
 baz
-dev.off()
-
+#dev.off()
+##################################################################
 #=========================
 # entire sample_type data
 #=========================
@@ -93,9 +90,24 @@ outfile_bp = paste0("boxplots_stats_", my_sample_type, ".pdf")
 output_boxplot_stats = paste0(outdir_plots, outfile_bp); output_boxplot_stats
 pdf(output_boxplot_stats, width=22, height=16)
 
-bar=doMyPlots(lf_fp_npa)
-#bang=cowplot::plot_grid(plotlist=bar, align = 'hv', ncol=7, nrow=5)
-#bang
-
-
-dev.off()
\ No newline at end of file
+#===========================
+# combine: ggpubr::ggarrange
+#===========================
+plots_npa = doMyPlots(lf_fp_npa)
+npa_plot = ggpubr::ggarrange(plotlist = plots_npa
+                             , align = "hv"
+                             , ncol = 7
+                             , nrow = 5
+                             , common.legend = T)
+#npa_plot
+npa_plot_annot = annotate_figure(npa_plot
+                , top = text_grob(my_sample_type, color = "purple", face = "bold", size = 14)
+                , bottom = text_grob("Mosaic data\nFlu positive adults (n=133)"
+                                   , color = "blue"
+                                   , hjust = 1, x = 0.98, face = "italic", size = 10)
+                , left = text_grob("Levels (pg/ml)", color = "black", rot = 90, size = 18)
+                #, right = "I'm done, thanks :-)!",
+                #, fig.lab = "Figure 1", fig.lab.face = "bold"
+)
+npa_plot_annot
+dev.off()