300 lines
No EOL
10 KiB
R
300 lines
No EOL
10 KiB
R
#!/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)
|
|
|
|
#=============
|
|
# Output: unpaired analysis of time for sam
|
|
#=============
|
|
stats_time_unpaired_sam = paste0(outdir_stats, "stats_time_unpaired_sam.csv")
|
|
#%%========================================================
|
|
# 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))
|
|
}
|
|
|
|
# check: satisfied!!!!
|
|
wilcox.test()
|
|
|
|
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))
|
|
}
|
|
|
|
# check: satisfied!!!!
|
|
wilcox.test()
|
|
|
|
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])
|
|
|
|
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.))
|
|
|
|
# reflect stats method correctly
|
|
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
|
|
|
|
# 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 = "sam"
|
|
|
|
# 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 ~ "."
|
|
, padjust_signif <=0.0001 ~ '****'
|
|
, padjust_signif <=0.001 ~ '***'
|
|
, padjust_signif <=0.01 ~ '**'
|
|
, padjust_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")
|
|
|
|
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()
|
|
}
|
|
|
|
#******************
|
|
# write output file
|
|
#******************
|
|
cat("UNpaired stats for groups will be:", stats_time_unpaired_sam)
|
|
write.csv(combined_unpaired_stats_f, stats_time_unpaired_sam, row.names = FALSE) |