From bbdd2d12e56ad3f8c66981913182a6911dc1619f Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Thu, 29 Oct 2020 10:45:49 +0000 Subject: [PATCH] changed df to adults df to extract relevant info --- data_extraction_formatting.R | 109 ++++++------ read_data.R | 26 ++- stats_unpaired_npa.R | 315 ---------------------------------- stats_unpaired_sam.R | 319 ----------------------------------- stats_unpaired_serum.R | 316 ---------------------------------- 5 files changed, 78 insertions(+), 1007 deletions(-) delete mode 100755 stats_unpaired_npa.R delete mode 100755 stats_unpaired_sam.R delete mode 100755 stats_unpaired_serum.R diff --git a/data_extraction_formatting.R b/data_extraction_formatting.R index ee164c3..7b99913 100755 --- a/data_extraction_formatting.R +++ b/data_extraction_formatting.R @@ -3,7 +3,7 @@ getwd() setwd('~/git/mosaic_2020/') getwd() ######################################################################## -# TASK: Extract relevant columns from mosaic data +# TASK: Extract relevant columns from mosaic adults data # sam # serum # npa @@ -14,17 +14,17 @@ getwd() source("read_data.R") # clear unnecessary variables -#rm() +rm(all_df) ######################################################################## #========= # sam #========= sam_regex = regex(".*_sam[1-3]{1}$", ignore_case = T) -sam_cols_i = str_extract(colnames(all_df), sam_regex) # not boolean -#sam_cols_b = colnames(all_df)%in%sam_cols_i # boolean +sam_cols_i = str_extract(colnames(adult_df), sam_regex) # not boolean +#sam_cols_b = colnames(adult_df)%in%sam_cols_i # boolean -sam_cols = colnames(all_df)[colnames(all_df)%in%sam_cols_i] +sam_cols = colnames(adult_df)[colnames(adult_df)%in%sam_cols_i] # this contains log columns + daysamp_samXX: omitting these sam_regex_log_days = regex("log|day.*_sam[1-3]{1}$", ignore_case = T, perl = T) @@ -48,7 +48,7 @@ cat("Extracting SAM cols + metadata_cols") if ( length(sam_cols_to_extract) == length(meta_data_cols) + length(sam_cols_clean) ){ cat("Extracing", length(sam_cols_to_extract), "columns for sam") - sam_df = all_df[, sam_cols_to_extract] + sam_df = adult_df[, sam_cols_to_extract] }else{ cat("FAIL: length mismatch" , "Expeceted to extract:", length(meta_data_cols) + length(sam_cols_clean), "columns" @@ -61,10 +61,10 @@ colnames_sam_df = colnames(sam_df); colnames_sam_df # serum #========= serum_regex = regex(".*_serum[1-3]{1}$", ignore_case = T) -serum_cols_i = str_extract(colnames(all_df), serum_regex) # not boolean -#serum_cols_b = colnames(all_df)%in%serum_cols_i # boolean +serum_cols_i = str_extract(colnames(adult_df), serum_regex) # not boolean +#serum_cols_b = colnames(adult_df)%in%serum_cols_i # boolean -serum_cols = colnames(all_df)[colnames(all_df)%in%serum_cols_i] +serum_cols = colnames(adult_df)[colnames(adult_df)%in%serum_cols_i] # this contains log columns + dayserump_serumXX: omitting these serum_regex_log_days = regex("log|day.*_serum[1-3]{1}$", ignore_case = T, perl = T) @@ -88,7 +88,7 @@ cat("Extracting SERUM cols + metadata_cols") if ( length(serum_cols_to_extract) == length(meta_data_cols) + length(serum_cols_clean) ){ cat("Extracing", length(serum_cols_to_extract), "columns for serum") - serum_df = all_df[, serum_cols_to_extract] + serum_df = adult_df[, serum_cols_to_extract] }else{ cat("FAIL: length mismatch" , "Expeceted to extract:", length(meta_data_cols) + length(serum_cols_clean), "columns" @@ -101,10 +101,10 @@ colnames_serum_df = colnames(serum_df); colnames_serum_df # npa #========= npa_regex = regex(".*_npa[1-3]{1}$", ignore_case = T) -npa_cols_i = str_extract(colnames(all_df), npa_regex) # not boolean -#npa_cols_b = colnames(all_df)%in%npa_cols_i # boolean +npa_cols_i = str_extract(colnames(adult_df), npa_regex) # not boolean +#npa_cols_b = colnames(adult_df)%in%npa_cols_i # boolean -npa_cols = colnames(all_df)[colnames(all_df)%in%npa_cols_i] +npa_cols = colnames(adult_df)[colnames(adult_df)%in%npa_cols_i] # this contains log columns + daynpap_npaXX: omitting these npa_regex_log_days = regex("log|day|vl_samptime|ct.*_npa[1-3]{1}$", ignore_case = T, perl = T) @@ -128,7 +128,7 @@ cat("Extracting NPA cols + metadata_cols") if ( length(npa_cols_to_extract) == length(meta_data_cols) + length(npa_cols_clean) ){ cat("Extracing", length(npa_cols_to_extract), "columns for npa") - npa_df = all_df[, npa_cols_to_extract] + npa_df = adult_df[, npa_cols_to_extract] }else{ cat("FAIL: length mismatch" , "Expeceted to extract:", length(meta_data_cols) + length(npa_cols_clean), "columns" @@ -137,8 +137,11 @@ if ( length(npa_cols_to_extract) == length(meta_data_cols) + length(npa_cols_cle colnames_npa_df = colnames(npa_df); colnames_npa_df +#============== +# quick checks +#============== colnames_check = as.data.frame(cbind(colnames_sam_df, colnames_serum_df, colnames_npa_df)) -tail(colnames_check) +tail(colnames_check) # gives a warning message due to differeing no. of rows for cbind! # put NA where a match doesn't exist # unmatched lengths @@ -168,12 +171,16 @@ quick_check = as.data.frame(cbind(metadata_all$mosaic , metadata_all$adult , metadata_all$age , metadata_all$obesity - , metadata_all$obese2)) + , metadata_all$obese2 + )) colnames(quick_check) = c("mosaic", "adult", "age", "obesity", "obese2") ########################################################################## # LF data ########################################################################## +cols_to_omit = c("adult", "obese2" + , "height", "height_unit", "weight" + , "weight_unit", "visual_est_bmi", "bmi_rating") #============== # lf data: sam @@ -181,19 +188,11 @@ colnames(quick_check) = c("mosaic", "adult", "age", "obesity", "obese2") str(sam_df) table(sam_df$obesity); table(sam_df$obese2) -sam_df_adults = sam_df[sam_df$adult == 1,] - -cols_to_omit = c("type" - #, "flustat" - #, "obesity" - #, "obese2" - , "height", "height_unit", "weight" - , "weight_unit", "visual_est_bmi", "bmi_rating") - -#sam_df_adults_clean = sam_df_adults[!cols_to_omit] +#sam_df_adults = sam_df[sam_df$adult == 1,] # resolved at source and only dealing wit age as adult +sam_df_adults = sam_df wf_cols = colnames(sam_df_adults)[!colnames(sam_df_adults)%in%cols_to_omit] -sam_df_adults_clean = sam_df_adults[wf_cols] +sam_wf = sam_df_adults[wf_cols] pivot_cols = meta_data_cols # subselect pivot_cols @@ -208,25 +207,25 @@ if (length(pivot_cols) == length(meta_data_cols) - length(cols_to_omit)){ quit() } -expected_rows_sam_lf = nrow(sam_df_adults_clean) * (length(sam_df_adults_clean) - length(pivot_cols)); expected_rows_sam_lf +expected_rows_sam_lf = nrow(sam_wf) * (length(sam_wf) - length(pivot_cols)); expected_rows_sam_lf # using regex: -sam_adults_lf = sam_df_adults_clean %>% +sam_lf = sam_wf %>% tidyr::pivot_longer(-all_of(pivot_cols) , names_to = c("mediator", "sample_type", "timepoint") , names_pattern = "(.*)_(.*)([1-3]{1})" , values_to = "value") if ( - (nrow(sam_adults_lf) == expected_rows_sam_lf) & (sum(table(is.na(sam_adults_lf$mediator))) == expected_rows_sam_lf) + (nrow(sam_lf) == expected_rows_sam_lf) & (sum(table(is.na(sam_lf$mediator))) == expected_rows_sam_lf) ) { cat(paste0("PASS: long format data has correct no. of rows and NA in mediator:" - , "\nNo. of rows: ", nrow(sam_adults_lf) - , "\nNo. of cols: ", ncol(sam_adults_lf))) + , "\nNo. of rows: ", nrow(sam_lf) + , "\nNo. of cols: ", ncol(sam_lf))) } else{ cat(paste0("FAIL:long format data has unexpected no. of rows or NAs in mediator" , "\nExpected no. of rows: ", expected_rows_sam_lf - , "\nGot: ", nrow(sam_adults_lf) + , "\nGot: ", nrow(sam_lf) , "\ncheck expected rows calculation!")) quit() } @@ -241,11 +240,11 @@ if ( str(serum_df) table(serum_df$obesity); table(serum_df$obese2) -serum_df_adults = serum_df[serum_df$adult == 1,] +#serum_df_adults = serum_df[serum_df$adult == 1,] # extract based on age +serum_df_adults = serum_df -#serum_df_adults_clean = serum_df_adults[!cols_to_omit] wf_cols = colnames(serum_df_adults)[!colnames(serum_df_adults)%in%cols_to_omit] -serum_df_adults_clean = serum_df_adults[wf_cols] +serum_wf = serum_df_adults[wf_cols] pivot_cols = meta_data_cols pivot_cols = meta_data_cols[!meta_data_cols%in%cols_to_omit];pivot_cols @@ -259,25 +258,25 @@ if (length(pivot_cols) == length(meta_data_cols) - length(cols_to_omit)){ quit() } -expected_rows_serum_lf = nrow(serum_df_adults_clean) * (length(serum_df_adults_clean) - length(pivot_cols)); expected_rows_serum_lf +expected_rows_serum_lf = nrow(serum_wf) * (length(serum_wf) - length(pivot_cols)); expected_rows_serum_lf # using regex: -serum_adults_lf = serum_df_adults_clean %>% +serum_lf = serum_wf %>% tidyr::pivot_longer(-all_of(pivot_cols) , names_to = c("mediator", "sample_type", "timepoint") , names_pattern = "(.*)_(.*)([1-3]{1})" , values_to = "value") if ( - (nrow(serum_adults_lf) == expected_rows_serum_lf) & (sum(table(is.na(serum_adults_lf$mediator))) == expected_rows_serum_lf) + (nrow(serum_lf) == expected_rows_serum_lf) & (sum(table(is.na(serum_lf$mediator))) == expected_rows_serum_lf) ) { cat(paste0("PASS: long format data has correct no. of rows and NA in mediator:" - , "\nNo. of rows: ", nrow(serum_adults_lf) - , "\nNo. of cols: ", ncol(serum_adults_lf))) + , "\nNo. of rows: ", nrow(serum_lf) + , "\nNo. of cols: ", ncol(serum_lf))) } else{ cat(paste0("FAIL:long format data has unexpected no. of rows or NAs in mediator" , "\nExpected no. of rows: ", expected_rows_serum_lf - , "\nGot: ", nrow(serum_adults_lf) + , "\nGot: ", nrow(serum_lf) , "\ncheck expected rows calculation!")) quit() } @@ -288,11 +287,11 @@ if ( str(npa_df) table(npa_df$obesity); table(npa_df$obese2) -npa_df_adults = npa_df[npa_df$adult == 1,] -#npa_df_adults_clean = npa_df_adults[!cols_to_omit] +#npa_df_adults = npa_df[npa_df$adult == 1,] # extract based on age +npa_df_adults = npa_df wf_cols = colnames(npa_df_adults)[!colnames(npa_df_adults)%in%cols_to_omit] -npa_df_adults_clean = npa_df_adults[wf_cols] +npa_wf = npa_df_adults[wf_cols] pivot_cols = meta_data_cols pivot_cols = meta_data_cols[!meta_data_cols%in%cols_to_omit];pivot_cols @@ -306,25 +305,25 @@ if (length(pivot_cols) == length(meta_data_cols) - length(cols_to_omit)){ quit() } -expected_rows_npa_lf = nrow(npa_df_adults_clean) * (length(npa_df_adults_clean) - length(pivot_cols)); expected_rows_npa_lf +expected_rows_npa_lf = nrow(npa_wf) * (length(npa_wf) - length(pivot_cols)); expected_rows_npa_lf # using regex: -npa_adults_lf = npa_df_adults_clean %>% +npa_lf = npa_wf %>% tidyr::pivot_longer(-all_of(pivot_cols) , names_to = c("mediator", "sample_type", "timepoint") , names_pattern = "(.*)_(.*)([1-3]{1})" , values_to = "value") if ( - (nrow(npa_adults_lf) == expected_rows_npa_lf) & (sum(table(is.na(npa_adults_lf$mediator))) == expected_rows_npa_lf) + (nrow(npa_lf) == expected_rows_npa_lf) & (sum(table(is.na(npa_lf$mediator))) == expected_rows_npa_lf) ) { cat(paste0("PASS: long format data has correct no. of rows and NA in mediator:" - , "\nNo. of rows: ", nrow(npa_adults_lf) - , "\nNo. of cols: ", ncol(npa_adults_lf))) + , "\nNo. of rows: ", nrow(npa_lf) + , "\nNo. of cols: ", ncol(npa_lf))) } else{ cat(paste0("FAIL:long format data has unexpected no. of rows or NAs in mediator" , "\nExpected no. of rows: ", expected_rows_npa_lf - , "\nGot: ", nrow(npa_adults_lf) + , "\nGot: ", nrow(npa_lf) , "\ncheck expected rows calculation!")) quit() } @@ -334,12 +333,14 @@ if ( rm(sam_regex, sam_regex_log_days, sam_cols, sam_cols_clean, sam_cols_i, sam_cols_to_extract, sam_cols_to_omit) rm(serum_regex, serum_regex_log_days, serum_cols, serum_cols_clean, serum_cols_i, serum_cols_to_extract, serum_cols_to_omit) rm(npa_regex, npa_regex_log_days, npa_cols, npa_cols_clean, npa_cols_i, npa_cols_to_extract, npa_cols_to_omit) -rm(all_df) +rm(adult_df) rm(colnames_check) -rm(i, j, expected_cols, start, wf_cols, extra_cols, cols_to_omit) +rm(i, j + #, expected_cols + , start, wf_cols, extra_cols, cols_to_omit) # rm not_clean dfs rm(sam_df_adults, serum_df_adults, npa_df_adults) -# rm df containing non-adults +# rm df rm(sam_df, serum_df, npa_df) diff --git a/read_data.R b/read_data.R index 82ca50f..9237058 100755 --- a/read_data.R +++ b/read_data.R @@ -26,9 +26,16 @@ all_df <- read.csv("/home/backup/MOSAIC/MEDIATOR_Data/master_file/Mosaic_master_ , fileEncoding = 'latin1') # meta data columns -meta_data_cols = c("mosaic", "gender", "age", "adult", "flustat", "type" - , "obesity", "obese2", "height", "height_unit", "weight" - , "weight_unit", "visual_est_bmi", "bmi_rating") +meta_data_cols = c("mosaic", "gender", "age" + , "adult" + , "flustat", "type" + , "obesity" + , "obese2" + , "height", "height_unit" + , "weight", "weight_unit" + , "ia_height_ftin", "ia_height_m", "ia_weight" + , "visual_est_bmi", "bmi_rating" + ) # check if these columns to select are present in the data meta_data_cols%in%colnames(all_df) @@ -36,6 +43,19 @@ all(meta_data_cols%in%colnames(all_df)) metadata_all = all_df[, meta_data_cols] + +#============== +# adult patients +#============== + +adult_df = all_df[all_df$age>=18,] + +if (table(adult_df$adult == 1)[[1]] == nrow(adult_df) ){ + cat ("PASS: adult df extracted successfully") +} else{ + cat ("FAIL: adult df number mismatch!") +} + #============ # hc #============ diff --git a/stats_unpaired_npa.R b/stats_unpaired_npa.R deleted file mode 100755 index 560ffb4..0000000 --- a/stats_unpaired_npa.R +++ /dev/null @@ -1,315 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd("~/git/mosaic_2020/") -getwd() -############################################################ -# TASK: unpaired (time) analysis of mediators: NPA -############################################################ -#============= -# Input -#============= -source("data_extraction_formatting.R") - -# clear variables -rm(sam_adults_lf, sam_df_adults_clean -, serum_adults_lf, serum_df_adults_clean) -rm(colnames_sam_df, expected_rows_sam_lf -, colnames_serum_df, expected_rows_serum_lf) - -rm(pivot_cols) - -my_sample_type = "npa" -#============= -# Output: unpaired analysis of time for npa -#============= -outfile_name = paste0("stats_time_unpaired_", my_sample_type, ".csv") -stats_time_unpaired = paste0(outdir_stats, outfile_name) -#%%======================================================== -# data assignment for stats -wf = npa_df_adults_clean -lf = npa_adults_lf -#%%======================================================== -table(lf$timepoint) -lf$timepoint = paste0("t", lf$timepoint) - -######################################################################## -# 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)] - -stats_un_t1$timepoint = "t1" - -stats_un_t1 = as.data.frame(stats_un_t1) -class(stats_un_t1) - -# calculate n_obs for each mediator -n_t1 = data.frame(table(lf_t1_comp$mediator)) -colnames(n_t1) = c("mediator", "n_obs") -n_t1$mediator = as.character(n_t1$mediator) - -# 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) -stats_un_t2$timepoint = "t2" - -stats_un_t2 = as.data.frame(stats_un_t2) -class(stats_un_t2) - -# calculate n_obs for each mediator -n_t2 = data.frame(table(lf_t2_comp$mediator)) -colnames(n_t2) = c("mediator", "n_obs") -n_t2$mediator = as.character(n_t2$mediator) - -# 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) - -stats_un_t3$timepoint = "t3" - -stats_un_t3 = as.data.frame(stats_un_t3) -class(stats_un_t3) - -# calculate n_obs for each mediator -n_t3 = data.frame(table(lf_t3_comp$mediator)) -colnames(n_t3) = c("mediator", "n_obs") -n_t3$mediator = as.character(n_t3$mediator) - -# 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" - , "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" - , "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 -#****************** -cat("UNpaired stats for groups will be:", stats_time_unpaired) -write.csv(combined_unpaired_stats_f, stats_time_unpaired, row.names = FALSE) diff --git a/stats_unpaired_sam.R b/stats_unpaired_sam.R deleted file mode 100755 index 0c624ac..0000000 --- a/stats_unpaired_sam.R +++ /dev/null @@ -1,319 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd("~/git/mosaic_2020/") -getwd() -############################################################ -# TASK: unpaired (time) analysis of mediators: SAM -############################################################ -#============= -# Input -#============= -source("data_extraction_formatting.R") - -# clear variables -rm(npa_adults_lf, npa_df_adults_clean -, serum_adults_lf, serum_df_adults_clean) -rm(colnames_npa_df, expected_rows_npa_lf -, colnames_serum_df, expected_rows_serum_lf) - -rm(pivot_cols) - -my_sample_type = "sam" -#============= -# Output: unpaired analysis of time for sam -#============= -outfile_name = paste0("stats_time_unpaired_", my_sample_type, ".csv") -stats_time_unpaired = paste0(outdir_stats, outfile_name) -#%%======================================================== -# data assignment for stats -wf = sam_df_adults_clean -lf = sam_adults_lf -#%%======================================================== -table(lf$timepoint) -lf$timepoint = paste0("t", lf$timepoint) - -######################################################################## -# 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)] - -stats_un_t1$timepoint = "t1" - -stats_un_t1 = as.data.frame(stats_un_t1) -class(stats_un_t1) - -# calculate n_obs for each mediator -n_t1 = data.frame(table(lf_t1_comp$mediator)) -colnames(n_t1) = c("mediator", "n_obs") -n_t1$mediator = as.character(n_t1$mediator) - -# 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) -stats_un_t2$timepoint = "t2" - -stats_un_t2 = as.data.frame(stats_un_t2) -class(stats_un_t2) - -# calculate n_obs for each mediator -n_t2 = data.frame(table(lf_t2_comp$mediator)) -colnames(n_t2) = c("mediator", "n_obs") -n_t2$mediator = as.character(n_t2$mediator) - -# 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) - -stats_un_t3$timepoint = "t3" - -stats_un_t3 = as.data.frame(stats_un_t3) -class(stats_un_t3) - -# calculate n_obs for each mediator -n_t3 = data.frame(table(lf_t3_comp$mediator)) -colnames(n_t3) = c("mediator", "n_obs") -n_t3$mediator = as.character(n_t3$mediator) - -# 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" - , "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" - , "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 -#****************** -cat("UNpaired stats for groups will be:", stats_time_unpaired) -write.csv(combined_unpaired_stats_f, stats_time_unpaired, row.names = FALSE) diff --git a/stats_unpaired_serum.R b/stats_unpaired_serum.R deleted file mode 100755 index 2d4b096..0000000 --- a/stats_unpaired_serum.R +++ /dev/null @@ -1,316 +0,0 @@ -#!/usr/bin/Rscript -getwd() -setwd("~/git/mosaic_2020/") -getwd() -############################################################ -# TASK: unpaired (time) analysis of mediators: serum -############################################################ -#============= -# Input -#============= -source("data_extraction_formatting.R") - -# clear variables -rm(sam_adults_lf, sam_df_adults_clean -, npa_adults_lf, npa_df_adults_clean) -rm(colnames_sam_df, expected_rows_sam_lf -, colnames_npa_df, expected_rows_npa_lf) - -rm(pivot_cols) - -my_sample_type = "serum" -#============= -# Output: unpaired analysis of time for serum -#============= -outfile_name = paste0("stats_time_unpaired_", my_sample_type, ".csv") -stats_time_unpaired = paste0(outdir_stats, outfile_name) -#%%======================================================== -# data assignment for stats -wf = serum_df_adults_clean -lf = serum_adults_lf -#%%======================================================== -table(lf$timepoint) -lf$timepoint = paste0("t", lf$timepoint) - -######################################################################## -# 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)] - -stats_un_t1$timepoint = "t1" - -stats_un_t1 = as.data.frame(stats_un_t1) -class(stats_un_t1) - -# calculate n_obs for each mediator -n_t1 = data.frame(table(lf_t1_comp$mediator)) -colnames(n_t1) = c("mediator", "n_obs") -n_t1$mediator = as.character(n_t1$mediator) - -# 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) -stats_un_t2$timepoint = "t2" - -stats_un_t2 = as.data.frame(stats_un_t2) -class(stats_un_t2) - -# calculate n_obs for each mediator -n_t2 = data.frame(table(lf_t2_comp$mediator)) -colnames(n_t2) = c("mediator", "n_obs") -n_t2$mediator = as.character(n_t2$mediator) - -# 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) - -stats_un_t3$timepoint = "t3" - -stats_un_t3 = as.data.frame(stats_un_t3) -class(stats_un_t3) - - -# calculate n_obs for each mediator -n_t3 = data.frame(table(lf_t3_comp$mediator)) -colnames(n_t3) = c("mediator", "n_obs") -n_t3$mediator = as.character(n_t3$mediator) - -# 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" - , "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" - , "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 -#****************** -cat("UNpaired stats for groups will be:", stats_time_unpaired) -write.csv(combined_unpaired_stats_f, stats_time_unpaired, row.names = FALSE)