diff --git a/read_data.R b/read_data.R index 97c8c78..82ca50f 100644 --- a/read_data.R +++ b/read_data.R @@ -12,10 +12,10 @@ maindir = "~/git/mosaic_2020/" outdir = paste0(maindir, "output/") ifelse(!dir.exists(outdir), dir.create(outdir), FALSE) -outdir_stats = "~/git/mosaic_2020/output/stats/" +outdir_stats = paste0(maindir, "output/stats/") ifelse(!dir.exists(outdir_stats), dir.create(outdir_stats), FALSE) -outdir_plots = "~/git/mosaic_2020/output/plots/" +outdir_plots = paste0(maindir, "output/plots") ifelse(!dir.exists(outdir_plots), dir.create(outdir_plots), FALSE) ######################################################################## # static file read: csv @@ -23,7 +23,7 @@ ifelse(!dir.exists(outdir_plots), dir.create(outdir_plots), FALSE) # all patients #============== all_df <- read.csv("/home/backup/MOSAIC/MEDIATOR_Data/master_file/Mosaic_master_file_from_stata.csv" - , fileEncoding='latin1') + , fileEncoding = 'latin1') # meta data columns meta_data_cols = c("mosaic", "gender", "age", "adult", "flustat", "type" @@ -42,4 +42,4 @@ metadata_all = all_df[, meta_data_cols] #hc_data<- read.csv("/home/backup/MOSAIC/MEDIATOR_Data/master_file/Mediators_for_HC.csv") #str(hc_data) #table(hc_data$Timepoint, hc_data$Sample) -######################################################################## +######################################################################## \ No newline at end of file diff --git a/stats_unpaired_flu_npa.R b/stats_unpaired_flu_npa.R index b771f59..4810d92 100644 --- a/stats_unpaired_flu_npa.R +++ b/stats_unpaired_flu_npa.R @@ -24,6 +24,7 @@ rm(pivot_cols) # Output: unpaired analysis of time for npa #============= stats_time_unpaired_flu_npa = paste0(outdir_stats, "flu_stats_time_unpaired_npa.csv") + #%%======================================================== # data assignment for stats wf = npa_df_adults_clean[npa_df_adults_clean$flustat == 1,] @@ -87,8 +88,8 @@ if (all(n_t1$mediator%in%stats_un_t1$mediator)) { , "\nncol:", ncol(stats_un_t1)) } -# check: satisfied!!!! -wilcox.test() +# 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) @@ -134,8 +135,8 @@ if (all(n_t2$mediator%in%stats_un_t2$mediator)) { , "\nncol:", ncol(stats_un_t2)) } -# check: satisfied!!!! -wilcox.test() +# 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) @@ -182,12 +183,13 @@ if (all(n_t3$mediator%in%stats_un_t3$mediator)) { , "\nncol:", ncol(stats_un_t3)) } -# check: satisfied!!!! -wilcox.test() +# 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 #============== @@ -224,46 +226,40 @@ if ( nrow(combined_unpaired_stats) == expected_rows && ncol(combined_unpaired_st quit() } -#=============================================================== +####################################################################### +#================= # formatting df -# delete unnecessary column +#================= +# delete: unnecessary column combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.)) -# reflect stats method correctly -combined_unpaired_stats$method +# add sample_type +combined_unpaired_stats$sample_type = "npa" + +# 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 -# replace "." in colnames with "_" -colnames(combined_unpaired_stats) -#names(combined_unpaired_stats) = gsub("\.", "_", names(combined_unpaired_stats)) # weird!!!! - -colnames(combined_unpaired_stats) = c("mediator" - , "group1" - , "group2" - , "p" - , "p_adj" - , "p_format" - , "p_signif" - , "method" - , "timepoint" - , "n_obs") - -colnames(combined_unpaired_stats) -combined_unpaired_stats$sample_type = "npa" - -# add an extra column for padjust_signif -combined_unpaired_stats$padjust_signif = round(combined_unpaired_stats$p_adj, digits = 2) - -# add appropriate symbols for padjust_signif +# 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')) - + , 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) @@ -275,10 +271,12 @@ my_col_order2 = c("mediator" , "group2" , "method" , "p" - , "p_format" - , "p_signif" - , "p_adj" - , "padjust_signif") + , "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...") @@ -293,7 +291,24 @@ if( length(my_col_order2) == ncol(combined_unpaired_stats) && (all(my_col_order2 , "\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" + , "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) + #****************** # write output file #******************