added fd corrected p-values for ks stats
This commit is contained in:
parent
f5f1e388c3
commit
365c322953
3 changed files with 333 additions and 176 deletions
|
|
@ -18,38 +18,26 @@ source("~/git/LSHTM_analysis/scripts/plotting/plotting_colnames.R")
|
|||
# Output
|
||||
#=============
|
||||
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
|
||||
outdir_stats = "~/git/LSHTM_analysis/scripts/plotting/plotting_thesis/stats/"
|
||||
outdir_stats = paste0(outdir_images,"stats/")
|
||||
|
||||
# ks test by lineage
|
||||
#ks_lineage = paste0(outdir, "/KS_lineage_all_muts.csv")
|
||||
|
||||
###########################
|
||||
# Data for stats
|
||||
# you need merged_df2 or merged_df2_comp
|
||||
# you need df2 or df2_comp
|
||||
# since this is one-many relationship
|
||||
# i.e the same SNP can belong to multiple lineages
|
||||
# using the _comp dataset means
|
||||
# we lose some muts and at this level, we should use
|
||||
# as much info as available, hence use df with NA
|
||||
###########################
|
||||
|
||||
# REASSIGNMENT
|
||||
my_df = merged_df2
|
||||
|
||||
# delete variables not required
|
||||
rm(merged_df2, merged_df2_comp, merged_df3, merged_df3_comp)
|
||||
rm(merged_df2_lig, merged_df2_comp_lig, merged_df3_lig, merged_df3_comp_lig, my_df_u, my_df_u_lig)
|
||||
df2 = merged_df2[, colnames(merged_df2)%in%plotting_cols]
|
||||
|
||||
# quick checks
|
||||
colnames(my_df)
|
||||
str(my_df)
|
||||
colnames(df2)
|
||||
str(df2)
|
||||
|
||||
# Ensure correct data type in columns to plot: need to be factor
|
||||
#is.factor(my_df$lineage)
|
||||
#my_df$lineage = as.factor(my_df$lineage)
|
||||
#is.factor(my_df$lineage)
|
||||
########################################################################
|
||||
table(my_df$lineage); str(my_df$lineage)
|
||||
table(df2$lineage); str(df2$lineage)
|
||||
|
||||
# subset only lineages1-4
|
||||
sel_lineages = c("L1"
|
||||
|
|
@ -58,29 +46,22 @@ sel_lineages = c("L1"
|
|||
, "L4")
|
||||
|
||||
# subset selected lineages
|
||||
df_lin = subset(my_df, subset = lineage %in% sel_lineages)
|
||||
df_lin = subset(df2, subset = lineage %in% sel_lineages)
|
||||
|
||||
table(df_lin$lineage)
|
||||
table(df_lin$sensitivity)
|
||||
|
||||
table(df_lin$lineage, df_lin$sensitivity)
|
||||
|
||||
#==============
|
||||
# dr_muts_col
|
||||
#==============
|
||||
#table(df_lin$mutation_info); str(df_lin$mutation_info)
|
||||
#ensure lineage is a factor
|
||||
#str(df_lin$lineage); str(df_lin$lineage_labels)
|
||||
#df_lin$lineage = as.factor(df_lin$lineage)
|
||||
#df_lin$lineage_labels = as.factor(df_lin$lineage)
|
||||
table(df_lin$lineage); table(df_lin$lineage_labels)]
|
||||
|
||||
# subset df with dr muts only
|
||||
#df_lin_dr = subset(df_lin, mutation_info == dr_muts_col)
|
||||
#table(df_lin_dr$mutation_info)
|
||||
|
||||
#==============
|
||||
# other_muts_col
|
||||
#==============
|
||||
#df_lin_other = subset(df_lin, mutation_info == other_muts_col)
|
||||
#table(df_lin_other$mutation_info)
|
||||
|
||||
#=======================================================================
|
||||
#==============================
|
||||
# Stats for average stability
|
||||
#=============================
|
||||
# individual: CHECKS
|
||||
lin1 = df_lin[df_lin$lineage == "L1",]$avg_stability_scaled
|
||||
lin2 = df_lin[df_lin$lineage == "L2",]$avg_stability_scaled
|
||||
|
|
@ -89,14 +70,14 @@ lin4 = df_lin[df_lin$lineage == "L4",]$avg_stability_scaled
|
|||
|
||||
ks.test(lin1, lin4)
|
||||
|
||||
|
||||
ks.test(df_lin$avg_stability_scaled[df_lin$lineage == "L2"]
|
||||
, df_lin$avg_stability_scaled[df_lin$lineage == "L3"])
|
||||
ks.test(df_lin$avg_stability_scaled[df_lin$lineage == "L1"]
|
||||
, df_lin$avg_stability_scaled[df_lin$lineage == "L4"])
|
||||
|
||||
#=======================================================================
|
||||
my_lineages = levels(factor(df_lin$lineage)); my_lineages
|
||||
#=======================================================================
|
||||
# Loop
|
||||
#0 : < 2.2e-16
|
||||
#=====================
|
||||
# Lineage 1 comparisons
|
||||
#=====================
|
||||
|
|
@ -131,7 +112,7 @@ for (i in my_lineages_comp_l1){
|
|||
, df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value
|
||||
|
||||
# print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval))
|
||||
l1_df$comparison = lineage_comp
|
||||
l1_df$comparison = lineage_comp
|
||||
l1_df$method = ks_method
|
||||
l1_df$ks_statistic = ks_statistic[[1]]
|
||||
l1_df$ks_pvalue = ks_pvalue
|
||||
|
|
@ -142,6 +123,33 @@ for (i in my_lineages_comp_l1){
|
|||
ks_df_l1 = rbind(ks_df_l1,l1_df)
|
||||
|
||||
}
|
||||
|
||||
ks_df_l1
|
||||
|
||||
# adjusted p-value: bonferroni
|
||||
ks_df_l1$p_adj_bonferroni = p.adjust(ks_df_l1$ks_pvalue, method = "bonferroni")
|
||||
ks_df_l1$signif_bon = ks_df_l1$p_adj_bonferroni
|
||||
ks_df_l1 = dplyr::mutate(ks_df_l1
|
||||
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
||||
, signif_bon <=0.0001 ~ '****'
|
||||
, signif_bon <=0.001 ~ '***'
|
||||
, signif_bon <=0.01 ~ '**'
|
||||
, signif_bon <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
# adjusted p-value:fdr
|
||||
ks_df_l1$p_adj_fdr = p.adjust(ks_df_l1$ks_pvalue, method = "fdr")
|
||||
ks_df_l1$signif_fdr = ks_df_l1$p_adj_fdr
|
||||
ks_df_l1 = dplyr::mutate(ks_df_l1
|
||||
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
||||
, signif_fdr <=0.0001 ~ '****'
|
||||
, signif_fdr <=0.001 ~ '***'
|
||||
, signif_fdr <=0.01 ~ '**'
|
||||
, signif_fdr <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
ks_df_l1
|
||||
|
||||
#####################################################################
|
||||
|
||||
#=====================
|
||||
|
|
@ -190,6 +198,31 @@ for (i in my_lineages_comp_l2){
|
|||
|
||||
}
|
||||
|
||||
|
||||
# adjusted p-value: bonferroni
|
||||
ks_df_l2$p_adj_bonferroni = p.adjust(ks_df_l2$ks_pvalue, method = "bonferroni")
|
||||
ks_df_l2$signif_bon = ks_df_l2$p_adj_bonferroni
|
||||
ks_df_l2 = dplyr::mutate(ks_df_l2
|
||||
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
||||
, signif_bon <=0.0001 ~ '****'
|
||||
, signif_bon <=0.001 ~ '***'
|
||||
, signif_bon <=0.01 ~ '**'
|
||||
, signif_bon <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
# adjusted p-value:fdr
|
||||
ks_df_l2$p_adj_fdr = p.adjust(ks_df_l2$ks_pvalue, method = "fdr")
|
||||
ks_df_l2$signif_fdr = ks_df_l2$p_adj_fdr
|
||||
ks_df_l2 = dplyr::mutate(ks_df_l2
|
||||
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
||||
, signif_fdr <=0.0001 ~ '****'
|
||||
, signif_fdr <=0.001 ~ '***'
|
||||
, signif_fdr <=0.01 ~ '**'
|
||||
, signif_fdr <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
ks_df_l2
|
||||
|
||||
#=====================
|
||||
# Lineage 3 comparisons
|
||||
#=====================
|
||||
|
|
@ -234,11 +267,33 @@ for (i in my_lineages_comp_l3){
|
|||
ks_df_l3 = rbind(ks_df_l3, l3_df)
|
||||
|
||||
}
|
||||
######################################################################################
|
||||
|
||||
# adjusted p-value: bonferroni
|
||||
ks_df_l3$p_adj_bonferroni = p.adjust(ks_df_l3$ks_pvalue, method = "bonferroni")
|
||||
ks_df_l3$signif_bon = ks_df_l3$p_adj_bonferroni
|
||||
ks_df_l3 = dplyr::mutate(ks_df_l3
|
||||
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
||||
, signif_bon <=0.0001 ~ '****'
|
||||
, signif_bon <=0.001 ~ '***'
|
||||
, signif_bon <=0.01 ~ '**'
|
||||
, signif_bon <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
# adjusted p-value:fdr
|
||||
ks_df_l3$p_adj_fdr = p.adjust(ks_df_l3$ks_pvalue, method = "fdr")
|
||||
ks_df_l3$signif_fdr = ks_df_l3$p_adj_fdr
|
||||
ks_df_l3 = dplyr::mutate(ks_df_l3
|
||||
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
||||
, signif_fdr <=0.0001 ~ '****'
|
||||
, signif_fdr <=0.001 ~ '***'
|
||||
, signif_fdr <=0.01 ~ '**'
|
||||
, signif_fdr <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
ks_df_l3
|
||||
|
||||
####################################################################
|
||||
# combine all 4 ks_dfs
|
||||
# combine all 3 ks_dfs
|
||||
n_dfs = 3
|
||||
if ( all.equal(nrow(ks_df_l1), nrow(ks_df_l2), nrow(ks_df_l3)) &&
|
||||
all.equal(ncol(ks_df_l1), ncol(ks_df_l2), ncol(ks_df_l3)) ){
|
||||
|
|
@ -263,12 +318,9 @@ if ( all.equal(nrow(ks_df_l1), nrow(ks_df_l2), nrow(ks_df_l3)) &&
|
|||
, "\nCheck hardcoded value of n_dfs")
|
||||
}
|
||||
|
||||
#--------------
|
||||
# formatting
|
||||
#--------------
|
||||
# add total_n number
|
||||
#ks_df_combined$total_samples_analysed = nrow(df_lin)
|
||||
|
||||
#----------------------------
|
||||
# ADD extra cols: formatting
|
||||
#----------------------------
|
||||
# adding pvalue_signif
|
||||
ks_df_combined$pvalue_signif = ks_df_combined$ks_pvalue
|
||||
str(ks_df_combined$pvalue_signif)
|
||||
|
|
@ -322,9 +374,31 @@ overall_RS_df$ks_pvalue = ks_pvalue
|
|||
overall_RS_df$n_samples = n_samples_all
|
||||
overall_RS_df$n_samples_total= n_samples_total
|
||||
|
||||
#--------------
|
||||
# formatting
|
||||
#--------------
|
||||
#----------------------------
|
||||
# ADD extra cols
|
||||
#----------------------------
|
||||
# adjusted p-value: bonferroni
|
||||
overall_RS_df$p_adj_bonferroni = p.adjust(overall_RS_df$ks_pvalue, method = "bonferroni")
|
||||
overall_RS_df$signif_bon = overall_RS_df$p_adj_bonferroni
|
||||
overall_RS_df = dplyr::mutate(overall_RS_df
|
||||
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
||||
, signif_bon <=0.0001 ~ '****'
|
||||
, signif_bon <=0.001 ~ '***'
|
||||
, signif_bon <=0.01 ~ '**'
|
||||
, signif_bon <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
# adjusted p-value:fdr
|
||||
overall_RS_df$p_adj_fdr = p.adjust(overall_RS_df$ks_pvalue, method = "fdr")
|
||||
overall_RS_df$signif_fdr = overall_RS_df$p_adj_fdr
|
||||
overall_RS_df = dplyr::mutate(overall_RS_df
|
||||
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
||||
, signif_fdr <=0.0001 ~ '****'
|
||||
, signif_fdr <=0.001 ~ '***'
|
||||
, signif_fdr <=0.01 ~ '**'
|
||||
, signif_fdr <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
# unadjusted p-values
|
||||
overall_RS_df$pvalue_signif = overall_RS_df$ks_pvalue
|
||||
overall_RS_df = dplyr::mutate(overall_RS_df
|
||||
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
|
||||
|
|
@ -333,8 +407,6 @@ overall_RS_df = dplyr::mutate(overall_RS_df
|
|||
, pvalue_signif <=0.01 ~ '**'
|
||||
, pvalue_signif <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
overall_RS_df
|
||||
#####################################################################
|
||||
if (all(colnames(ks_df_combined_f) == colnames(overall_RS_df))){
|
||||
|
||||
|
|
@ -349,20 +421,12 @@ ks_df_combined_f2 = rbind(ks_df_combined_f, overall_RS_df)
|
|||
ks_df_combined_f2$ks_comp_type = "between_lineages"
|
||||
ks_df_combined_f2$gene_name = tolower(gene)
|
||||
|
||||
# #==============================
|
||||
# # write output file: KS test
|
||||
# #===============================
|
||||
# Out_lineage_bwL = paste0(outdir_stats
|
||||
# , tolower(gene)
|
||||
# , "_ks_lineage_bw.csv")
|
||||
#
|
||||
# cat("Output of KS test bt lineage:", Out_lineage_bwL)
|
||||
# write.csv(ks_df_combined_f2, Out_lineage_bwL, row.names = FALSE)
|
||||
ks_df_combined_f2
|
||||
|
||||
###########################################################################
|
||||
#=======================
|
||||
#=================================
|
||||
# Within lineage R vs S: MANUAL
|
||||
#=======================
|
||||
#=================================
|
||||
lin1 = df_lin[df_lin$lineage == my_lin1,]
|
||||
ks.test(lin1$avg_stability_scaled[lin1$sensitivity == "R"]
|
||||
, lin1$avg_stability_scaled[lin1$sensitivity == "S"])
|
||||
|
|
@ -426,35 +490,49 @@ for (i in c(my_lin1
|
|||
}
|
||||
all_within_lin_df
|
||||
|
||||
all_within_lin_df$pvalue_signif = all_within_lin_df$ks_pvalue
|
||||
str(all_within_lin_df$pvalue_signif)
|
||||
|
||||
#----------------------------
|
||||
# ADD extra cols: formatting
|
||||
#----------------------------
|
||||
# adjusted p-value: bonferroni
|
||||
all_within_lin_df$p_adj_bonferroni = p.adjust(all_within_lin_df$ks_pvalue, method = "bonferroni")
|
||||
all_within_lin_df$signif_bon = all_within_lin_df$p_adj_bonferroni
|
||||
all_within_lin_df = dplyr::mutate(all_within_lin_df
|
||||
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
|
||||
, pvalue_signif <=0.0001 ~ '****'
|
||||
, pvalue_signif <=0.001 ~ '***'
|
||||
, pvalue_signif <=0.01 ~ '**'
|
||||
, pvalue_signif <0.05 ~ '*'
|
||||
, signif_bon = case_when(signif_bon == 0.05 ~ "."
|
||||
, signif_bon <=0.0001 ~ '****'
|
||||
, signif_bon <=0.001 ~ '***'
|
||||
, signif_bon <=0.01 ~ '**'
|
||||
, signif_bon <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
all_within_lin_df
|
||||
# adjusted p-value:fdr
|
||||
all_within_lin_df$p_adj_fdr = p.adjust(all_within_lin_df$ks_pvalue, method = "fdr")
|
||||
all_within_lin_df$signif_fdr = all_within_lin_df$p_adj_fdr
|
||||
all_within_lin_df = dplyr::mutate(all_within_lin_df
|
||||
, signif_fdr = case_when(signif_fdr == 0.05 ~ "."
|
||||
, signif_fdr <=0.0001 ~ '****'
|
||||
, signif_fdr <=0.001 ~ '***'
|
||||
, signif_fdr <=0.01 ~ '**'
|
||||
, signif_fdr <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
# ADD extra cols
|
||||
# unadjusted p-value
|
||||
all_within_lin_df$pvalue_signif = all_within_lin_df$ks_pvalue
|
||||
str(all_within_lin_df$pvalue_signif)
|
||||
all_within_lin_df = dplyr::mutate(all_within_lin_df
|
||||
, pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
|
||||
, pvalue_signif <=0.0001 ~ '****'
|
||||
, pvalue_signif <=0.001 ~ '***'
|
||||
, pvalue_signif <=0.01 ~ '**'
|
||||
, pvalue_signif <0.05 ~ '*'
|
||||
, TRUE ~ 'ns'))
|
||||
|
||||
# ADD info cols
|
||||
all_within_lin_df$ks_comp_type = "within_lineages"
|
||||
all_within_lin_df$gene_name = tolower(gene)
|
||||
|
||||
# #--------------------
|
||||
# # write output file: KS test within grpup
|
||||
# #----------------------
|
||||
# Out_ks_withinL = paste0(outdir_stats
|
||||
# , tolower(gene)
|
||||
# , "_ks_lineage_within.csv")
|
||||
# cat("Output of KS test within lineage:",Out_ks_withinL )
|
||||
# write.csv(all_within_lin_df, Out_ks_withinL, row.names = FALSE)
|
||||
|
||||
##################################################################
|
||||
|
||||
if (all(colnames(ks_df_combined_f2) == colnames(Out_ks_withinL))){
|
||||
if (all(colnames(ks_df_combined_f2) == colnames(all_within_lin_df))){
|
||||
|
||||
cat("\nPASS:combining KS test results")
|
||||
}else{
|
||||
|
|
@ -469,5 +547,6 @@ ks_df_combined_all = rbind(ks_df_combined_f2, all_within_lin_df)
|
|||
Out_ks_all = paste0(outdir_stats
|
||||
, tolower(gene)
|
||||
, "_ks_lineage_all_comp.csv")
|
||||
|
||||
cat("Output of KS test all comparisons:", Out_ks_all )
|
||||
write.csv(ks_df_combined_all, Out_ks_all, row.names = FALSE)
|
||||
write.csv(ks_df_combined_all, Out_ks_all, row.names = FALSE)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue