From 8015dd4d6a4c45a8c871fa76934ebd8a8f2561f3 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Fri, 30 Oct 2020 20:56:15 +0000 Subject: [PATCH] fixed colours as well --- boxplot_with_stats.R | 101 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 92 insertions(+), 9 deletions(-) diff --git a/boxplot_with_stats.R b/boxplot_with_stats.R index 7736230..039534e 100755 --- a/boxplot_with_stats.R +++ b/boxplot_with_stats.R @@ -9,7 +9,7 @@ getwd() #============= # Input #============= -#source("boxplot_linear.R") +source("boxplot_linear.R") ####################################################### med_names = c("eotaxin3", "il12p70", "itac", "il13") lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,] @@ -73,7 +73,8 @@ stat_npa3 <- lf_test %>% group_by(timepoint, mediator) %>% wilcox_test(value ~ obesity, paired = F) %>% add_significance("p") -stat_npa3 +stat_npa3 +stat_npa3$p_format = round(stat_npa3$p, 3) stat_npa3 <- stat_npa3 %>% add_xy_position(x = "timepoint", dodge = 0.8) @@ -108,6 +109,8 @@ 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 @@ -115,24 +118,25 @@ p3 = ggplot(lf_test)+ geom_boxplot(aes(x = timepoint , vjust = 0) , axis.title.x = element_blank() , axis.title.y = element_blank() - , legend.position = "none" + , strip.text.x = element_text(size = 20) + , 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") - -p3_build = ggbuild(p3) - + 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} {p.signif}" + , 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 @@ -150,3 +154,82 @@ get_facet_ymax <- function (x){ } 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" + #, "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_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) + , 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" + #, 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()