separated checks from flu stats and added bonferrroni p
This commit is contained in:
parent
8e06f382c0
commit
fe382dbf89
2 changed files with 62 additions and 47 deletions
|
|
@ -12,10 +12,10 @@ maindir = "~/git/mosaic_2020/"
|
||||||
outdir = paste0(maindir, "output/")
|
outdir = paste0(maindir, "output/")
|
||||||
ifelse(!dir.exists(outdir), dir.create(outdir), FALSE)
|
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)
|
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)
|
ifelse(!dir.exists(outdir_plots), dir.create(outdir_plots), FALSE)
|
||||||
########################################################################
|
########################################################################
|
||||||
# static file read: csv
|
# static file read: csv
|
||||||
|
|
@ -23,7 +23,7 @@ ifelse(!dir.exists(outdir_plots), dir.create(outdir_plots), FALSE)
|
||||||
# all patients
|
# all patients
|
||||||
#==============
|
#==============
|
||||||
all_df <- read.csv("/home/backup/MOSAIC/MEDIATOR_Data/master_file/Mosaic_master_file_from_stata.csv"
|
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 columns
|
||||||
meta_data_cols = c("mosaic", "gender", "age", "adult", "flustat", "type"
|
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")
|
#hc_data<- read.csv("/home/backup/MOSAIC/MEDIATOR_Data/master_file/Mediators_for_HC.csv")
|
||||||
#str(hc_data)
|
#str(hc_data)
|
||||||
#table(hc_data$Timepoint, hc_data$Sample)
|
#table(hc_data$Timepoint, hc_data$Sample)
|
||||||
########################################################################
|
########################################################################
|
||||||
|
|
@ -24,6 +24,7 @@ rm(pivot_cols)
|
||||||
# Output: unpaired analysis of time for npa
|
# Output: unpaired analysis of time for npa
|
||||||
#=============
|
#=============
|
||||||
stats_time_unpaired_flu_npa = paste0(outdir_stats, "flu_stats_time_unpaired_npa.csv")
|
stats_time_unpaired_flu_npa = paste0(outdir_stats, "flu_stats_time_unpaired_npa.csv")
|
||||||
|
|
||||||
#%%========================================================
|
#%%========================================================
|
||||||
# data assignment for stats
|
# data assignment for stats
|
||||||
wf = npa_df_adults_clean[npa_df_adults_clean$flustat == 1,]
|
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))
|
, "\nncol:", ncol(stats_un_t1))
|
||||||
}
|
}
|
||||||
|
|
||||||
# check: satisfied!!!!
|
# add bonferroni adjustment as well
|
||||||
wilcox.test()
|
stats_un_t1$p_adj_bonferroni = p.adjust(stats_un_t1$p, method = "bonferroni")
|
||||||
|
|
||||||
rm(n_t1)
|
rm(n_t1)
|
||||||
rm(lf_t1_comp)
|
rm(lf_t1_comp)
|
||||||
|
|
@ -134,8 +135,8 @@ if (all(n_t2$mediator%in%stats_un_t2$mediator)) {
|
||||||
, "\nncol:", ncol(stats_un_t2))
|
, "\nncol:", ncol(stats_un_t2))
|
||||||
}
|
}
|
||||||
|
|
||||||
# check: satisfied!!!!
|
# add bonferroni adjustment as well
|
||||||
wilcox.test()
|
stats_un_t2$p_adj_bonferroni = p.adjust(stats_un_t2$p, method = "bonferroni")
|
||||||
|
|
||||||
rm(n_t2)
|
rm(n_t2)
|
||||||
rm(lf_t2_comp)
|
rm(lf_t2_comp)
|
||||||
|
|
@ -182,12 +183,13 @@ if (all(n_t3$mediator%in%stats_un_t3$mediator)) {
|
||||||
, "\nncol:", ncol(stats_un_t3))
|
, "\nncol:", ncol(stats_un_t3))
|
||||||
}
|
}
|
||||||
|
|
||||||
# check: satisfied!!!!
|
# add bonferroni adjustment as well
|
||||||
wilcox.test()
|
stats_un_t3$p_adj_bonferroni = p.adjust(stats_un_t3$p, method = "bonferroni")
|
||||||
|
|
||||||
rm(n_t3)
|
rm(n_t3)
|
||||||
rm(lf_t3_comp)
|
rm(lf_t3_comp)
|
||||||
|
|
||||||
|
###########################################################################################
|
||||||
#==============
|
#==============
|
||||||
# Rbind these dfs
|
# Rbind these dfs
|
||||||
#==============
|
#==============
|
||||||
|
|
@ -224,46 +226,40 @@ if ( nrow(combined_unpaired_stats) == expected_rows && ncol(combined_unpaired_st
|
||||||
quit()
|
quit()
|
||||||
}
|
}
|
||||||
|
|
||||||
#===============================================================
|
#######################################################################
|
||||||
|
#=================
|
||||||
# formatting df
|
# formatting df
|
||||||
# delete unnecessary column
|
#=================
|
||||||
|
# delete: unnecessary column
|
||||||
combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.))
|
combined_unpaired_stats = subset(combined_unpaired_stats, select = -c(.y.))
|
||||||
|
|
||||||
# reflect stats method correctly
|
# add sample_type
|
||||||
combined_unpaired_stats$method
|
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 = gsub("Wilcoxon", "Wilcoxon_unpaired", combined_unpaired_stats$method)
|
||||||
combined_unpaired_stats$method = "wilcoxon unpaired"
|
combined_unpaired_stats$method = "wilcoxon unpaired"
|
||||||
combined_unpaired_stats$method
|
combined_unpaired_stats$method
|
||||||
|
|
||||||
# replace "." in colnames with "_"
|
# add an extra column for padjust_signif: my_adjust_method
|
||||||
colnames(combined_unpaired_stats)
|
combined_unpaired_stats$padjust_signif = combined_unpaired_stats$p.adj
|
||||||
#names(combined_unpaired_stats) = gsub("\.", "_", names(combined_unpaired_stats)) # weird!!!!
|
# add appropriate symbols for padjust_signif: my_adjust_method
|
||||||
|
|
||||||
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
|
|
||||||
combined_unpaired_stats = dplyr::mutate(combined_unpaired_stats, padjust_signif = case_when(padjust_signif == 0.05 ~ "."
|
combined_unpaired_stats = dplyr::mutate(combined_unpaired_stats, padjust_signif = case_when(padjust_signif == 0.05 ~ "."
|
||||||
, padjust_signif <=0.0001 ~ '****'
|
, padjust_signif <=0.0001 ~ '****'
|
||||||
, padjust_signif <=0.001 ~ '***'
|
, padjust_signif <=0.001 ~ '***'
|
||||||
, padjust_signif <=0.01 ~ '**'
|
, padjust_signif <=0.01 ~ '**'
|
||||||
, padjust_signif <0.05 ~ '*'
|
, padjust_signif <0.05 ~ '*'
|
||||||
, TRUE ~ 'ns'))
|
, 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
|
# reorder columns
|
||||||
print("preparing to reorder columns...")
|
print("preparing to reorder columns...")
|
||||||
colnames(combined_unpaired_stats)
|
colnames(combined_unpaired_stats)
|
||||||
|
|
@ -275,10 +271,12 @@ my_col_order2 = c("mediator"
|
||||||
, "group2"
|
, "group2"
|
||||||
, "method"
|
, "method"
|
||||||
, "p"
|
, "p"
|
||||||
, "p_format"
|
, "p.format"
|
||||||
, "p_signif"
|
, "p.signif"
|
||||||
, "p_adj"
|
, "p.adj"
|
||||||
, "padjust_signif")
|
, "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))) ){
|
if( length(my_col_order2) == ncol(combined_unpaired_stats) && (all(my_col_order2%in%colnames(combined_unpaired_stats))) ){
|
||||||
print("PASS: Reordering columns...")
|
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)))
|
, "\nGot:", length(my_col_order2)))
|
||||||
quit()
|
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
|
# write output file
|
||||||
#******************
|
#******************
|
||||||
|
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue