added boxplot_stat_log.R to generate graphs with log data and stats

This commit is contained in:
Tanushree Tunstall 2020-11-02 14:28:52 +00:00
parent 550e8b2c50
commit b8ff7995ac
6 changed files with 249 additions and 10 deletions

View file

@ -28,7 +28,7 @@ fp_npa = length(unique(lf_fp_npa$mosaic)); fp_npa
cat("\nPlotting boxplots with stats for:", my_sample_npa cat("\nPlotting boxplots with stats for:", my_sample_npa
, "\n========================================================\n") , "\n========================================================\n")
plots_npa = doMyPlots(lf_fp_npa) plots_npa = doMyPlotsStats(lf_fp_npa)
npa_plot = ggpubr::ggarrange(plotlist = plots_npa npa_plot = ggpubr::ggarrange(plotlist = plots_npa
, align = "hv" , align = "hv"
, ncol = 7 , ncol = 7
@ -94,7 +94,7 @@ fp_serum = length(unique(lf_fp_serum$mosaic)); fp_serum
cat("\nPlotting boxplots with stats for:", my_sample_serum cat("\nPlotting boxplots with stats for:", my_sample_serum
, "\n========================================================\n") , "\n========================================================\n")
plots_serum = doMyPlots(lf_fp_serum) plots_serum = doMyPlotsStats(lf_fp_serum)
serum_plot = ggpubr::ggarrange(plotlist = plots_serum serum_plot = ggpubr::ggarrange(plotlist = plots_serum
, align = "hv" , align = "hv"
, ncol = 7 , ncol = 7

View file

@ -5,7 +5,7 @@ getwd()
############################################################ ############################################################
# TASK: boxplots with stats # TASK: boxplots with stats
############################################################ ############################################################
doMyPlots <- function(df) { doMyPlotsStats <- function(df) {
mediators = levels(as.factor(df$mediator)) mediators = levels(as.factor(df$mediator))
plots <- list() plots <- list()
@ -61,4 +61,5 @@ doMyPlots <- function(df) {
return(plots) return(plots)
} }
############################################################### ###############################################################

View file

@ -22,8 +22,16 @@ pdf(output_boxplot_stats, width=22, height=16)
#med_names = c("eotaxin3", "il12p70", "itac", "il13") #med_names = c("eotaxin3", "il12p70", "itac", "il13")
#lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,] #lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,]
#baz=cowplot::plot_grid(plotlist=foo, align = 'hv', ncol=2, nrow=2) baz=cowplot::plot_grid(plotlist=foo, align = 'hv', ncol=2, nrow=2)
foo_annot = annotate_figure(foo baz
sam_plot = ggpubr::ggarrange(plotlist = foo
, align = "hv"
, ncol = 7
, nrow = 5
, common.legend = T)
sam_plot
foo_annot = annotate_figure(sam_plot
, top = text_grob(my_sample_name , top = text_grob(my_sample_name
, color = "blue" , color = "blue"
, face = "bold" , face = "bold"
@ -41,3 +49,80 @@ foo_annot = annotate_figure(foo
foo_annot foo_annot
dev.off() dev.off()
################################################################## ##################################################################
# ind components of function
my_sample_name = "SAM"
to_remove = c("ifna2a"
, "il29")
lf_fp_sam2 = lf_fp_sam[!lf_fp_sam$mediator%in%to_remove,]
lf_fp_sam3 = lf_fp_sam[lf_fp_sam$mediator%in%to_remove,]
df_test = lf_fp_sam3
df_test = lf_fp_sam2
table(df_test$mediator)
p = ggplot(df_test)+ geom_boxplot(aes(x = timepoint
, y = value
, color = obesity
#, palette = c("#00BFC4", "#F8766D")
))+
scale_colour_manual(values=c("#00BFC4", "#F8766D")) +
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()
, legend.position = "none"
, plot.subtitle = element_text(size = 20, hjust = 0.5)
, plot.title = element_text(size = 20, hjust = 0.5)) +
labs(title = "TEST")
p
#--------
# stats
#---------
to_remove = c("ifna2a"
, "il29")
lf_fp_sam3 = lf_fp_sam[lf_fp_sam$mediator%in%to_remove,]
lf_fp_sam4 = lf_fp_sam3[lf_fp_sam3$timepoint == "t3",]
table(lf_fp_sam4$timepoint)
wilcox_test(value~obesity, data = lf_fp_sam4)
table(lf_fp_sam3$mediator, lf_fp_sam3$timepoint)
time = c("t3")
lf_fp_sam3 = lf_fp_sam3[!lf_fp_sam3$timepoint%in%time,]
table(lf_fp_sam3$mediator, lf_fp_sam3$timepoint)
med = c("il29")
time2 = c("t2")
lf_fp_sam3 = lf_fp_sam3[!lf_fp_sam3$mediator%in%med && lf_fp_sam3$timepoint%in%time2,]
table(lf_fp_sam3$mediator, lf_fp_sam3$timepoint)
df_test = lf_fp_sam3
table(df_test$mediator)
stat_df <- df_test %>%
group_by(timepoint, mediator) %>%
wilcox_test(value ~ obesity, paired = F) %>%
add_significance("p")
stat_df$p_format = round(stat_df$p, digits = 3)
stat_df <- stat_df %>%
add_xy_position(x = "timepoint", dodge = 0.8)
p2 = p + stat_pvalue_manual(stat_df
#, y.position = max_y
, label = "{p_format} {p.signif}"
, hide.ns = T
, tip.length = 0)+
scale_y_continuous(expand = expansion(mult = c(0.05, 0.25)))
p2

153
boxplot_stat_log.R Normal file
View file

@ -0,0 +1,153 @@
#!/usr/bin/Rscript
getwd()
setwd("~/git/mosaic_2020/")
getwd()
###############################################################
#=============
# Input
#=============
source("boxplot_stat_function.R")
source("plot_data.R")
#=============
# Output
#=============
outfile_bp = paste0("boxplots_stats_npa_serum_LOG", ".pdf")
output_boxplot_stats = paste0(outdir_plots, outfile_bp); output_boxplot_stats
pdf(output_boxplot_stats, width=22, height=16)
###############################################################
#===============================
# data assignment for plots
#===============================
#----------------------------------------
#-------
# NPA
#-------
# vl data juggling with 0!
vl_temp_fix = lf_fp_npa[lf_fp_npa$mediator=="vl_pfu_ul",]
vl_temp_fix$value[1:20]
vl_temp_fix$value[vl_temp_fix$value == 0] <- 1
vl_temp_fix$value[1:20]
orig_nrow = nrow(lf_fp_npa); orig_nrow; table(lf_fp_npa$mediator)
lf_fp_npa = lf_fp_npa[!lf_fp_npa$mediator=="vl_pfu_ul",]
temp_nrow = nrow(lf_fp_npa); temp_nrow; table(lf_fp_npa$mediator)
# add vl_temp_fix back
lf_fp_npa = rbind(lf_fp_npa, vl_temp_fix)
revised_nrow = nrow(lf_fp_npa); revised_nrow
orig_nrow == revised_nrow
head(lf_fp_npa$value)
# log all values and reassign
lf_fp_npa$value = log10(lf_fp_npa$value)
head(lf_fp_npa$value)
vl_temp_fix2 = lf_fp_npa[lf_fp_npa$mediator=="vl_pfu_ul",]
vl_temp_fix2$value[1:20]
#------------------------------------------
my_sample_npa = "NPA"
fp_npa = length(unique(lf_fp_npa$mosaic)); fp_npa
cat("\nPlotting boxplots with stats for:", my_sample_npa
, "\n========================================================\n")
plots_npa = doMyPlotsStats(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_npa
, color = "blue"
, face = "bold"
, size = 14)
, bottom = text_grob(paste0("Mosaic data\nFlu positive adults (n=", fp_npa, ")")
, color = "blue"
, hjust = 1
, x = 0.98
, face = "italic"
, size = 10)
, left = text_grob("Levels (pg/ml)"
, color = "black"
, rot = 90
, size = 18))
npa_plot_annot
#dev.off()
#-------------
# SAM
#-------------
#my_sample_sam = "SAM"
#fp_sam = length(unique(lf_fp_sam$mosaic)); fp_sam
#cat("\nPlotting boxplots with stats for:", my_sample_sam
# , "\n========================================================\n")
#plots_sam = doMyPlots(lf_fp_sam)
#sam_plot = ggpubr::ggarrange(plotlist = plots_sam
# , align = "hv"
# , ncol = 7
# , nrow = 5
# , common.legend = T)
#sam_plot
#sam_plot_annot = annotate_figure(sam_plot
# , top = text_grob(my_sample_sam
# , color = "blue"
# , face = "bold"
# , size = 14)
# , bottom = text_grob(paste0("Mosaic data\nFlu positive adults (n=", fp_sam, ")")
# , color = "blue"
# , hjust = 1
# , x = 0.98
# , face = "italic"
# , size = 10)
# , left = text_grob("Levels (pg/ml)"
# , color = "black"
# , rot = 90
# , size = 18))
#sam_plot_annot
#dev.off()
#------------------------------------------
#-------------
# SERUM
#-------------
head(lf_fp_serum$value)
# reassign
lf_fp_serum$value = log10(lf_fp_serum$value)
head(lf_fp_serum$value)
my_sample_serum = "serum"
fp_serum = length(unique(lf_fp_serum$mosaic)); fp_serum
cat("\nPlotting boxplots with stats for:", my_sample_serum
, "\n========================================================\n")
plots_serum = doMyPlotsStats(lf_fp_serum)
serum_plot = ggpubr::ggarrange(plotlist = plots_serum
, align = "hv"
, ncol = 7
, nrow = 5
, common.legend = T)
#serum_plot
serum_plot_annot = annotate_figure(serum_plot
, top = text_grob(my_sample_serum
, color = "blue"
, face = "bold"
, size = 14)
, bottom = text_grob(paste0("Mosaic data\nFlu positive adults (n=", fp_serum, ")")
, color = "blue"
, hjust = 1
, x = 0.98
, face = "italic"
, size = 10)
, left = text_grob("Levels (pg/ml)"
, color = "black"
, rot = 90
, size = 18))
serum_plot_annot
dev.off()
###############################################################

View file

@ -9,7 +9,7 @@ getwd()
#============= #=============
# Input # Input
#============= #=============
source("boxplot_linear.R") source("plot_data.R")
####################################################### #######################################################
med_names = c("eotaxin3", "il12p70", "itac", "il13") med_names = c("eotaxin3", "il12p70", "itac", "il13")
lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,] lf_test = lf_fp_npa[lf_fp_npa$mediator%in%med_names,]

View file

@ -21,7 +21,7 @@ metadata_all$mosaic[metadata_all$adult==1 & metadata_all$age<=18]
#----------- #-----------
# npa # npa
#----------- #-----------
#wf_fp_npa = npa_wf[npa_wf$flustat == 1,] wf_fp_npa = npa_wf[npa_wf$flustat == 1,]
lf_fp_npa = npa_lf[npa_lf$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 = paste0("t", lf_fp_npa$timepoint)
lf_fp_npa$timepoint = as.factor(lf_fp_npa$timepoint) lf_fp_npa$timepoint = as.factor(lf_fp_npa$timepoint)
@ -35,7 +35,7 @@ table(lf_fp_npa$mediator)
#----------- #-----------
# sam # sam
#----------- #-----------
#wf_fp_sam = samm_wf[samm_wf$flustat == 1,] wf_fp_sam = sam_wf[sam_wf$flustat == 1,]
lf_fp_sam = sam_lf[sam_lf$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 = paste0("t", lf_fp_sam$timepoint)
lf_fp_sam$timepoint = as.factor(lf_fp_sam$timepoint) lf_fp_sam$timepoint = as.factor(lf_fp_sam$timepoint)
@ -49,7 +49,7 @@ table(lf_fp_sam$mediator)
#----------- #-----------
# serum # serum
#----------- #-----------
#wf_fp_serum = serum_wf[serum_wf$flustat == 1,] wf_fp_serum = serum_wf[serum_wf$flustat == 1,]
lf_fp_serum = serum_lf[serum_lf$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 = paste0("t", lf_fp_serum$timepoint)
lf_fp_serum$timepoint = as.factor(lf_fp_serum$timepoint) lf_fp_serum$timepoint = as.factor(lf_fp_serum$timepoint)