added unpaired stats for asthma and severe for sam, serum and npa
This commit is contained in:
parent
9acc36c8de
commit
28a356fc19
8 changed files with 2741 additions and 0 deletions
183
asthma/check_stats.R
Executable file
183
asthma/check_stats.R
Executable file
|
|
@ -0,0 +1,183 @@
|
|||
#!/usr/bin/Rscript
|
||||
getwd()
|
||||
setwd("~/git/mosaic_2020/")
|
||||
getwd()
|
||||
########################################################################
|
||||
# TASK: check stats
|
||||
########################################################################
|
||||
# load libraries, packages and local imports
|
||||
source("Header_TT.R")
|
||||
|
||||
#=============
|
||||
# Input
|
||||
#=============
|
||||
source("data_extraction_mediators.R")
|
||||
#source("read_data.R")
|
||||
|
||||
########################################################################
|
||||
#=============
|
||||
# Data extraction for non-severe
|
||||
#=============
|
||||
table(fp_adults_ics$asthma)
|
||||
table(fp_adults_ics$obesity)
|
||||
|
||||
df = fp_adults_ics[fp_adults_ics$asthma == 1,]
|
||||
table(df$obesity)
|
||||
|
||||
df_ob = df[df$obesity == 1,]
|
||||
df_nob = df[df$obesity == 0,]
|
||||
|
||||
rm(npa_lf, npa_wf, sam_lf, sam_wf, serum_lf, serum_wf, fp_adults_ics)
|
||||
|
||||
#============
|
||||
# npa
|
||||
#============
|
||||
# il8: t2
|
||||
wilcox.test(df_ob$il8_serum2, df_nob$il8_serum2)
|
||||
wilcox_test(il8_serum2 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il8_serum2)
|
||||
na_c = sum(is.na(foo$df.il8_serum2))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il8_serum2, df$obesity, median, na.rm = T)
|
||||
tapply(df$il8_serum2, df$obesity, summary, na.rm = T)
|
||||
|
||||
# il8_2: t2
|
||||
wilcox.test(df_ob$il8_2_serum2, df_nob$il8_2_serum2)
|
||||
wilcox_test(il8_2_serum2 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il8_2_serum2)
|
||||
na_c = sum(is.na(foo$df.il8_2_serum2))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il8_2_serum2, df$obesity, median, na.rm = T)
|
||||
tapply(df$il8_2_serum2, df$obesity, summary, na.rm = T)
|
||||
|
||||
|
||||
# il4: t3
|
||||
wilcox.test(df_ob$il4_serum3, df_nob$il4_serum3)
|
||||
wilcox_test(il4_serum3 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il4_serum3)
|
||||
na_c = sum(is.na(foo$df.il4_serum3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il4_serum3, df$obesity, median, na.rm = T)
|
||||
tapply(df$il4_serum3, df$obesity, summary, na.rm = T)
|
||||
|
||||
# ifnb: t3
|
||||
wilcox.test(df_ob$ifnb_serum3, df_nob$ifnb_serum3)
|
||||
wilcox_test(ifnb_serum3 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$ifnb_serum3)
|
||||
na_c = sum(is.na(foo$df.ifnb_serum3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$ifnb_serum3, df$obesity, median, na.rm = T)
|
||||
tapply(df$ifnb_serum3, df$obesity, summary, na.rm = T)
|
||||
|
||||
#============
|
||||
# sam
|
||||
#============
|
||||
# tarc: t3
|
||||
wilcox.test(df_ob$tarc_sam3
|
||||
, df_nob$tarc_sam3)
|
||||
|
||||
foo = data.frame(df$obesity, df$tarc_sam3)
|
||||
na_c = sum(is.na(foo$df.tarc_sam3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$tarc_sam3, df$obesity, median, na.rm = T)
|
||||
|
||||
|
||||
# il13: t3
|
||||
wilcox.test(df_ob$il13_sam3
|
||||
, df_nob$il13_sam3)
|
||||
|
||||
foo = data.frame(df$obesity, df$il13_sam3)
|
||||
na_c = sum(is.na(foo$df.il13_sam3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il13_sam3, df$obesity, median, na.rm = T)
|
||||
|
||||
|
||||
# eotaxin: t1
|
||||
wilcox.test(df_ob$eotaxin_sam1
|
||||
, df_nob$eotaxin_sam1)
|
||||
|
||||
foo = data.frame(df$obesity, df$eotaxin_sam1)
|
||||
na_c = sum(is.na(foo$df.eotaxin_sam1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$eotaxin_sam1, df$obesity, median, na.rm = T)
|
||||
|
||||
|
||||
#============
|
||||
# serum
|
||||
#============
|
||||
# ifna2a: t1
|
||||
wilcox.test(df_ob$ifna2a_serum1,
|
||||
df_nob$ifna2a_serum1)
|
||||
|
||||
foo = data.frame(df$obesity, df$ifna2a_serum1)
|
||||
na_c = sum(is.na(foo$df.ifna2a_serum1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$ifna2a_serum1, df$obesity, median, na.rm = T)
|
||||
|
||||
|
||||
# mip17: t1
|
||||
wilcox.test(df_ob$mip17_serum1,
|
||||
df_nob$mip17_serum1)
|
||||
|
||||
foo = data.frame(df$obesity, df$mip17_serum1)
|
||||
na_c = sum(is.na(foo$df.mip17_serum1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$mip17_serum1, df$obesity, median, na.rm = T)
|
||||
|
||||
# il1: t1
|
||||
wilcox.test(df_ob$il1_serum1
|
||||
, df_nob$il1_serum1)
|
||||
|
||||
foo = data.frame(df$obesity, df$il1_serum1)
|
||||
na_c = sum(is.na(foo$df.il1_serum1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il1_serum1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il1_serum1, df$obesity, summary, na.rm = T)
|
||||
boxplot(df_ob$il1_serum1, df_nob$il1_serum1)
|
||||
|
||||
|
||||
# ifna2a: t2
|
||||
wilcox.test(df_ob$ifna2a_serum2,
|
||||
df_nob$ifna2a_serum2)
|
||||
|
||||
foo = data.frame(df$obesity, df$ifna2a_serum2)
|
||||
na_c = sum(is.na(foo$df.ifna2a_serum2))
|
||||
nrow(df) - na_c
|
||||
|
||||
|
||||
tapply(df$ifna2a_serum2, df$obesity, median, na.rm = T)
|
||||
|
||||
# ifna2a: t3
|
||||
wilcox.test(df_ob$ifna2a_serum3,
|
||||
df_nob$ifna2a_serum3)
|
||||
|
||||
foo = data.frame(df$obesity, df$ifna2a_serum3)
|
||||
na_c = sum(is.na(foo$df.ifna2a_serum3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$ifna2a_serum3, df$obesity, median, na.rm = T)
|
||||
|
||||
|
||||
# ifnb: t3
|
||||
wilcox.test(df_ob$ifnb_serum3
|
||||
, df_nob$ifnb_serum3)
|
||||
|
||||
foo = data.frame(df$obesity, df$ifnb_serum3)
|
||||
na_c = sum(is.na(foo$df.ifnb_serum3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$ifnb_serum3, df$obesity, median, na.rm = T)
|
||||
386
asthma/flu_a_stats_unpaired_npa.R
Executable file
386
asthma/flu_a_stats_unpaired_npa.R
Executable file
|
|
@ -0,0 +1,386 @@
|
|||
#!/usr/bin/Rscript
|
||||
getwd()
|
||||
setwd("~/git/mosaic_2020/")
|
||||
getwd()
|
||||
############################################################
|
||||
# TASK: unpaired (time) analysis of mediators:
|
||||
# sample type: NPA
|
||||
# data: Flu positive adult patients
|
||||
# group: obesity
|
||||
############################################################
|
||||
my_sample_type = "npa"
|
||||
|
||||
#=============
|
||||
# Input
|
||||
#=============
|
||||
#source("data_extraction_mediators.R")
|
||||
source("plot_data_a.R")
|
||||
|
||||
# check: copd and asthma conflict
|
||||
table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
|
||||
|
||||
#=============
|
||||
# Output
|
||||
#=============
|
||||
outfile_name = paste0("flu_stats_time_unpaired_a_", my_sample_type, ".csv")
|
||||
flu_stats_time_unpaired = paste0(outdir_stats_a, outfile_name)
|
||||
flu_stats_time_unpaired
|
||||
|
||||
# quick checks
|
||||
table(wf_fp_npa$T1_resp_score)
|
||||
table(lf_fp_npa$mediator, lf_fp_npa$timepoint, lf_fp_npa$T1_resp_score)
|
||||
table(wf_fp_npa$obesity)
|
||||
|
||||
#===============================
|
||||
# data assignment for stats
|
||||
#================================
|
||||
wf = wf_fp_npa
|
||||
lf = lf_fp_npa
|
||||
|
||||
identical(wf, wf_fp_npa)
|
||||
identical(lf, lf_fp_npa)
|
||||
########################################################################
|
||||
# clear variables
|
||||
rm(lf_fp_sam, lf_fp_serum
|
||||
, wf_fp_sam, wf_fp_serum)
|
||||
|
||||
# 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))
|
||||
}
|
||||
|
||||
# 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)
|
||||
write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
|
||||
395
asthma/flu_a_stats_unpaired_sam.R
Executable file
395
asthma/flu_a_stats_unpaired_sam.R
Executable file
|
|
@ -0,0 +1,395 @@
|
|||
#!/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_mediators.R")
|
||||
source("plot_data_a.R")
|
||||
|
||||
# check: copd and asthma conflict
|
||||
table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
|
||||
|
||||
#=============
|
||||
# Output
|
||||
#=============
|
||||
outfile_name = paste0("flu_stats_time_unpaired_a_", my_sample_type, ".csv")
|
||||
flu_stats_time_unpaired = paste0(outdir_stats_a, outfile_name)
|
||||
flu_stats_time_unpaired
|
||||
|
||||
# quick checks
|
||||
table(wf_fp_sam$T1_resp_score)
|
||||
table(lf_fp_sam$mediator, lf_fp_sam$timepoint, lf_fp_sam$T1_resp_score)
|
||||
table(wf_fp_sam$obesity)
|
||||
|
||||
#===============================
|
||||
# data assignment for stats
|
||||
#================================
|
||||
wf = wf_fp_sam
|
||||
lf = lf_fp_sam
|
||||
|
||||
identical(wf, wf_fp_sam)
|
||||
identical(lf, lf_fp_sam)
|
||||
|
||||
########################################################################
|
||||
# clear variables
|
||||
rm(lf_fp_npa, lf_fp_serum
|
||||
, wf_fp_npa, wf_fp_serum)
|
||||
|
||||
# sanity checks
|
||||
table(lf$timepoint)
|
||||
length(unique(lf$mosaic))
|
||||
|
||||
#if (table(lf$flustat) == table(sam_lf$flustat)[[2]]){
|
||||
# cat("Analysing Flu positive patients for:", my_sample_type)
|
||||
#}else{
|
||||
# cat("FAIL: problem with subsetting data for:", my_sample_type)
|
||||
# quit()
|
||||
#}
|
||||
|
||||
########################################################################
|
||||
# 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)
|
||||
write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
|
||||
376
asthma/flu_a_stats_unpaired_serum.R
Executable file
376
asthma/flu_a_stats_unpaired_serum.R
Executable file
|
|
@ -0,0 +1,376 @@
|
|||
#!/usr/bin/Rscript
|
||||
getwd()
|
||||
setwd("~/git/mosaic_2020/")
|
||||
getwd()
|
||||
############################################################
|
||||
# TASK: unpaired (time) analysis of mediators: serum
|
||||
############################################################
|
||||
my_sample_type = "serum"
|
||||
|
||||
#=============
|
||||
# Input
|
||||
#=============
|
||||
#source("data_extraction_mediators.R")
|
||||
source("plot_data_a.R")
|
||||
|
||||
# check: copd and asthma conflict
|
||||
table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
|
||||
|
||||
#=============
|
||||
# Output
|
||||
#=============
|
||||
outfile_name = paste0("flu_stats_time_unpaired_a_", my_sample_type, ".csv")
|
||||
flu_stats_time_unpaired = paste0(outdir_stats_a, outfile_name)
|
||||
flu_stats_time_unpaired
|
||||
|
||||
# quick checks
|
||||
table(wf_fp_serum$T1_resp_score)
|
||||
table(lf_fp_serum$mediator, lf_fp_npa$timepoint, lf_fp_npa$T1_resp_score)
|
||||
table(wf_fp_serum$obesity)
|
||||
|
||||
#===============================
|
||||
# data assignment for stats
|
||||
#================================
|
||||
wf = wf_fp_serum
|
||||
lf = lf_fp_serum
|
||||
|
||||
identical(wf, wf_fp_serum)
|
||||
identical(lf, lf_fp_serum)
|
||||
########################################################################
|
||||
# clear variables
|
||||
rm(lf_fp_sam, lf_fp_npa
|
||||
, wf_fp_sam, wf_fp_npa)
|
||||
########################################################################
|
||||
# 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))
|
||||
}
|
||||
|
||||
# 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("UNpaired stats for groups will be:", flu_stats_time_unpaired)
|
||||
write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
|
||||
244
severe/check_stats.R
Executable file
244
severe/check_stats.R
Executable file
|
|
@ -0,0 +1,244 @@
|
|||
#!/usr/bin/Rscript
|
||||
getwd()
|
||||
setwd("~/git/mosaic_2020/")
|
||||
getwd()
|
||||
########################################################################
|
||||
# TASK: check stats
|
||||
########################################################################
|
||||
# load libraries, packages and local imports
|
||||
source("Header_TT.R")
|
||||
|
||||
#=============
|
||||
# Input
|
||||
#=============
|
||||
source("data_extraction_mediators.R")
|
||||
#source("read_data.R")
|
||||
|
||||
########################################################################
|
||||
#=============
|
||||
# Data extraction for non-severe
|
||||
#=============
|
||||
table(fp_adults_ics$T1_resp_score)
|
||||
table(fp_adults_ics$obesity)
|
||||
|
||||
df = fp_adults_ics[fp_adults_ics$T1_resp_score==3,]
|
||||
table(df$obesity)
|
||||
|
||||
df_ob = df[df$obesity == 1,]
|
||||
df_nob = df[df$obesity == 0,]
|
||||
|
||||
rm(npa_lf, npa_wf, sam_lf, sam_wf, serum_lf, serum_wf, fp_adults_ics)
|
||||
|
||||
#============
|
||||
# npa
|
||||
#============
|
||||
# eotaxin: t1
|
||||
wilcox.test(df_ob$eotaxin_npa1, df_nob$eotaxin_npa1)
|
||||
wilcox_test(eotaxin_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$eotaxin_npa1)
|
||||
na_c = sum(is.na(foo$df.eotaxin_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$eotaxin_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$eotaxin_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
# eotaxin: t1
|
||||
wilcox.test(df_ob$eotaxin3_npa1, df_nob$eotaxin3_npa1)
|
||||
wilcox_test(eotaxin3_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$eotaxin3_npa1)
|
||||
na_c = sum(is.na(foo$df.eotaxin3_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$eotaxin3_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$eotaxin3_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
# il1: t1
|
||||
wilcox.test(df_ob$il1_npa1, df_nob$il1_npa1)
|
||||
wilcox_test(il1_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il1_npa1)
|
||||
na_c = sum(is.na(foo$df.il1_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il1_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il1_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
# il10: t1
|
||||
wilcox.test(df_ob$il10_npa1, df_nob$il10_npa1)
|
||||
wilcox_test(il10_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il10_npa1)
|
||||
na_c = sum(is.na(foo$df.il10_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il10_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il10_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
# il12p70: t1
|
||||
wilcox.test(df_ob$il12p70_npa1, df_nob$il12p70_npa1)
|
||||
wilcox_test(il12p70_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il12p70_npa1)
|
||||
na_c = sum(is.na(foo$df.il12p70_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il12p70_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il12p70_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
# il17: t1
|
||||
wilcox.test(df_ob$il17_npa1, df_nob$il17_npa1)
|
||||
wilcox_test(il17_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il17_npa1)
|
||||
na_c = sum(is.na(foo$df.il7_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il17_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il7_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
# il5: t1
|
||||
wilcox.test(df_ob$il5_npa1, df_nob$il5_npa1)
|
||||
wilcox_test(il5_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il5_npa1)
|
||||
na_c = sum(is.na(foo$df.il5_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il5_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il5_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
|
||||
# il6: t1
|
||||
wilcox.test(df_ob$il6_npa1, df_nob$il6_npa1)
|
||||
wilcox_test(il6_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il6_npa1)
|
||||
na_c = sum(is.na(foo$df.il6_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il6_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il6_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
# mcp4: t1
|
||||
wilcox.test(df_ob$mcp4_npa1, df_nob$mcp4_npa1)
|
||||
wilcox_test(mcp4_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$mcp4_npa1)
|
||||
na_c = sum(is.na(foo$df.mcp4_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$mcp4_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$mcp4_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
|
||||
# mip17: t1
|
||||
wilcox.test(df_ob$mip17_npa1, df_nob$mip17_npa1)
|
||||
wilcox_test(mip17_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$mip17_npa1)
|
||||
na_c = sum(is.na(foo$df.mip17_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$mip17_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$mip17_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
|
||||
# neopterin: t1
|
||||
wilcox.test(df_ob$neopterin_npa1, df_nob$neopterin_npa1)
|
||||
wilcox_test(neopterin_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$neopterin_npa1)
|
||||
na_c = sum(is.na(foo$df.neopterin_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$neopterin_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$neopterin_npa1, df$obesity, summary, na.rm = T)
|
||||
|
||||
# tnfr1: t1
|
||||
wilcox.test(df_ob$tnfr1_npa1, df_nob$tnfr1_npa1)
|
||||
wilcox_test(tnfr1_npa1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$neopterin_npa1)
|
||||
na_c = sum(is.na(foo$df.tnfr1_npa1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$tnfr1_npa1, df$obesity, median, na.rm = T)
|
||||
tapply(df$tnfr1_npa1 , df$obesity, summary, na.rm = T)
|
||||
|
||||
#============
|
||||
# sam
|
||||
#============
|
||||
# il5: t1
|
||||
wilcox.test(df_ob$il5_sam1, df_nob$il5_sam1)
|
||||
wilcox_test(il5_sam1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il5_sam1)
|
||||
na_c = sum(is.na(foo$df.il5_sam1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il5_sam1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il5_sam1 , df$obesity, summary, na.rm = T)
|
||||
|
||||
#============
|
||||
# serum
|
||||
#============
|
||||
# il17: t1
|
||||
wilcox.test(df_ob$il17_serum1, df_nob$il17_serum1)
|
||||
wilcox_test(il17_serum1 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il17_serum1)
|
||||
na_c = sum(is.na(foo$df.il17_serum1))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il17_serum1, df$obesity, median, na.rm = T)
|
||||
tapply(df$il17_serum1 , df$obesity, summary, na.rm = T)
|
||||
|
||||
|
||||
# il17: t3
|
||||
wilcox.test(df_ob$il17_serum3, df_nob$il17_serum3)
|
||||
wilcox_test(il17_serum3 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$il17_serum3)
|
||||
na_c = sum(is.na(foo$df.il17_serum3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$il17_serum3, df$obesity, median, na.rm = T)
|
||||
tapply(df$il17_serum3 , df$obesity, summary, na.rm = T)
|
||||
|
||||
# mcp4: t3
|
||||
wilcox.test(df_ob$mcp4_serum3, df_nob$mcp4_serum3)
|
||||
wilcox_test(mcp4_serum3 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$mcp4_serum3)
|
||||
na_c = sum(is.na(foo$df.mcp4_serum3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$mcp4_serum3, df$obesity, median, na.rm = T)
|
||||
tapply(df$mcp4_serum3, df$obesity, summary, na.rm = T)
|
||||
|
||||
|
||||
# eotaxin3: t3
|
||||
wilcox.test(df_ob$eotaxin3_serum3, df_nob$eotaxin3_serum3)
|
||||
wilcox_test(eotaxin3_serum3 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$eotaxin3_serum3)
|
||||
na_c = sum(is.na(foo$df.eotaxin3_serum3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$eotaxin3_serum3, df$obesity, median, na.rm = T)
|
||||
tapply(df$eotaxin3_serum3, df$obesity, summary, na.rm = T)
|
||||
|
||||
|
||||
# rantes: t3
|
||||
wilcox.test(df_ob$rantes_serum3, df_nob$rantes_serum3)
|
||||
wilcox_test(rantes_serum3 ~ obesity, data = df)
|
||||
|
||||
foo = data.frame(df$obesity, df$rantes_serum3)
|
||||
na_c = sum(is.na(foo$df.rantes_serum3))
|
||||
nrow(df) - na_c
|
||||
|
||||
tapply(df$rantes_serum3, df$obesity, median, na.rm = T)
|
||||
tapply(df$rantes_serum3, df$obesity, summary, na.rm = T)
|
||||
|
||||
386
severe/flu_s_stats_unpaired_npa.R
Executable file
386
severe/flu_s_stats_unpaired_npa.R
Executable file
|
|
@ -0,0 +1,386 @@
|
|||
#!/usr/bin/Rscript
|
||||
getwd()
|
||||
setwd("~/git/mosaic_2020/")
|
||||
getwd()
|
||||
############################################################
|
||||
# TASK: unpaired (time) analysis of mediators:
|
||||
# sample type: NPA
|
||||
# data: Flu positive adult patients
|
||||
# group: obesity
|
||||
############################################################
|
||||
my_sample_type = "npa"
|
||||
|
||||
#=============
|
||||
# Input
|
||||
#=============
|
||||
#source("data_extraction_mediators.R")
|
||||
source("plot_data_s.R")
|
||||
|
||||
# check: copd and asthma conflict
|
||||
table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
|
||||
|
||||
#=============
|
||||
# Output
|
||||
#=============
|
||||
outfile_name = paste0("flu_stats_time_unpaired_s_", my_sample_type, ".csv")
|
||||
flu_stats_time_unpaired = paste0(outdir_stats_s, outfile_name)
|
||||
flu_stats_time_unpaired
|
||||
|
||||
# quick checks
|
||||
table(wf_fp_npa$T1_resp_score)
|
||||
table(lf_fp_npa$mediator, lf_fp_npa$timepoint, lf_fp_npa$T1_resp_score)
|
||||
|
||||
#===============================
|
||||
# data assignment for stats
|
||||
#================================
|
||||
wf = wf_fp_npa
|
||||
lf = lf_fp_npa
|
||||
|
||||
identical(wf, wf_fp_npa)
|
||||
identical(lf, lf_fp_npa)
|
||||
########################################################################
|
||||
# clear variables
|
||||
rm(lf_fp_sam, lf_fp_serum
|
||||
, wf_fp_sam, wf_fp_serum)
|
||||
|
||||
# sanity checks
|
||||
table(lf$timepoint)
|
||||
length(unique(lf$mosaic))
|
||||
table(wf$obesity)
|
||||
|
||||
########################################################################
|
||||
# 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))
|
||||
}
|
||||
|
||||
# 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)
|
||||
write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
|
||||
395
severe/flu_s_stats_unpaired_sam.R
Executable file
395
severe/flu_s_stats_unpaired_sam.R
Executable file
|
|
@ -0,0 +1,395 @@
|
|||
#!/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_mediators.R")
|
||||
source("plot_data_s.R")
|
||||
|
||||
# check: copd and asthma conflict
|
||||
table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
|
||||
|
||||
#=============
|
||||
# Output
|
||||
#=============
|
||||
outfile_name = paste0("flu_stats_time_unpaired_s_", my_sample_type, ".csv")
|
||||
flu_stats_time_unpaired = paste0(outdir_stats_s, outfile_name)
|
||||
flu_stats_time_unpaired
|
||||
|
||||
# quick checks
|
||||
table(wf_fp_sam$T1_resp_score)
|
||||
table(lf_fp_sam$mediator, lf_fp_sam$timepoint, lf_fp_sam$T1_resp_score)
|
||||
table(wf$obesity)
|
||||
|
||||
#===============================
|
||||
# data assignment for stats
|
||||
#================================
|
||||
wf = wf_fp_sam
|
||||
lf = lf_fp_sam
|
||||
|
||||
identical(wf, wf_fp_sam)
|
||||
identical(lf, lf_fp_sam)
|
||||
|
||||
########################################################################
|
||||
# clear variables
|
||||
rm(lf_fp_npa, lf_fp_serum
|
||||
, wf_fp_npa, wf_fp_serum)
|
||||
|
||||
# sanity checks
|
||||
table(lf$timepoint)
|
||||
length(unique(lf$mosaic))
|
||||
|
||||
#if (table(lf$flustat) == table(sam_lf$flustat)[[2]]){
|
||||
# cat("Analysing Flu positive patients for:", my_sample_type)
|
||||
#}else{
|
||||
# cat("FAIL: problem with subsetting data for:", my_sample_type)
|
||||
# quit()
|
||||
#}
|
||||
|
||||
########################################################################
|
||||
# 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)
|
||||
write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
|
||||
376
severe/flu_s_stats_unpaired_serum.R
Executable file
376
severe/flu_s_stats_unpaired_serum.R
Executable file
|
|
@ -0,0 +1,376 @@
|
|||
#!/usr/bin/Rscript
|
||||
getwd()
|
||||
setwd("~/git/mosaic_2020/")
|
||||
getwd()
|
||||
############################################################
|
||||
# TASK: unpaired (time) analysis of mediators: serum
|
||||
############################################################
|
||||
my_sample_type = "serum"
|
||||
|
||||
#=============
|
||||
# Input
|
||||
#=============
|
||||
#source("data_extraction_mediators.R")
|
||||
source("plot_data_s.R")
|
||||
|
||||
# check: copd and asthma conflict
|
||||
table(fp_adults_ics$ia_exac_copd==1 & fp_adults_ics$asthma == 1)
|
||||
|
||||
#=============
|
||||
# Output
|
||||
#=============
|
||||
outfile_name = paste0("flu_stats_time_unpaired_s_", my_sample_type, ".csv")
|
||||
flu_stats_time_unpaired = paste0(outdir_stats_s, outfile_name)
|
||||
flu_stats_time_unpaired
|
||||
|
||||
# quick checks
|
||||
table(wf_fp_npa$T1_resp_score)
|
||||
table(lf_fp_npa$mediator, lf_fp_npa$timepoint, lf_fp_npa$T1_resp_score)
|
||||
table(wf$obesity)
|
||||
|
||||
#===============================
|
||||
# data assignment for stats
|
||||
#================================
|
||||
wf = wf_fp_serum
|
||||
lf = lf_fp_serum
|
||||
|
||||
identical(wf, wf_fp_serum)
|
||||
identical(lf, lf_fp_serum)
|
||||
########################################################################
|
||||
# clear variables
|
||||
rm(lf_fp_sam, lf_fp_npa
|
||||
, wf_fp_sam, wf_fp_npa)
|
||||
########################################################################
|
||||
# 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))
|
||||
}
|
||||
|
||||
# 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("UNpaired stats for groups will be:", flu_stats_time_unpaired)
|
||||
write.csv(combined_unpaired_stats_f, flu_stats_time_unpaired, row.names = FALSE)
|
||||
Loading…
Add table
Add a link
Reference in a new issue