From af334b559e86a5a2895e0f94281a40504a5d6f2d Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Fri, 30 Oct 2020 17:23:40 +0000 Subject: [PATCH] added boxplot_stat_function.R --- boxplot_stat_function.R | 95 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 95 insertions(+) create mode 100644 boxplot_stat_function.R diff --git a/boxplot_stat_function.R b/boxplot_stat_function.R new file mode 100644 index 0000000..43b994b --- /dev/null +++ b/boxplot_stat_function.R @@ -0,0 +1,95 @@ +#!/usr/bin/Rscript +getwd() +setwd("~/git/mosaic_2020/") +getwd() +############################################################ +# TASK: boxplots at T1 +# FIXME: currently not rendering, problem with NAs for stats? +############################################################ +doMyPlots <- function(x, title = "TEST") { + mediators = levels(as.factor(x$mediator)) + + plots <- list() + + for (i in mediators) { + cat("Plotting:", i, "\n") + single=x[x$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("#00BFC4", "#F8766D")) + + p2 = p2 +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_ + , 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 + , x = "" + , y = "Levels") + + # stats + 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) + + # add stats to boxplot + 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 + } + return(plots) + + } + +############################################################### +#============ +# test +#============ +source("boxplot_linear.R") + +#============= +# Output: +#============= +outfile_bp = paste0("boxplots_TEST_stats_all", ".pdf") +output_boxplot_stats = paste0(outdir_plots, outfile_bp); output_boxplot_stats + +pdf(output_boxplot_stats, width=32, height=18) + +med_names = c("eotaxin3", "il12p70", "itac", "il13") +#lf_test = lf_fp_npa[lf_fp_npa$mediator == "il12p70",] +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 + +# entire sample_type data +#bar=doMyPlots(lf_fp_npa) +#bang=cowplot::plot_grid(plotlist=bar, align = 'hv', ncol=7, nrow=5) +#bang + +dev.off()