#!/usr/bin/env Rscript ######################################################### # TASK: KS test for PS/DUET lineage distributions #======================================================================= #!/usr/bin/env Rscript #source("~/git/LSHTM_analysis/config/alr.R") source("~/git/LSHTM_analysis/config/embb.R") #source("~/git/LSHTM_analysis/config/katg.R") #source("~/git/LSHTM_analysis/config/gid.R") #source("~/git/LSHTM_analysis/config/pnca.R") #source("~/git/LSHTM_analysis/config/rpob.R") # get plottting dfs source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") source("~/git/LSHTM_analysis/scripts/plotting/plotting_colnames.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)