From 37f8cd619eb46a06c104f891048506b43cc2c326 Mon Sep 17 00:00:00 2001
From: Tanushree Tunstall <tanu@tunstall.in>
Date: Fri, 27 Nov 2020 18:22:00 +0000
Subject: [PATCH] added stats for non_severe flu stats in resp dir

---
 non_severe/flu_ns_stats_unpaired_npa.R   | 395 ++++++++++++++++++++++
 non_severe/flu_ns_stats_unpaired_sam.R   | 397 +++++++++++++++++++++++
 non_severe/flu_ns_stats_unpaired_serum.R | 392 ++++++++++++++++++++++
 3 files changed, 1184 insertions(+)
 create mode 100755 non_severe/flu_ns_stats_unpaired_npa.R
 create mode 100755 non_severe/flu_ns_stats_unpaired_sam.R
 create mode 100755 non_severe/flu_ns_stats_unpaired_serum.R

diff --git a/non_severe/flu_ns_stats_unpaired_npa.R b/non_severe/flu_ns_stats_unpaired_npa.R
new file mode 100755
index 0000000..7535923
--- /dev/null
+++ b/non_severe/flu_ns_stats_unpaired_npa.R
@@ -0,0 +1,395 @@
+#!/usr/bin/Rscript   
+getwd()
+setwd("~/git/mosaic_2020/")
+getwd()
+############################################################
+# TASK: unpaired (time) analysis of mediators:
+# sample type: NPA
+# data: Flu positive adult patients
+# group: obesity
+############################################################
+my_sample_type = "npa"
+
+#=============
+# Input
+#=============
+source("data_extraction_mediators.R")
+
+# check: copd and asthma conflict
+table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
+
+#=============
+# Output
+#=============
+outfile_name = paste0("flu_stats_time_unpaired_ns_", my_sample_type, ".csv")
+flu_stats_time_unpaired = paste0(outdir_stats_ns, outfile_name)
+flu_stats_time_unpaired
+
+# quick checks
+table(npa_wf$T1_resp_score)
+table(npa_lf$mediator, npa_lf$timepoint, npa_lf$T1_resp_score)
+table(npa_lf$T1_resp_score!=3)
+
+#===============================
+# Data assignment for stats: fp, non_severe patients
+#================================
+wf = npa_wf[npa_wf$T1_resp_score != 3,]
+lf = npa_lf[npa_lf$T1_resp_score != 3,]
+lf$timepoint = paste0("t", lf$timepoint)
+lf = lf[!lf$mediator == "vitd",]
+
+########################################################################
+# clear variables
+rm(sam_lf, sam_wf
+, serum_lf, serum_wf) 
+rm(colnames_sam_df, expected_rows_sam_lf
+, colnames_serum_df, expected_rows_serum_lf)
+
+rm(pivot_cols)
+
+# sanity checks
+table(lf$timepoint)
+length(unique(lf$mosaic))
+
+#if (table(lf$flustat) == table(npa_lf$flustat)[[2]]){
+#  cat("Analysing Flu positive patients for:", my_sample_type)
+#}else{
+#  cat("FAIL: problem with subsetting data for:", my_sample_type)
+#  quit()
+#}
+########################################################################
+# Unpaired stats at each timepoint b/w groups: wilcoxon UNpaired analysis 
+# with correction
+########################################################################
+# with adjustment: fdr and BH are identical
+my_adjust_method = "BH"
+
+#==============
+# unpaired: t1
+#==============
+lf_t1 = lf[lf$timepoint == "t1",]
+sum(is.na(lf_t1$value))
+
+foo = lf_t1[which(is.na(lf_t1$value)),]
+#ci = which(is.na(lf_t1$value))
+#lf_t1_comp = lf_t1[-ci,]
+
+lf_t1_comp = lf_t1[-which(is.na(lf_t1$value)),]
+stats_un_t1 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t1 
+                            , data = lf_t1_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+
+foo$mosaic[!unique(foo$mosaic)%in%unique(lf_t1_comp$mosaic)]
+
+# add timepoint and convert to df
+stats_un_t1$timepoint = "t1"
+stats_un_t1 = as.data.frame(stats_un_t1)
+class(stats_un_t1)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t1
+#----------------------------------------
+#n_t1 = data.frame(table(lf_t1_comp$mediator))
+n_t1_all = data.frame(table(lf_t1$mediator))
+colnames(n_t1_all) = c("mediator", "n_obs")
+n_t1_all$mediator = as.character(n_t1_all$mediator)
+
+n_t1_comp = data.frame(table(lf_t1_comp$mediator))
+colnames(n_t1_comp) = c("mediator", "n_obs_complete")
+n_t1_comp$mediator = as.character(n_t1_comp$mediator)
+
+merge_cols = intersect(names(n_t1_all), names(n_t1_comp)); merge_cols
+n_t1= merge(n_t1_all, n_t1_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#===================================
+merging_cols = intersect(names(stats_un_t1), names(n_t1)); merging_cols
+if (all(n_t1$mediator%in%stats_un_t1$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t1)
+      , "\nncol:", ncol(stats_un_t1))
+}else{
+  nf = n_t1$mediator[!n_t1$mediator%in%stats_un_t1$mediator]
+  stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t1)
+      , "\nncol:", ncol(stats_un_t1))
+}
+
+# add bonferroni adjustment as well
+stats_un_t1$p_adj_bonferroni = p.adjust(stats_un_t1$p, method = "bonferroni")
+
+rm(n_t1)
+rm(lf_t1_comp)
+
+#==============
+# unpaired: t2
+#==============
+lf_t2 = lf[lf$timepoint == "t2",]
+lf_t2_comp = lf_t2[-which(is.na(lf_t2$value)),]
+                   
+stats_un_t2 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t2 
+                            , data = lf_t2_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+
+# add timepoint and convert to df
+stats_un_t2$timepoint = "t2"
+stats_un_t2 = as.data.frame(stats_un_t2)
+class(stats_un_t2)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t2
+#----------------------------------------
+#n_t2 = data.frame(table(lf_t2_comp$mediator))
+n_t2_all = data.frame(table(lf_t2$mediator))
+colnames(n_t2_all) = c("mediator", "n_obs")
+n_t2_all$mediator = as.character(n_t2_all$mediator)
+
+n_t2_comp = data.frame(table(lf_t2_comp$mediator))
+colnames(n_t2_comp) = c("mediator", "n_obs_complete")
+n_t2_comp$mediator = as.character(n_t2_comp$mediator)
+
+merge_cols = intersect(names(n_t2_all), names(n_t2_comp)); merge_cols
+n_t2= merge(n_t2_all, n_t2_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
+merging_cols = intersect(names(stats_un_t2), names(n_t2)); merging_cols
+if (all(n_t2$mediator%in%stats_un_t2$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t2)
+      , "\nncol:", ncol(stats_un_t2))
+}else{
+  nf = n_t2$mediator[!n_t2$mediator%in%stats_un_t2$mediator]
+  stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t2)
+      , "\nncol:", ncol(stats_un_t2))
+}
+
+# add bonferroni adjustment as well
+stats_un_t2$p_adj_bonferroni = p.adjust(stats_un_t2$p, method = "bonferroni")
+
+rm(n_t2)
+rm(lf_t2_comp)
+
+#==============
+# unpaired: t3
+#==============
+lf_t3 = lf[lf$timepoint == "t3",]
+lf_t3_comp = lf_t3[-which(is.na(lf_t3$value)),]
+
+stats_un_t3 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t3 
+                            , data = lf_t3_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+
+# add timepoint and convert to df
+stats_un_t3$timepoint = "t3"
+stats_un_t3 = as.data.frame(stats_un_t3)
+class(stats_un_t3)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t3
+#----------------------------------------
+#n_t3 = data.frame(table(lf_t3_comp$mediator))
+n_t3_all = data.frame(table(lf_t3$mediator))
+colnames(n_t3_all) = c("mediator", "n_obs")
+n_t3_all$mediator = as.character(n_t3_all$mediator)
+
+n_t3_comp = data.frame(table(lf_t3_comp$mediator))
+colnames(n_t3_comp) = c("mediator", "n_obs_complete")
+n_t3_comp$mediator = as.character(n_t3_comp$mediator)
+
+merge_cols = intersect(names(n_t3_all), names(n_t3_comp)); merge_cols
+n_t3= merge(n_t3_all, n_t3_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
+merging_cols = intersect(names(stats_un_t3), names(n_t3)); merging_cols
+if (all(n_t3$mediator%in%stats_un_t3$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t3)
+      , "\nncol:", ncol(stats_un_t3))
+}else{
+  nf = n_t3$mediator[!n_t3$mediator%in%stats_un_t3$mediator]
+  stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t3)
+      , "\nncol:", ncol(stats_un_t3))
+}
+
+# add bonferroni adjustment as well
+stats_un_t3$p_adj_bonferroni = p.adjust(stats_un_t3$p, method = "bonferroni")
+
+rm(n_t3)
+rm(lf_t3_comp)
+
+########################################################################
+#==============
+# Rbind these dfs
+#==============
+str(stats_un_t1);str(stats_un_t2); str(stats_un_t3)
+
+n_dfs = 3
+
+if ( all.equal(nrow(stats_un_t1), nrow(stats_un_t2), nrow(stats_un_t3)) && 
+     all.equal(ncol(stats_un_t1), ncol(stats_un_t2), ncol(stats_un_t3)) ) {
+  expected_rows = nrow(stats_un_t1) * n_dfs
+  expected_cols = ncol(stats_un_t1)
+  print("PASS: expected_rows and cols variables generated for downstream sanity checks")
+}else{
+  cat("FAIL: dfs have different no. of rows and cols"
+      , "\nCheck harcoded value of n_dfs"
+      , "\nexpected_rows and cols could not be generated")
+  quit()
+}
+
+if ( all.equal(colnames(stats_un_t1), colnames(stats_un_t2), colnames(stats_un_t3)) ){
+  print("PASS: colnames match. Rbind the 3 dfs...")
+  combined_unpaired_stats = rbind(stats_un_t1, stats_un_t2, stats_un_t3)
+} else{
+  cat("FAIL: cannot combined dfs. Colnames don't match!")
+  quit()
+}
+
+if ( nrow(combined_unpaired_stats) == expected_rows && ncol(combined_unpaired_stats) == expected_cols ){
+  cat("PASS: combined_df has expected dimension"
+      , "\nNo. of rows in combined_df:", nrow(combined_unpaired_stats)
+      , "\nNo. of cols in combined_df:", ncol(combined_unpaired_stats) )
+}else{
+  cat("FAIL: combined_df dimension mismatch")
+  quit()
+}
+
+#######################################################################
+#=================
+# formatting df
+#=================
+# delete: unnecessary column
+combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.))
+
+# add sample_type
+cat("Adding sample type info as a column", my_sample_type, "...")
+combined_unpaired_stats$sample_type = my_sample_type 
+
+# add: reflect stats method correctly i.e paired or unpaired
+# incase there are NA due to LLODs, the gsub won't work!
+#combined_unpaired_stats$method = gsub("Wilcoxon", "Wilcoxon_unpaired", combined_unpaired_stats$method)
+combined_unpaired_stats$method = "wilcoxon unpaired"
+combined_unpaired_stats$method
+
+# add an extra column for padjust_signif: my_adjust_method
+combined_unpaired_stats$padjust_signif = combined_unpaired_stats$p.adj
+# add appropriate symbols for padjust_signif: my_adjust_method
+combined_unpaired_stats  = dplyr::mutate(combined_unpaired_stats, padjust_signif = case_when(padjust_signif == 0.05 ~ "."
+                                                                                             , padjust_signif <=0.0001 ~ '****'
+                                                                                             , padjust_signif <=0.001 ~ '***'
+                                                                                             , padjust_signif <=0.01 ~ '**'
+                                                                                             , padjust_signif <0.05 ~ '*'
+                                                                                             , TRUE ~ 'ns'))
+# add an extra column for p_bon_signif
+combined_unpaired_stats$p_bon_signif = combined_unpaired_stats$p_adj_bonferroni
+# add appropriate symbols for p_bon_signif
+combined_unpaired_stats  = dplyr::mutate(combined_unpaired_stats, p_bon_signif = case_when(p_bon_signif == 0.05 ~ "."
+                                                                                             , p_bon_signif <=0.0001 ~ '****'
+                                                                                             , p_bon_signif <=0.001 ~ '***'
+                                                                                             , p_bon_signif <=0.01 ~ '**'
+                                                                                             , p_bon_signif <0.05 ~ '*'
+                                                                                             , TRUE ~ 'ns'))
+# reorder columns
+print("preparing to reorder columns...")
+colnames(combined_unpaired_stats)
+my_col_order2 = c("mediator"
+                  , "timepoint"
+                  , "sample_type"
+                  , "n_obs"
+                  , "n_obs_complete"
+                  , "group1"
+                  , "group2"
+                  , "method"
+                  , "p"
+                  , "p.format"
+                  , "p.signif"
+                  , "p.adj"
+                  , "padjust_signif"
+                  , "p_adj_bonferroni"
+                  , "p_bon_signif")
+
+if( length(my_col_order2) == ncol(combined_unpaired_stats) && (all(my_col_order2%in%colnames(combined_unpaired_stats))) ){
+  print("PASS: Reordering columns...")
+  combined_unpaired_stats_f = combined_unpaired_stats[, my_col_order2]
+  print("Successful: column reordering")
+  print("formatted df called:'combined_unpaired_stats_f'")
+  cat('\nformatted df has the following dimensions\n')
+  print(dim(combined_unpaired_stats_f ))
+} else{
+  cat(paste0("FAIL:Cannot reorder columns, length mismatch"
+             , "\nExpected column order for: ", ncol(combined_unpaired_stats)
+             , "\nGot:", length(my_col_order2)))
+  quit()
+}  
+# assign nice column names like replace "." with "_"
+colnames(combined_unpaired_stats_f) = c("mediator"
+                                        , "timepoint" 
+                                        , "sample_type"
+                                        , "n_obs" 
+                                        , "n_obs_complete"
+                                        , "group1"
+                                        , "group2"
+                                        , "method" 
+                                        , "p" 
+                                        , "p_format"
+                                        , "p_signif"
+                                        , paste0("p_adj_fdr_", my_adjust_method)
+                                        , paste0("p_", my_adjust_method, "_signif")
+                                        , "p_adj_bonferroni"
+                                        , "p_bon_signif")
+
+colnames(combined_unpaired_stats_f)
+
+#---------------
+# quick summary
+#---------------
+# count how many meds are significant
+n_sig  = length(combined_unpaired_stats_f$mediator[combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)])
+sig_meds = combined_unpaired_stats_f[(combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)),]
+
+sig_meds$med_time = paste0(sig_meds$mediator, "@", sig_meds$timepoint)
+
+cat("\nTotal no. of statistically significant mediators in", toupper(my_sample_type)
+    , "are:", n_sig
+    , "\nThese are:", sig_meds$med_time)
+
+########################################################################
+#******************
+# write output file
+#******************
+cat("\nUNpaired stats for groups will be:", flu_stats_time_unpaired)
+write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
diff --git a/non_severe/flu_ns_stats_unpaired_sam.R b/non_severe/flu_ns_stats_unpaired_sam.R
new file mode 100755
index 0000000..902ea9f
--- /dev/null
+++ b/non_severe/flu_ns_stats_unpaired_sam.R
@@ -0,0 +1,397 @@
+#!/usr/bin/Rscript   
+getwd()
+setwd("~/git/mosaic_2020/")
+getwd()
+############################################################
+# TASK: unpaired (time) analysis of mediators:
+# sample type: SAM
+# data: Flu positive adult patients
+# group: obesity
+############################################################
+my_sample_type = "sam"
+
+#=============
+# Input
+#=============
+source("data_extraction_mediators.R")
+
+# check: copd and asthma conflict
+table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
+
+#=============
+# Output
+#=============
+outfile_name = paste0("flu_stats_time_unpaired_ns_", my_sample_type, ".csv")
+flu_stats_time_unpaired = paste0(outdir_stats_ns, outfile_name)
+flu_stats_time_unpaired
+
+# quick checks
+table(sam_wf$T1_resp_score)
+table(sam_lf$mediator, sam_lf$timepoint, sam_lf$T1_resp_score)
+table(sam_lf$T1_resp_score!=3)
+
+#===============================
+# Data assignment for stats: fp, non_severe patients
+#================================
+wf = sam_wf[sam_wf$T1_resp_score != 3,]
+lf = sam_lf[sam_lf$T1_resp_score != 3,]
+lf$timepoint = paste0("t", lf$timepoint)
+lf = lf[!lf$mediator == "vitd",]
+
+########################################################################
+# clear variables
+rm(npa_lf, npa_wf
+, serum_lf, serum_wf) 
+rm(colnames_npa_df, expected_rows_npa_lf
+, colnames_serum_df, expected_rows_serum_lf) 
+
+rm(pivot_cols)
+
+# sanity checks
+table(lf$timepoint)
+length(unique(lf$mosaic))
+
+#if (table(lf$flustat) == table(sam_lf$flustat)[[2]]){
+#  cat("Analysing Flu positive patients for:", my_sample_type)
+#}else{
+#  cat("FAIL: problem with subsetting data for:", my_sample_type)
+#  quit()
+#}
+
+########################################################################
+# Unpaired stats at each timepoint b/w groups: wilcoxon UNpaired analysis 
+# with correction
+#######################################################################
+# with adjustment: fdr and BH are identical
+my_adjust_method = "BH"
+
+#==============
+# unpaired: t1
+#==============
+lf_t1 = lf[lf$timepoint == "t1",]
+sum(is.na(lf_t1$value))
+
+foo = lf_t1[which(is.na(lf_t1$value)),]
+#ci = which(is.na(lf_t1$value))
+#lf_t1_comp = lf_t1[-ci,]
+
+lf_t1_comp = lf_t1[-which(is.na(lf_t1$value)),]
+stats_un_t1 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t1 
+                            , data = lf_t1_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+
+foo$mosaic[!unique(foo$mosaic)%in%unique(lf_t1_comp$mosaic)]
+
+# add timepoint and convert to df
+stats_un_t1$timepoint = "t1"
+stats_un_t1 = as.data.frame(stats_un_t1)
+class(stats_un_t1)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t1
+#----------------------------------------
+#n_t1 = data.frame(table(lf_t1_comp$mediator))
+n_t1_all = data.frame(table(lf_t1$mediator))
+colnames(n_t1_all) = c("mediator", "n_obs")
+n_t1_all$mediator = as.character(n_t1_all$mediator)
+
+n_t1_comp = data.frame(table(lf_t1_comp$mediator))
+colnames(n_t1_comp) = c("mediator", "n_obs_complete")
+n_t1_comp$mediator = as.character(n_t1_comp$mediator)
+
+merge_cols = intersect(names(n_t1_all), names(n_t1_comp)); merge_cols
+n_t1= merge(n_t1_all, n_t1_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#===================================
+merging_cols = intersect(names(stats_un_t1), names(n_t1)); merging_cols
+if (all(n_t1$mediator%in%stats_un_t1$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t1)
+      , "\nncol:", ncol(stats_un_t1))
+}else{
+  nf = n_t1$mediator[!n_t1$mediator%in%stats_un_t1$mediator]
+  stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t1)
+      , "\nncol:", ncol(stats_un_t1))
+}
+
+# add bonferroni adjustment as well
+stats_un_t1$p_adj_bonferroni = p.adjust(stats_un_t1$p, method = "bonferroni")
+
+rm(n_t1)
+rm(lf_t1_comp)
+
+#==============
+# unpaired: t2
+#==============
+lf_t2 = lf[lf$timepoint == "t2",]
+lf_t2_comp = lf_t2[-which(is.na(lf_t2$value)),]
+                   
+stats_un_t2 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t2 
+                            , data = lf_t2_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+# add timepoint and convert to df                           
+stats_un_t2$timepoint = "t2"
+stats_un_t2 = as.data.frame(stats_un_t2)
+class(stats_un_t2)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t2
+#----------------------------------------
+#n_t2 = data.frame(table(lf_t2_comp$mediator))
+n_t2_all = data.frame(table(lf_t2$mediator))
+colnames(n_t2_all) = c("mediator", "n_obs")
+n_t2_all$mediator = as.character(n_t2_all$mediator)
+
+n_t2_comp = data.frame(table(lf_t2_comp$mediator))
+colnames(n_t2_comp) = c("mediator", "n_obs_complete")
+n_t2_comp$mediator = as.character(n_t2_comp$mediator)
+
+merge_cols = intersect(names(n_t2_all), names(n_t2_comp)); merge_cols
+n_t2= merge(n_t2_all, n_t2_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#===================================
+merging_cols = intersect(names(stats_un_t2), names(n_t2)); merging_cols
+if (all(n_t2$mediator%in%stats_un_t2$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t2)
+      , "\nncol:", ncol(stats_un_t2))
+}else{
+  nf = n_t2$mediator[!n_t2$mediator%in%stats_un_t2$mediator]
+  stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t2)
+      , "\nncol:", ncol(stats_un_t2))
+}
+
+# add bonferroni adjustment as well
+stats_un_t2$p_adj_bonferroni = p.adjust(stats_un_t2$p, method = "bonferroni")
+
+rm(n_t2)
+rm(lf_t2_comp)
+
+#==============
+# unpaired: t3
+#==============
+lf_t3 = lf[lf$timepoint == "t3",]
+lf_t3_comp = lf_t3[-which(is.na(lf_t3$value)),]
+
+stats_un_t3 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t3 
+                            , data = lf_t3_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+# add timepoint and convert to df                            
+stats_un_t3$timepoint = "t3"
+stats_un_t3 = as.data.frame(stats_un_t3)
+class(stats_un_t3)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t3
+#----------------------------------------
+#n_t3 = data.frame(table(lf_t3_comp$mediator))
+n_t3_all = data.frame(table(lf_t3$mediator))
+colnames(n_t3_all) = c("mediator", "n_obs")
+n_t3_all$mediator = as.character(n_t3_all$mediator)
+
+n_t3_comp = data.frame(table(lf_t3_comp$mediator))
+colnames(n_t3_comp) = c("mediator", "n_obs_complete")
+n_t3_comp$mediator = as.character(n_t3_comp$mediator)
+
+merge_cols = intersect(names(n_t3_all), names(n_t3_comp)); merge_cols
+n_t3= merge(n_t3_all, n_t3_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
+merging_cols = intersect(names(stats_un_t3), names(n_t3)); merging_cols
+if (all(n_t3$mediator%in%stats_un_t3$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t3)
+      , "\nncol:", ncol(stats_un_t3))
+}else{
+  nf = n_t3$mediator[!n_t3$mediator%in%stats_un_t3$mediator]
+  stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t3)
+      , "\nncol:", ncol(stats_un_t3))
+}
+
+# check: satisfied!!!!
+# FIXME: supply the col name automatically?
+wilcox.test(wf$ifna2a_sam3[wf$obesity == 1], wf$ifna2a_sam3[wf$obesity == 0])
+
+# add bonferroni adjustment as well
+stats_un_t3$p_adj_bonferroni = p.adjust(stats_un_t3$p, method = "bonferroni")
+
+rm(n_t3)
+rm(lf_t3_comp)
+########################################################################
+#=================
+# Rbind these dfs
+#=================
+str(stats_un_t1);str(stats_un_t2); str(stats_un_t3)
+
+n_dfs = 3
+
+if ( all.equal(nrow(stats_un_t1), nrow(stats_un_t2), nrow(stats_un_t3)) && 
+     all.equal(ncol(stats_un_t1), ncol(stats_un_t2), ncol(stats_un_t3)) ) {
+  expected_rows = nrow(stats_un_t1) * n_dfs
+  expected_cols = ncol(stats_un_t1)
+  print("PASS: expected_rows and cols variables generated for downstream sanity checks")
+}else{
+  cat("FAIL: dfs have different no. of rows and cols"
+      , "\nCheck harcoded value of n_dfs"
+      , "\nexpected_rows and cols could not be generated")
+  quit()
+}
+
+if ( all.equal(colnames(stats_un_t1), colnames(stats_un_t2), colnames(stats_un_t3)) ){
+  print("PASS: colnames match. Rbind the 3 dfs...")
+  combined_unpaired_stats = rbind(stats_un_t1, stats_un_t2, stats_un_t3)
+} else{
+  cat("FAIL: cannot combined dfs. Colnames don't match!")
+  quit()
+}
+
+if ( nrow(combined_unpaired_stats) == expected_rows && ncol(combined_unpaired_stats) == expected_cols ){
+  cat("PASS: combined_df has expected dimension"
+      , "\nNo. of rows in combined_df:", nrow(combined_unpaired_stats)
+      , "\nNo. of cols in combined_df:", ncol(combined_unpaired_stats) )
+}else{
+  cat("FAIL: combined_df dimension mismatch")
+  quit()
+}
+
+#######################################################################
+#=================
+# formatting df
+#=================
+# delete: unnecessary column
+combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.))
+
+# add sample_type
+cat("Adding sample type info as a column", my_sample_type, "...")
+combined_unpaired_stats$sample_type = my_sample_type 
+
+# add: reflect stats method correctly i.e paired or unpaired
+# incase there are NA due to LLODs, the gsub won't work!
+#combined_unpaired_stats$method = gsub("Wilcoxon", "Wilcoxon_unpaired", combined_unpaired_stats$method)
+combined_unpaired_stats$method = "wilcoxon unpaired"
+combined_unpaired_stats$method
+
+# add an extra column for padjust_signif: my_adjust_method
+combined_unpaired_stats$padjust_signif = combined_unpaired_stats$p.adj
+# add appropriate symbols for padjust_signif: my_adjust_method
+combined_unpaired_stats  = dplyr::mutate(combined_unpaired_stats, padjust_signif = case_when(padjust_signif == 0.05 ~ "."
+                                                                                             , padjust_signif <=0.0001 ~ '****'
+                                                                                             , padjust_signif <=0.001 ~ '***'
+                                                                                             , padjust_signif <=0.01 ~ '**'
+                                                                                             , padjust_signif <0.05 ~ '*'
+                                                                                             , TRUE ~ 'ns'))
+# add an extra column for p_bon_signif
+combined_unpaired_stats$p_bon_signif = combined_unpaired_stats$p_adj_bonferroni
+# add appropriate symbols for p_bon_signif
+combined_unpaired_stats  = dplyr::mutate(combined_unpaired_stats, p_bon_signif = case_when(p_bon_signif == 0.05 ~ "."
+                                                                                             , p_bon_signif <=0.0001 ~ '****'
+                                                                                             , p_bon_signif <=0.001 ~ '***'
+                                                                                             , p_bon_signif <=0.01 ~ '**'
+                                                                                             , p_bon_signif <0.05 ~ '*'
+                                                                                             , TRUE ~ 'ns'))
+# reorder columns
+print("preparing to reorder columns...")
+colnames(combined_unpaired_stats)
+my_col_order2 = c("mediator"
+                  , "timepoint"
+                  , "sample_type"
+                  , "n_obs"
+                  , "n_obs_complete"
+                  , "group1"
+                  , "group2"
+                  , "method"
+                  , "p"
+                  , "p.format"
+                  , "p.signif"
+                  , "p.adj"
+                  , "padjust_signif"
+                  , "p_adj_bonferroni"
+                  , "p_bon_signif")
+
+if( length(my_col_order2) == ncol(combined_unpaired_stats) && (all(my_col_order2%in%colnames(combined_unpaired_stats))) ){
+  print("PASS: Reordering columns...")
+  combined_unpaired_stats_f = combined_unpaired_stats[, my_col_order2]
+  print("Successful: column reordering")
+  print("formatted df called:'combined_unpaired_stats_f'")
+  cat('\nformatted df has the following dimensions\n')
+  print(dim(combined_unpaired_stats_f ))
+} else{
+  cat(paste0("FAIL:Cannot reorder columns, length mismatch"
+             , "\nExpected column order for: ", ncol(combined_unpaired_stats)
+             , "\nGot:", length(my_col_order2)))
+  quit()
+}  
+# assign nice column names like replace "." with "_"
+colnames(combined_unpaired_stats_f) = c("mediator"
+                                        , "timepoint" 
+                                        , "sample_type"
+                                        , "n_obs" 
+                                        , "n_obs_complete"
+                                        , "group1"
+                                        , "group2"
+                                        , "method" 
+                                        , "p" 
+                                        , "p_format"
+                                        , "p_signif"
+                                        , paste0("p_adj_fdr_", my_adjust_method)
+                                        , paste0("p_", my_adjust_method, "_signif")
+                                        , "p_adj_bonferroni"
+                                        , "p_bon_signif")
+
+colnames(combined_unpaired_stats_f)
+
+#---------------
+# quick summary
+#---------------
+# count how many meds are significant
+n_sig  = length(combined_unpaired_stats_f$mediator[combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)])
+sig_meds = combined_unpaired_stats_f[(combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)),]
+
+sig_meds$med_time = paste0(sig_meds$mediator, "@", sig_meds$timepoint)
+
+cat("\nTotal no. of statistically significant mediators in", toupper(my_sample_type)
+    , "are:", n_sig
+    , "\nThese are:", sig_meds$med_time)
+
+#######################################################################
+#******************
+# write output file
+#******************
+cat("\nUNpaired stats for groups will be:", flu_stats_time_unpaired)
+#write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
diff --git a/non_severe/flu_ns_stats_unpaired_serum.R b/non_severe/flu_ns_stats_unpaired_serum.R
new file mode 100755
index 0000000..a9dd416
--- /dev/null
+++ b/non_severe/flu_ns_stats_unpaired_serum.R
@@ -0,0 +1,392 @@
+#!/usr/bin/Rscript   
+getwd()
+setwd("~/git/mosaic_2020/")
+getwd()
+############################################################
+# TASK: unpaired (time) analysis of mediators: serum
+############################################################
+my_sample_type = "serum"
+
+#=============
+# Input
+#=============
+source("data_extraction_mediators.R")
+
+# check: copd and asthma conflict
+table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
+
+#=============
+# Output
+#=============
+outfile_name = paste0("flu_stats_time_unpaired_ns_", my_sample_type, ".csv")
+flu_stats_time_unpaired = paste0(outdir_stats_ns, outfile_name)
+flu_stats_time_unpaired
+
+# quick checks
+table(serum_wf$T1_resp_score)
+table(serum_lf$mediator,serum_lf$timepoint, serum_lf$T1_resp_score)
+table(serum_lf$T1_resp_score!=3)
+
+#===============================
+# data assignment for stats
+#================================
+wf = serum_wf[serum_wf$T1_resp_score != 3,]
+lf = serum_lf[serum_lf$T1_resp_score != 3,]
+lf$timepoint = paste0("t", lf$timepoint)
+
+########################################################################
+# clear variables
+rm(sam_lf, sam_wf
+, npa_lf, npa_wf) 
+rm(colnames_sam_df, expected_rows_sam_lf
+, colnames_npa_df, expected_rows_npa_lf)
+
+rm(pivot_cols)
+
+# sanity checks
+table(lf$timepoint)
+
+if (table(lf$flustat) == table(serum_lf$flustat)[[2]]){
+  cat("Analysing Flu positive patients for:", my_sample_type)
+}else{
+  cat("FAIL: problem with subsetting data for:", my_sample_type)
+  quit()
+}
+
+########################################################################
+# Unpaired stats at each timepoint b/w groups: wilcoxon UNpaired analysis
+# with correction
+#######################################################################
+# with adjustment: fdr and BH are identical
+my_adjust_method = "BH"
+
+#==============
+# unpaired: t1
+#==============
+lf_t1 = lf[lf$timepoint == "t1",]
+sum(is.na(lf_t1$value))
+
+foo = lf_t1[which(is.na(lf_t1$value)),]
+ci = which(is.na(lf_t1$value))
+#lf_t1_comp = lf_t1[-ci,]
+
+lf_t1_comp = lf_t1[-which(is.na(lf_t1$value)),]
+stats_un_t1 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t1 
+                            , data = lf_t1_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+
+foo$mosaic[!unique(foo$mosaic)%in%unique(lf_t1_comp$mosaic)]
+
+# add timepoint and convert to df
+stats_un_t1$timepoint = "t1"
+stats_un_t1 = as.data.frame(stats_un_t1)
+class(stats_un_t1)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t1
+#----------------------------------------
+#n_t1 = data.frame(table(lf_t1_comp$mediator))
+n_t1_all = data.frame(table(lf_t1$mediator))
+colnames(n_t1_all) = c("mediator", "n_obs")
+n_t1_all$mediator = as.character(n_t1_all$mediator)
+
+n_t1_comp = data.frame(table(lf_t1_comp$mediator))
+colnames(n_t1_comp) = c("mediator", "n_obs_complete")
+n_t1_comp$mediator = as.character(n_t1_comp$mediator)
+
+merge_cols = intersect(names(n_t1_all), names(n_t1_comp)); merge_cols
+n_t1= merge(n_t1_all, n_t1_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
+merging_cols = intersect(names(stats_un_t1), names(n_t1)); merging_cols
+if (all(n_t1$mediator%in%stats_un_t1$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t1)
+      , "\nncol:", ncol(stats_un_t1))
+}else{
+  nf = n_t1$mediator[!n_t1$mediator%in%stats_un_t1$mediator]
+  stats_un_t1 = merge(stats_un_t1, n_t1, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t1)
+      , "\nncol:", ncol(stats_un_t1))
+}
+
+# add bonferroni adjustment as well
+stats_un_t1$p_adj_bonferroni = p.adjust(stats_un_t1$p, method = "bonferroni")
+
+rm(n_t1)
+rm(lf_t1_comp)
+
+#==============
+# unpaired: t2
+#==============
+lf_t2 = lf[lf$timepoint == "t2",]
+lf_t2_comp = lf_t2[-which(is.na(lf_t2$value)),]
+                   
+stats_un_t2 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t2 
+                            , data = lf_t2_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+# add timepoint and convert to df
+stats_un_t2$timepoint = "t2"
+stats_un_t2 = as.data.frame(stats_un_t2)
+class(stats_un_t2)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t2
+#----------------------------------------
+#n_t2 = data.frame(table(lf_t2_comp$mediator))
+n_t2_all = data.frame(table(lf_t2$mediator))
+colnames(n_t2_all) = c("mediator", "n_obs")
+n_t2_all$mediator = as.character(n_t2_all$mediator)
+
+n_t2_comp = data.frame(table(lf_t2_comp$mediator))
+colnames(n_t2_comp) = c("mediator", "n_obs_complete")
+n_t2_comp$mediator = as.character(n_t2_comp$mediator)
+
+merge_cols = intersect(names(n_t2_all), names(n_t2_comp)); merge_cols
+n_t2= merge(n_t2_all, n_t2_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
+merging_cols = intersect(names(stats_un_t2), names(n_t2)); merging_cols
+if (all(n_t2$mediator%in%stats_un_t2$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t2)
+      , "\nncol:", ncol(stats_un_t2))
+}else{
+  nf = n_t2$mediator[!n_t2$mediator%in%stats_un_t2$mediator]
+  stats_un_t2 = merge(stats_un_t2, n_t2, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t2)
+      , "\nncol:", ncol(stats_un_t2))
+}
+
+# add bonferroni adjustment as well
+stats_un_t2$p_adj_bonferroni = p.adjust(stats_un_t2$p, method = "bonferroni")
+
+rm(n_t2)
+rm(lf_t2_comp)
+
+#==============
+# unpaired: t3
+#==============
+lf_t3 = lf[lf$timepoint == "t3",]
+lf_t3_comp = lf_t3[-which(is.na(lf_t3$value)),]
+
+stats_un_t3 = compare_means(value~obesity
+                            , group.by = "mediator"
+                            #, data = lf_t3 
+                            , data = lf_t3_comp
+                            , paired = FALSE
+                            , p.adjust.method = my_adjust_method)
+
+# add timepoint and convert to df
+stats_un_t3$timepoint = "t3"
+stats_un_t3 = as.data.frame(stats_un_t3)
+class(stats_un_t3)
+
+#----------------------------------------
+# calculate n_obs for each mediator: t3
+#----------------------------------------
+#n_t3 = data.frame(table(lf_t3_comp$mediator))
+n_t3_all = data.frame(table(lf_t3$mediator))
+colnames(n_t3_all) = c("mediator", "n_obs")
+n_t3_all$mediator = as.character(n_t3_all$mediator)
+
+n_t3_comp = data.frame(table(lf_t3_comp$mediator))
+colnames(n_t3_comp) = c("mediator", "n_obs_complete")
+n_t3_comp$mediator = as.character(n_t3_comp$mediator)
+
+merge_cols = intersect(names(n_t3_all), names(n_t3_comp)); merge_cols
+n_t3 = merge(n_t3_all, n_t3_comp, by = merge_cols, all = T)
+
+#==================================
+# Merge: merge stats + n_obs df
+#==================================
+merging_cols = intersect(names(stats_un_t3), names(n_t3)); merging_cols
+if (all(n_t3$mediator%in%stats_un_t3$mediator)) {
+  cat("PASS: merging stats and n_obs on column/s:", merging_cols)
+  stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all = T)
+  cat("\nsuccessfull merge:"
+      , "\nnrow:", nrow(stats_un_t3)
+      , "\nncol:", ncol(stats_un_t3))
+}else{
+  nf = n_t3$mediator[!n_t3$mediator%in%stats_un_t3$mediator]
+  stats_un_t3 = merge(stats_un_t3, n_t3, by = merging_cols, all.y = T)
+  cat("\nMerged with caution:"
+      , "\nnrows mismatch:", nf
+      , "not found in stats possibly due to all obs being LLODs"
+      , "\nintroduced NAs for:", nf
+      , "\nnrow:", nrow(stats_un_t3)
+      , "\nncol:", ncol(stats_un_t3))
+}
+
+# add bonferroni adjustment as well
+stats_un_t3$p_adj_bonferroni = p.adjust(stats_un_t3$p, method = "bonferroni")
+
+rm(n_t3)
+rm(lf_t3_comp)
+########################################################################
+#==============
+# Rbind these dfs
+#==============
+str(stats_un_t1);str(stats_un_t2); str(stats_un_t3)
+
+n_dfs = 3
+
+if ( all.equal(nrow(stats_un_t1), nrow(stats_un_t2), nrow(stats_un_t3)) && 
+     all.equal(ncol(stats_un_t1), ncol(stats_un_t2), ncol(stats_un_t3)) ) {
+  expected_rows = nrow(stats_un_t1) * n_dfs
+  expected_cols = ncol(stats_un_t1)
+  print("PASS: expected_rows and cols variables generated for downstream sanity checks")
+}else{
+  cat("FAIL: dfs have different no. of rows and cols"
+      , "\nCheck harcoded value of n_dfs"
+      , "\nexpected_rows and cols could not be generated")
+  quit()
+}
+
+if ( all.equal(colnames(stats_un_t1), colnames(stats_un_t2), colnames(stats_un_t3)) ){
+  print("PASS: colnames match. Rbind the 3 dfs...")
+  combined_unpaired_stats = rbind(stats_un_t1, stats_un_t2, stats_un_t3)
+} else{
+  cat("FAIL: cannot combined dfs. Colnames don't match!")
+  quit()
+} 
+
+if ( nrow(combined_unpaired_stats) == expected_rows && ncol(combined_unpaired_stats) == expected_cols ){
+  cat("PASS: combined_df has expected dimension"
+      , "\nNo. of rows in combined_df:", nrow(combined_unpaired_stats)
+      , "\nNo. of cols in combined_df:", ncol(combined_unpaired_stats) )
+}else{
+  cat("FAIL: combined_df dimension mismatch")
+  quit()
+}
+
+#######################################################################
+#=================
+# formatting df
+#=================
+# delete: unnecessary column
+combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.))
+
+# add sample_type
+cat("Adding sample type info as a column", my_sample_type, "...")
+combined_unpaired_stats$sample_type = my_sample_type
+
+# add: reflect stats method correctly i.e paired or unpaired
+# incase there are NA due to LLODs, the gsub won't work!
+#combined_unpaired_stats$method = gsub("Wilcoxon", "Wilcoxon_unpaired", combined_unpaired_stats$method)
+combined_unpaired_stats$method = "wilcoxon unpaired"
+combined_unpaired_stats$method
+
+# add an extra column for padjust_signif: my_adjust_method
+combined_unpaired_stats$padjust_signif = combined_unpaired_stats$p.adj
+# add appropriate symbols for padjust_signif: my_adjust_method
+combined_unpaired_stats  = dplyr::mutate(combined_unpaired_stats, padjust_signif = case_when(padjust_signif == 0.05 ~ "."
+                                                                                             , padjust_signif <=0.0001 ~ '****'
+                                                                                             , padjust_signif <=0.001 ~ '***'
+                                                                                             , padjust_signif <=0.01 ~ '**'
+                                                                                             , padjust_signif <0.05 ~ '*'
+                                                                                             , TRUE ~ 'ns'))
+# add an extra column for p_bon_signif
+combined_unpaired_stats$p_bon_signif = combined_unpaired_stats$p_adj_bonferroni
+# add appropriate symbols for p_bon_signif
+combined_unpaired_stats  = dplyr::mutate(combined_unpaired_stats, p_bon_signif = case_when(p_bon_signif == 0.05 ~ "."
+                                                                                             , p_bon_signif <=0.0001 ~ '****'
+                                                                                             , p_bon_signif <=0.001 ~ '***'
+                                                                                             , p_bon_signif <=0.01 ~ '**'
+                                                                                             , p_bon_signif <0.05 ~ '*'
+                                                                                             , TRUE ~ 'ns'))
+# reorder columns
+print("preparing to reorder columns...")
+colnames(combined_unpaired_stats)
+my_col_order2 = c("mediator"
+                  , "timepoint"
+                  , "sample_type"
+                  , "n_obs"
+                  , "n_obs_complete"
+                  , "group1"
+                  , "group2"
+                  , "method"
+                  , "p"
+                  , "p.format"
+                  , "p.signif"
+                  , "p.adj"
+                  , "padjust_signif"
+                  , "p_adj_bonferroni"
+                  , "p_bon_signif")
+
+if( length(my_col_order2) == ncol(combined_unpaired_stats) && (all(my_col_order2%in%colnames(combined_unpaired_stats))) ){
+  print("PASS: Reordering columns...")
+  combined_unpaired_stats_f = combined_unpaired_stats[, my_col_order2]
+  print("Successful: column reordering")
+  print("formatted df called:'combined_unpaired_stats_f'")
+  cat('\nformatted df has the following dimensions\n')
+  print(dim(combined_unpaired_stats_f ))
+} else{
+  cat(paste0("FAIL:Cannot reorder columns, length mismatch"
+             , "\nExpected column order for: ", ncol(combined_unpaired_stats)
+             , "\nGot:", length(my_col_order2)))
+  quit()
+}  
+# assign nice column names like replace "." with "_"
+colnames(combined_unpaired_stats_f) = c("mediator"
+                                        , "timepoint" 
+                                        , "sample_type"
+                                        , "n_obs" 
+                                        , "n_obs_complete"
+                                        , "group1"
+                                        , "group2"
+                                        , "method" 
+                                        , "p" 
+                                        , "p_format"
+                                        , "p_signif"
+                                        , paste0("p_adj_fdr_", my_adjust_method)
+                                        , paste0("p_", my_adjust_method, "_signif")
+                                        , "p_adj_bonferroni"
+                                        , "p_bon_signif")
+
+colnames(combined_unpaired_stats_f)
+
+#---------------
+# quick summary
+#---------------
+# count how many meds are significant
+n_sig  = length(combined_unpaired_stats_f$mediator[combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)])
+sig_meds = combined_unpaired_stats_f[(combined_unpaired_stats_f$p_signif<0.05 & !is.na(combined_unpaired_stats_f$p_signif<0.05)),]
+
+sig_meds$med_time = paste0(sig_meds$mediator, "@", sig_meds$timepoint)
+
+cat("\nTotal no. of statistically significant mediators in", toupper(my_sample_type)
+    , "are:", n_sig
+    , "\nThese are:", sig_meds$med_time)
+
+########################################################################
+#******************
+# write output file
+#******************
+cat("UNpaired stats for groups will be:", flu_stats_time_unpaired)
+#write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
+
+
+