diff --git a/scripts/plotting/plotting_thesis/rpob/rpob_ks_test_all_PS.R b/scripts/plotting/plotting_thesis/rpob/rpob_ks_test_all_PS.R new file mode 100755 index 0000000..8d18d99 --- /dev/null +++ b/scripts/plotting/plotting_thesis/rpob/rpob_ks_test_all_PS.R @@ -0,0 +1,545 @@ +#!/usr/bin/env Rscript +######################################################### +# TASK: KS test for PS/DUET lineage distributions +#======================================================================= +#!/usr/bin/env Rscript +#source("~/git/LSHTM_analysis/config/rpob.R") +# get plottting dfs +#source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") + +#============= +# Output +#============= +outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/") +outdir_stats = paste0(outdir_images,"stats/") + +# ks test by lineage +#ks_lineage = paste0(outdir, "/KS_lineage_all_muts.csv") + +########################### +# Data for stats +# you need df2 or df2_comp +# since this is one-many relationship +########################### + +# REASSIGNMENT +df2 = merged_df2[, colnames(merged_df2)%in%plotting_cols] + +# quick checks +colnames(df2) +str(df2) + +######################################################################## +table(df2$lineage); str(df2$lineage) + +# subset only lineages1-4 +sel_lineages = c("L1" + , "L2" + , "L3" + , "L4") + +# subset selected 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) + +#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) + +#============================== +# 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 +lin3 = df_lin[df_lin$lineage == "L3",]$avg_stability_scaled +lin4 = df_lin[df_lin$lineage == "L4",]$avg_stability_scaled + +ks.test(lin1, lin4) + +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 +#===================== +my_lin1 = "L1" +#my_lineages_comp_l1 = c("L2", "L3", "L4") +my_lineages_comp_l1 = my_lineages[-match(my_lin1, my_lineages)] + +ks_df_l1 = data.frame() + +for (i in my_lineages_comp_l1){ + cat(i) + l1_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA + , n_samples = NA + , n_samples_total = NA) + + lineage_comp = paste0(my_lin1, " vs ", i) + n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin1,]) + n_samples_i = nrow(df_lin[df_lin$lineage == i,]) + n_samples_all = paste0(my_lin1, "(", n_samples_lin, ")" + , ", " + , i, "(", n_samples_i, ")") + + n_samples_total = n_samples_lin + n_samples_i + + ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1] + , df_lin$avg_stability_scaled[df_lin$lineage == i])$method + + ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1] + , df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic + + ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin1] + , 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$method = ks_method + l1_df$ks_statistic = ks_statistic[[1]] + l1_df$ks_pvalue = ks_pvalue + l1_df$n_samples = n_samples_all + l1_df$n_samples_total= n_samples_total + + print(l1_df) + 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 + +##################################################################### + +#===================== +# Lineage 2 comparisons +#===================== +my_lin2 = "L2" +#my_lineages_comp_l2 = c("L1", lineage3", "L4") +my_lineages_comp_l2 = my_lineages[-match(my_lin2, my_lineages)] + +ks_df_l2 = data.frame() + +for (i in my_lineages_comp_l2){ + cat(i) + l2_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA + , n_samples = NA + , n_samples_total = NA) + + lineage_comp = paste0(my_lin2, " vs ", i) + n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin2,]) + n_samples_i = nrow(df_lin[df_lin$lineage == i,]) + n_samples_all = paste0(my_lin2, "(", n_samples_lin, ")" + , ", " + , i, "(", n_samples_i, ")") + + n_samples_total = n_samples_lin + n_samples_i + + ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2] + , df_lin$avg_stability_scaled[df_lin$lineage == i])$method + + ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2] + , df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic + + ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin2] + , df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value + + # print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval)) + l2_df$comparison = lineage_comp + l2_df$method = ks_method + l2_df$ks_statistic = ks_statistic[[1]] + l2_df$ks_pvalue = ks_pvalue + l2_df$n_samples = n_samples_all + l2_df$n_samples_total = n_samples_total + + print(l2_df) + ks_df_l2 = rbind(ks_df_l2, l2_df) + +} + + +# 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 +#===================== +my_lin3 = "L3" +#my_lineages_comp_l3 = c("L1", lineage2", "L4") +my_lineages_comp_l3 = my_lineages[-match(my_lin3, my_lineages)] + +ks_df_l3 = data.frame() + +for (i in my_lineages_comp_l3){ + cat(i) + l3_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA + , n_samples = NA + , n_samples_total = NA) + + lineage_comp = paste0(my_lin3, " vs ", i) + n_samples_lin = nrow(df_lin[df_lin$lineage == my_lin3,]) + n_samples_i = nrow(df_lin[df_lin$lineage == i,]) + n_samples_all = paste0(my_lin3, "(", n_samples_lin, ")" + , ", " + , i, "(", n_samples_i, ")") + n_samples_total = n_samples_lin + n_samples_i + + ks_method = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3] + , df_lin$avg_stability_scaled[df_lin$lineage == i])$method + + ks_statistic = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3] + , df_lin$avg_stability_scaled[df_lin$lineage == i])$statistic + + ks_pvalue = ks.test(df_lin$avg_stability_scaled[df_lin$lineage == my_lin3] + , df_lin$avg_stability_scaled[df_lin$lineage == i])$p.value + + # print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval)) + l3_df$comparison = lineage_comp + l3_df$method = ks_method + l3_df$ks_statistic = ks_statistic[[1]] + l3_df$ks_pvalue = ks_pvalue + l3_df$n_samples = n_samples_all + l3_df$n_samples_total = n_samples_total + + print(l3_df) + 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 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)) ){ + cat("\nPASS: Calculating expected rows and cols for sanity checks on combined_dfs") + expected_rows = nrow(ks_df_l1) * n_dfs + expected_cols = ncol(ks_df_l1) + ks_df_combined = rbind(ks_df_l1, ks_df_l2, ks_df_l3) + if ( nrow(ks_df_combined) == expected_rows && ncol(ks_df_combined) == expected_cols ){ + cat("\nPASS: combined df successfully created" + , "\nnrow combined_df:", nrow(ks_df_combined) + , "\nncol combined_df:", ncol(ks_df_combined)) + } + else{ + cat("\nFAIL: Dim mismatch" + , "\nExpected rows:", expected_rows + , "\nGot:", nrow(ks_df_combined) + , "\nExpected cols:", expected_cols + , "\nGot:", ncol(ks_df_combined)) + } +}else{ + cat("\nFAIL: Could not generate expected_rows and/or expected_cols" + , "\nCheck hardcoded value of n_dfs") +} + +#---------------------------- +# ADD extra cols: formatting +#---------------------------- +# adding pvalue_signif +ks_df_combined$pvalue_signif = ks_df_combined$ks_pvalue +str(ks_df_combined$pvalue_signif) + +ks_df_combined = dplyr::mutate(ks_df_combined + , 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')) + +# Remove duplicates +rows_to_remove = c("L2 vs L1", "L3 vs L1", "L3 vs L2") + +ks_df_combined_f = ks_df_combined[-match(rows_to_remove, ks_df_combined$comparison),] +################################################################################# +#================================ +# R vs S distribution: Overall +#================================= +df_lin_R = df_lin[df_lin$sensitivity == "R",] +df_lin_S = df_lin[df_lin$sensitivity == "S",] + +overall_RS_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA + , n_samples = NA + , n_samples_total = NA) + +comp_type = "overall R vs S distributions" + +n_samples_R = nrow(df_lin_R) +n_samples_S = nrow(df_lin_S) +n_samples_all = paste0("R(n=", n_samples_R, ")" , " ", "S(n=", n_samples_S, ")") +n_samples_total = n_samples_R+n_samples_S + +ks.test(df_lin_R$avg_stability_scaled + , df_lin_S$avg_stability_scaled) + +ks_method = ks.test(df_lin_R$avg_stability_scaled + , df_lin_S$avg_stability_scaled)$method + +ks_statistic = ks.test(df_lin_R$avg_stability_scaled + , df_lin_S$avg_stability_scaled)$statistic + +ks_pvalue = ks.test(df_lin_R$avg_stability_scaled + , df_lin_S$avg_stability_scaled)$p.value + +overall_RS_df$comparison = comp_type +overall_RS_df$method = ks_method +overall_RS_df$ks_statistic = ks_statistic +overall_RS_df$ks_pvalue = ks_pvalue +overall_RS_df$n_samples = n_samples_all +overall_RS_df$n_samples_total= n_samples_total + +#---------------------------- +# 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 ~ "." + , pvalue_signif <=0.0001 ~ '****' + , pvalue_signif <=0.001 ~ '***' + , pvalue_signif <=0.01 ~ '**' + , pvalue_signif <0.05 ~ '*' + , TRUE ~ 'ns')) +##################################################################### +if (all(colnames(ks_df_combined_f) == colnames(overall_RS_df))){ + + cat("\nPASS:combining KS test results") +}else{ + stop("\nAbort: Cannot combine results for KS test") +} + +ks_df_combined_f2 = rbind(ks_df_combined_f, overall_RS_df) + +# ADD extra cols +ks_df_combined_f2$ks_comp_type = "between_lineages" +ks_df_combined_f2$gene_name = tolower(gene) + +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"]) + +lin2 = df_lin[df_lin$lineage == my_lin2,] +ks.test(lin2$avg_stability_scaled[lin2$sensitivity == "R"] + , lin2$avg_stability_scaled[lin2$sensitivity == "S"]) + + +lin3 = df_lin[df_lin$lineage == "L3",] +ks.test(lin3$avg_stability_scaled[lin3$sensitivity == "R"] + , lin3$avg_stability_scaled[lin3$sensitivity == "S"]) + +lin4 = df_lin[df_lin$lineage == "L4",] +ks.test(lin4$avg_stability_scaled[lin4$sensitivity == "R"] + , lin4$avg_stability_scaled[lin4$sensitivity == "S"]) +################################################################### +#======================= +# Within lineage R vs S: LOOP +#======================= +all_within_lin_df = data.frame() + +for (i in c(my_lin1 + , my_lin2 + , my_lin3 + )){ + #cat(i) + + within_lin_df = data.frame(comparison = NA, method = NA, ks_statistic = NA, ks_pvalue = NA + , n_samples = NA + , n_samples_total = NA) + + within_lineage_comp = paste0(i, ": R vs S") + my_lin_df = df_lin[df_lin$lineage == i,] + + lin_R_n = nrow(my_lin_df[my_lin_df$sensitivity == "R",]) + lin_S_n = nrow(my_lin_df[my_lin_df$sensitivity == "S",]) + lin_total_n = lin_R_n+lin_S_n + n_samples_all = paste0("R(n=", lin_R_n, ")" , ", ", "S(n=", lin_S_n, ")") + + ks_method = ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"] + , my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$method + + ks_statistic = ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"] + , my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$statistic + + ks_pvalue =ks.test(my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "R"] + , my_lin_df$avg_stability_scaled[my_lin_df$sensitivity == "S"])$p.value + + # print(c(lineage_comp, ks_method, ks_statistic[[1]], ks_pval)) + within_lin_df$comparison = within_lineage_comp + within_lin_df$method = ks_method + within_lin_df$ks_statistic = ks_statistic[[1]] + within_lin_df$ks_pvalue = ks_pvalue + within_lin_df$n_samples = n_samples_all + within_lin_df$n_samples_total=lin_total_n + + print(my_lin_df) + all_within_lin_df = rbind(all_within_lin_df, within_lin_df) + +} +all_within_lin_df + +#---------------------------- +# 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 + , 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 +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')) + +# 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) + +################################################################## + +if (all(colnames(ks_df_combined_f2) == colnames(all_within_lin_df))){ + + cat("\nPASS:combining KS test results") +}else{ + stop("\nAbort: Cannot combine results for KS test") +} + +ks_df_combined_all = rbind(ks_df_combined_f2, all_within_lin_df) + +#-------------------- +# write output file: KS test within grpup +#---------------------- +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)