#!/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_formatting_non_asthmatics.R") # check: adult variable and age variable discrepancy! metadata_all$mosaic[metadata_all$adult==1 & metadata_all$age<=18] #============= # Output #============= outfile_name = paste0("flu_stats_time_unpaired_NA_", my_sample_type, ".csv") flu_stats_time_unpaired_na = paste0(outdir_stats, outfile_name) #=============================== # data assignment for stats #================================ wf = sam_wf[sam_wf$flustat == 1,] lf = sam_lf[sam_lf$flustat == 1,] 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)) ######################################################################## # 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_na) #write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired_na, row.names = FALSE)