#!/usr/bin/Rscript getwd() setwd("~/git/mosaic_2020/") getwd() ############################################################ # TASK: boxplots at T1 # FIXME: currently not rendering, problem with NAs for stats? ############################################################ #============= # Input #============= source("boxplot_linear.R") ####################################################### med_names = c("eotaxin3", "il12p70", "itac", "il13") lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,] mediators = levels(as.factor(lf_test$mediator)) mediators = levels(as.factor(lf_fp_npa$mediator)) lf_test = lf_fp_npa plots <- list() for (i in mediators) { single=lf_test[lf_test$mediator==i,] max_y = max(single$value, na.rm = T) cat("Plotting:", i, "max_y:", max_y, "\n") p2 = 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.x = element_blank() , axis.text.y = element_text(size = 15 , angle = 0 , hjust = 1 , vjust = 0) , axis.title.x = element_blank() , axis.title.y = element_blank() , legend.position = "none" , plot.subtitle = element_text(size = 20, hjust = 0.5) , plot.title = element_text(size = 20, hjust = 0.5)) + labs(title = i #, subtitle = "test2" ) stat_npa2 <- single %>% group_by(timepoint, mediator) %>% wilcox_test(value ~ obesity, paired = F) %>% add_significance("p") stat_npa2 stat_npa2 <- stat_npa2 %>% add_xy_position(x = "timepoint", dodge = 0.8) p2 = p2 + stat_pvalue_manual(stat_npa2 #, y.position = max_y , label = "{p} {p.signif}" , hide.ns=T , tip.length = 0)+ scale_y_continuous(expand = expansion(mult = c(0.05, 0.25))) + plots[[i]] <- p2 } #============================ # combine: cowplot_plot_grid #============================ #cowplot::plot_grid(plotlist=plots, align = 'hv', ncol=2, nrow=2) pdf("test.pdf", width = 22, height = 16) cowplot::plot_grid(plotlist=plots, align = 'hv', ncol=7, nrow=5) dev.off() #=========================== # combine: ggpubr::ggarrange #=========================== pdf("test2.pdf", width = 22, height = 16) npa_plot<- ggpubr::ggarrange(plotlist = plots, align = "hv" , ncol = 7 , nrow = 5 , common.legend = T) npa_plot annotate_figure(npa_plot, top = text_grob("NPA", 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" ) ################################################################## #======= # facet #======= #-------- # wilcox stats #-------- stat_npa3 <- lf_test %>% group_by(timepoint, mediator) %>% wilcox_test(value ~ obesity, paired = F) %>% add_significance("p") stat_npa3 stat_npa3$p_format = round(stat_npa3$p, 3) stat_npa3 <- stat_npa3 %>% add_xy_position(x = "timepoint", dodge = 0.8) head(stat_npa3) #-------- # summary stats #-------- my_summary = lf_test %>% group_by(timepoint, mediator) %>% get_summary_stats(value) my_summary my_max_df = subset(my_summary , select = c("mediator", "timepoint" #, "obesity" , "max")) head(my_max_df); head(stat_npa3) #----------------------------- # merge my_max and stat_npa3 #----------------------------- head(my_max_df) merging_cols = intersect(names(stat_npa3), names(my_max_df)); merging_cols stat_npa3_v2 = merge(stat_npa3, my_max_df, by = merging_cols, all.x = T) stat_npa3_v2$my_y_pos = (stat_npa3_v2$max)*1.2 stat_npa3_v2$my_y_pos head(stat_npa3_v2); head(stat_npa3) p3 = ggplot(lf_test)+ geom_boxplot(aes(x = timepoint , y = value , color = obesity))+ scale_colour_manual(values=c("blue", "red")) + 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_blank() , axis.title.y = element_blank() , strip.text.x = element_text(size = 20,margin = margin(0.05,0,0.07,0, "cm")) , legend.position = "top" , legend.title = element_text(color = "black", size = 20) , legend.text = element_text(size = 15) , plot.subtitle = element_text(size = 20, hjust = 0.5) , plot.title = element_text(size = 20, hjust = 0.5)) + labs(title ="NPA")+ facet_wrap(~mediator, scales = "free", nrow = 2, ncol = 2) #p3 p4 = p3 + stat_pvalue_manual(stat_npa3_v2 , y.position = "my_y_pos" #, y.position = 50 , label = "{p_format} {p.signif}" , hide.ns=T , tip.length = 0)+ scale_y_continuous(expand = expansion(mult = c(0.05, 0.25))) p4 #======================================= p3_build = ggplot_build(p3) p3_build$layout$panel_scales_y get_facet_ymax <- function (x){ ret <- x$layout$panel_scales_y y_max = NULL for (i in 1:length(ret)){ print(i) y_max_i <- max(x$layout$panel_scales_y[[i]]$range$range) y_max = c(y_max, y_max_i) } return(y_max) } y_max_l = get_facet_ymax(p3_build); y_max_l #=============================================== # facet wrap on actual data pdf("boxplot_stats.pdf", width = 20, height = 15) #======= # facet #======= #-------- # wilcox stats #-------- stat_npa3 <- lf_fp_npa %>% group_by(timepoint, mediator) %>% wilcox_test(value ~ obesity, paired = F) %>% add_significance("p") stat_npa3 stat_npa3$p_format = round(stat_npa3$p, 3) stat_npa3 <- stat_npa3 %>% add_xy_position(x = "timepoint", dodge = 0.8) head(stat_npa3) #-------- # summary stats #-------- my_summary = lf_fp_npa %>% group_by(timepoint, mediator) %>% get_summary_stats(value) my_summary my_max_df = subset(my_summary , select = c("mediator", "timepoint", "max")) head(my_max_df); head(stat_npa3) #----------------------------- # merge my_max and stat_npa3 #----------------------------- head(my_max_df) merging_cols = intersect(names(stat_npa3), names(my_max_df)); merging_cols stat_npa3_v2 = merge(stat_npa3, my_max_df, by = merging_cols, all.x = T) stat_npa3_v2$my_y_pos = (stat_npa3_v2$max)*1.2 stat_npa3_v2$my_y_pos head(stat_npa3_v2); head(stat_npa3) p3 = ggplot(lf_fp_npa)+ geom_boxplot(aes(x = timepoint , y = value , color = obesity))+ scale_colour_manual(values=c("blue", "red")) + #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_blank() , axis.title.y = element_blank() , strip.text.x = element_text(size = 20, margin = margin(2,0,2,0, "cm")) , legend.position = "top" , legend.title = element_text(color = "black", size = 20) , legend.text = element_text(size = 15) , plot.subtitle = element_text(size = 20, hjust = 0.5) , plot.title = element_text(size = 20, hjust = 0.5)) + labs(title ="NPA")+ facet_wrap(~mediator, scales = "free", nrow = 5, ncol = 7) #p3 p4 = p3 + stat_pvalue_manual(stat_npa3_v2 , y.position = "my_y_pos" , step.increase = 0.08 #, y.position = 50 , label = "{p_format} {p.signif}" , hide.ns=T , tip.length = 0) #+ #scale_y_continuous(expand = expansion(mult = c(0.05, 0.25))) p4 dev.off()