From c1441dc0d6603375cede6ae6f03be713212cf428 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Tue, 21 Feb 2023 20:41:59 +0000 Subject: [PATCH] import thesis plots from Misc --- .../embb/embb_ORandSNP_results.R | 4 +- .../embb/embb_lineage_diff_sensitivities.R | 4 +- .../pnca/basic_barplots_pnca_layout.R | 4 +- .../plotting_thesis/thesis-figs-alr.R | 378 +++++++++++++++++ .../plotting_thesis/thesis-figs-embb.R | 266 ++++++++++++ .../plotting_thesis/thesis-figs-gid.R | 378 +++++++++++++++++ .../plotting_thesis/thesis-figs-katg.R | 383 +++++++++++++++++ .../plotting_thesis/thesis-figs-pnca.R | 271 ++++++++++++ .../plotting_thesis/thesis-figs-rpob.R | 386 ++++++++++++++++++ 9 files changed, 2068 insertions(+), 6 deletions(-) create mode 100644 scripts/plotting/plotting_thesis/thesis-figs-alr.R create mode 100644 scripts/plotting/plotting_thesis/thesis-figs-embb.R create mode 100644 scripts/plotting/plotting_thesis/thesis-figs-gid.R create mode 100644 scripts/plotting/plotting_thesis/thesis-figs-katg.R create mode 100644 scripts/plotting/plotting_thesis/thesis-figs-pnca.R create mode 100644 scripts/plotting/plotting_thesis/thesis-figs-rpob.R diff --git a/scripts/plotting/plotting_thesis/embb/embb_ORandSNP_results.R b/scripts/plotting/plotting_thesis/embb/embb_ORandSNP_results.R index e563807..26145de 100644 --- a/scripts/plotting/plotting_thesis/embb/embb_ORandSNP_results.R +++ b/scripts/plotting/plotting_thesis/embb/embb_ORandSNP_results.R @@ -1,6 +1,6 @@ #!/usr/bin/env Rscript -#source("~/git/LSHTM_analysis/config/katg.R") -#source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +source("~/git/LSHTM_analysis/config/embb.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") #======= # output diff --git a/scripts/plotting/plotting_thesis/embb/embb_lineage_diff_sensitivities.R b/scripts/plotting/plotting_thesis/embb/embb_lineage_diff_sensitivities.R index 732d7ae..b7e4b91 100644 --- a/scripts/plotting/plotting_thesis/embb/embb_lineage_diff_sensitivities.R +++ b/scripts/plotting/plotting_thesis/embb/embb_lineage_diff_sensitivities.R @@ -1,8 +1,8 @@ #============= # Data: Input #============== -#source("~/git/LSHTM_analysis/config/embb.R") -source("~/git/LSHTM_analysis/config/gid.R") +source("~/git/LSHTM_analysis/config/embb.R") +#source("~/git/LSHTM_analysis/config/gid.R") source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") diff --git a/scripts/plotting/plotting_thesis/pnca/basic_barplots_pnca_layout.R b/scripts/plotting/plotting_thesis/pnca/basic_barplots_pnca_layout.R index ceeeb8b..5b112c0 100644 --- a/scripts/plotting/plotting_thesis/pnca/basic_barplots_pnca_layout.R +++ b/scripts/plotting/plotting_thesis/pnca/basic_barplots_pnca_layout.R @@ -1,8 +1,8 @@ #============= # Data: Input #============== -#source("~/git/LSHTM_analysis/config/pnca.R") -#source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +source("~/git/LSHTM_analysis/config/pnca.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") source("~/git/LSHTM_analysis/scripts/plotting/plotting_thesis/pnca/basic_barplots_pnca.R") source("~/git/LSHTM_analysis/scripts/plotting/plotting_thesis/pnca/pe_sens_site_count_pnca.R") diff --git a/scripts/plotting/plotting_thesis/thesis-figs-alr.R b/scripts/plotting/plotting_thesis/thesis-figs-alr.R new file mode 100644 index 0000000..b85c8dd --- /dev/null +++ b/scripts/plotting/plotting_thesis/thesis-figs-alr.R @@ -0,0 +1,378 @@ +# +# small_df3 = merged_df3[merged_df3['position'] <= 670 & merged_df3['position'] >= 300,] +# +# small_df3 = merged_df3[merged_df3['position'] <= 451 & merged_df3['position'] >= 300,] +# +# +# mutable_df3 = cbind(small_df3) +# plot_df = cbind(small_df3) + +source("~/git/LSHTM_analysis/config/alr.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +outdir="/home/pub/Work/LSHTM/Thesis_Plots/" + +# add "pos_count" position count column +merged_df3=merged_df3 %>% dplyr::add_count(position) +merged_df3$pos_count=merged_df3$n +merged_df3$n=NULL + +# re-sort the dataframe according to position count +sorted_df = cbind(merged_df3) +sorted_df = sorted_df %>% arrange(pos_count) + +# determine some split points + +# further split up our sorted dataframe into low/mid/high for stability plotting +#snpcount_high_df = sorted_df[mid_snpcount:nrow(sorted_df),] +#snpcount_mid_df = sorted_df[low_snpcount:mid_snpcount,] +#snpcount_low_df = sorted_df[1:low_snpcount,] + +snpcount_high_df = sorted_df[sorted_df$pos_count>1,] +snpcount_low_df = sorted_df[sorted_df$pos_count<2,] + +low_half_snpcount = round(nrow(snpcount_low_df)/2) + +snpcount_low_low_df = snpcount_low_df[1:low_half_snpcount,] +snpcount_low_high_df = snpcount_low_df[low_half_snpcount:nrow(snpcount_low_df),] + +#### Put sorted-by-pos_count bp_stability_hmaps here #### +#...if they look better :-) + +length(merged_df3$position) + +# manually picked numbers due to irregular position distribution +low = 170 +centre=300 +#low=round(max(merged_df3$position)/3) +#centre=low*2 + +# Generate heat bar for whole DF +heatbar_df = generate_distance_colour_map(merged_df3) + +heat_bar_legend = get_legend( + ggplot(merged_df3, aes(factor(ligand_distance), y=0)) + + theme(legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = ligand_distance, + + ) + ) + + scale_colour_gradientn( + name="Ligand\nDistance (Å)", + colours = magma(5, direction = -1)) +) +# Position Tile Legend +# STR #00ff00 +tile_colour_map = c("DCS" = "#00ff00" + ,"PLP" = "#000080") + +tile_legend=get_legend( + + ggplot(tile_map, aes(factor(tile),y=0 + , colour=tile_colour + , fill=tile_colour))+ + geom_tile() + + theme(legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt")) + + scale_colour_manual(name=NULL + #, values = tile_map$tile_colour + , values=tile_colour_map) + + scale_fill_manual(name=NULL + #,values=tile_map$tile_colour + , values = tile_colour_map) +) + +# SAV Tile Red/Blue Legend + +temp_snp_df=merged_df3%>%dplyr::add_count(position) + +snp_tile_legend = get_legend( + ggplot(temp_snp_df, aes(factor(avg_stability_scaled), y=0)) + + theme(#legend.title = element_blank(), + legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = avg_stability_scaled, + + ) + ) + + scale_colour_gradientn( + name="Protein Stability\n(ΔΔG Kcal/mol)", colours = rev( + c("#007ba4", # rev() because FUCK + "#00bfc4", + "#ffffff", + "#ffd9d0", + "#f8766d", + "#ae3121") + ) + ) +) +# Logop Chemistry Colours +temp_logo_p_snp=cbind(merged_df3) +tab_mt = table(temp_logo_p_snp[["mutant_type"]], temp_logo_p_snp[["position"]]) +tab_mt = unclass(tab_mt) +logoplot_chemistry_legend=get_legend( + ggplot() + + geom_logo( + tab_mt + , method = 'custom' + , col_scheme = "chemistry" + , seq_type = 'aa') + theme( + legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt") + ) +) + +stab=function(...){bp_stability_hmap( + small_df3, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "", + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4], + y_max_override=max(merged_df3$pos_count), + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + ... +) +} + + + +# LOW +small_df3 = merged_df3[merged_df3['position'] < low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting LOW Rows: ", nrow(small_df3))) + + +logo_p_snps_low = LogoPlotSnps( + mutable_df3, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +# needs to be: active_aa_pos AND odds ratio >10 +logo_or_low = LogoPlotCustomH( + small_df3, + y_axis_increment=10, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + rm_empty_y=FALSE, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_low = stab() + + +# MID +small_df3 = merged_df3[merged_df3['position'] <= centre & merged_df3['position'] >= low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting MID Rows: ", nrow(small_df3))) + +logo_p_snps_centre = LogoPlotSnps( + mutable_df3, + y_lab = "", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +logo_or_centre = LogoPlotCustomH( + small_df3, y_axis_increment=10, + rm_empty_y=FALSE, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_centre = stab() + +# HIGH +small_df3 = merged_df3[merged_df3['position'] > centre,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting HIGH Rows: ", nrow(small_df3))) + + +logo_p_snps_high = LogoPlotSnps( + mutable_df3, + y_lab="" , + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +logo_or_high = LogoPlotCustomH( + small_df3, + y_axis_increment=10, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + rm_empty_y=FALSE, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + + +stability_high = stab() +#### "sorted" stability plots #### + +# do stability plots for later assembly +stability_high_reorder = bp_stability_hmap( + snpcount_high_df, + reorder_position = T, + p_title = NULL, + my_ylab = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] + +) +stability_centre_reorder = bp_stability_hmap( + snpcount_low_high_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "", + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_low_reorder = bp_stability_hmap( + snpcount_low_low_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = NULL, + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] + +) + +png(paste0(outdir,"alr_average_stability-sorted.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high_reorder, stability_centre_reorder, stability_low_reorder,ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"alr_average_stability-unsorted.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high, stability_centre, stability_low,ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"alr_snp.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_p_snps_low, logo_p_snps_centre, logo_p_snps_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"alr_or.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_or_low, logo_or_centre, logo_or_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() diff --git a/scripts/plotting/plotting_thesis/thesis-figs-embb.R b/scripts/plotting/plotting_thesis/thesis-figs-embb.R new file mode 100644 index 0000000..8d6c9c7 --- /dev/null +++ b/scripts/plotting/plotting_thesis/thesis-figs-embb.R @@ -0,0 +1,266 @@ +# +# small_df3 = merged_df3[merged_df3['position'] <= 670 & merged_df3['position'] >= 300,] +# +# small_df3 = merged_df3[merged_df3['position'] <= 451 & merged_df3['position'] >= 300,] +# +# +# mutable_df3 = cbind(small_df3) +# plot_df = cbind(small_df3) + +source("~/git/LSHTM_analysis/config/embb.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +outdir="/home/pub/Work/LSHTM/Thesis_Plots/" + +# add "pos_count" position count column +merged_df3=merged_df3 %>% dplyr::add_count(position) +merged_df3$pos_count=merged_df3$n +merged_df3$n=NULL + +# re-sort the dataframe according to position count +sorted_df = cbind(merged_df3) +sorted_df = sorted_df %>% arrange(pos_count) + +# determine some split points + +# further split up our sorted dataframe into low/mid/high for stability plotting +# snpcount_high_df = sorted_df[mid_snpcount:nrow(sorted_df),] +# snpcount_mid_df = sorted_df[low_snpcount:mid_snpcount,] +# snpcount_low_df = sorted_df[1:low_snpcount,] + +snpcount_high_df = sorted_df[sorted_df$pos_count>1,] +snpcount_low_df = sorted_df[sorted_df$pos_count<2,] + +low_half_snpcount = round(nrow(snpcount_low_df)/2) + +snpcount_low_low_df = snpcount_low_df[1:low_half_snpcount,] +snpcount_low_high_df = snpcount_low_df[low_half_snpcount:nrow(snpcount_low_df),] + + +# do stability plots for later assembly +stability_high = bp_stability_hmap(snpcount_high_df, + reorder_position = T, + p_title = NULL, + my_ylab = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3 +) +stability_centre = bp_stability_hmap(snpcount_low_high_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "SAV Count", + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3 +) +stability_low = bp_stability_hmap(snpcount_low_low_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = NULL, + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3 +) + +low=round(max(merged_df3$position)/3) +centre=low*2 + +# Generate heat bar for whole DF +heatbar_df = generate_distance_colour_map(merged_df3) + +heat_bar_legend = get_legend( + ggplot(merged_df3, aes(factor(ligand_distance), y=0)) + + theme(legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = ligand_distance, + + ) + ) + + scale_colour_gradientn(name="Ligand\nDistance (Å)", colours = magma(5, direction = -1)) +) +# Position Tile Legend + +tile_colour_map = c("EMB" = "green" + ,"DPA" = "darkslategrey" + , "CDL" = "navyblue" + , "Ca" = "purple") + +tile_legend=get_legend( + + ggplot(tile_map, aes(factor(tile),y=0 + , colour=tile_colour + , fill=tile_colour))+ + geom_tile() + + theme(legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt")) + + scale_colour_manual(name=NULL + #, values = tile_map$tile_colour + , values=tile_colour_map) + + scale_fill_manual(name=NULL + #,values=tile_map$tile_colour + , values = tile_colour_map) +) + +# SAV Tile Red/Blue Legend + +temp_snp_df=merged_df3%>%dplyr::add_count(position) + +snp_tile_legend = get_legend( + ggplot(temp_snp_df, aes(factor(avg_stability_scaled), y=0)) + + theme(#legend.title = element_blank(), + legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = avg_stability_scaled, + + ) + ) + + scale_colour_gradientn(name="Protein Stability\n(ΔΔG Kcal/mol)", colours = rev(c("#007ba4", # rev() because FUCK + "#00bfc4", + "#ffffff", + "#ffd9d0", + "#f8766d", + "#ae3121"))) +) +# Logop Chemistry Colours +temp_logo_p_snp=cbind(merged_df3) +tab_mt = table(temp_logo_p_snp[["mutant_type"]], temp_logo_p_snp[["position"]]) +tab_mt = unclass(tab_mt) +logoplot_chemistry_legend=get_legend( + ggplot() + + geom_logo(tab_mt + , method = 'custom' + , col_scheme = "chemistry" + , seq_type = 'aa') + theme(legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt")) +) + +# LOW +small_df3 = merged_df3[merged_df3['position'] < low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting LOW Rows: ", nrow(small_df3))) + + +logo_p_snps_low = LogoPlotSnps(mutable_df3, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3 +) +# needs to be: active_aa_pos AND odds ratio >10 +logo_or_low = LogoPlotCustomH(small_df3, + y_axis_increment=10, + y_lab="", + rm_empty_y=T, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3 +) + + +# MID +small_df3 = merged_df3[merged_df3['position'] <= centre & merged_df3['position'] >= low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting MID Rows: ", nrow(small_df3))) +logo_p_snps_centre = LogoPlotSnps(mutable_df3, y_lab = "SAV Count") + +logo_or_centre = LogoPlotCustomH(small_df3, + y_axis_increment=10, + rm_empty_y=T, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3 +) + +# HIGH +small_df3 = merged_df3[merged_df3['position'] > centre,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting HIGH Rows: ", nrow(small_df3))) + + +logo_p_snps_high = LogoPlotSnps(mutable_df3, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3 +) + +logo_or_high = LogoPlotCustomH(small_df3, + y_axis_increment=10, + y_lab="", + rm_empty_y=T, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3 +) + +png(paste0(outdir,"embb_average_stability.png"), height = 8.25, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high, stability_centre, stability_low,ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"embb_snp.png"), height = 8.25, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_p_snps_low, logo_p_snps_centre, logo_p_snps_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"embb_or.png"), height = 8.25, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_or_low, logo_or_centre, logo_or_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() diff --git a/scripts/plotting/plotting_thesis/thesis-figs-gid.R b/scripts/plotting/plotting_thesis/thesis-figs-gid.R new file mode 100644 index 0000000..1a781fb --- /dev/null +++ b/scripts/plotting/plotting_thesis/thesis-figs-gid.R @@ -0,0 +1,378 @@ +# +# small_df3 = merged_df3[merged_df3['position'] <= 670 & merged_df3['position'] >= 300,] +# +# small_df3 = merged_df3[merged_df3['position'] <= 451 & merged_df3['position'] >= 300,] +# +# +# mutable_df3 = cbind(small_df3) +# plot_df = cbind(small_df3) + +source("~/git/LSHTM_analysis/config/gid.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +outdir="/home/pub/Work/LSHTM/Thesis_Plots/" + +# add "pos_count" position count column +merged_df3=merged_df3 %>% dplyr::add_count(position) +merged_df3$pos_count=merged_df3$n +merged_df3$n=NULL + +# re-sort the dataframe according to position count +sorted_df = cbind(merged_df3) +sorted_df = sorted_df %>% arrange(pos_count) + +# determine some split points + +# further split up our sorted dataframe into low/mid/high for stability plotting +#snpcount_high_df = sorted_df[mid_snpcount:nrow(sorted_df),] +#snpcount_mid_df = sorted_df[low_snpcount:mid_snpcount,] +#snpcount_low_df = sorted_df[1:low_snpcount,] + +snpcount_high_df = sorted_df[sorted_df$pos_count>1,] +snpcount_low_df = sorted_df[sorted_df$pos_count<2,] + +low_half_snpcount = round(nrow(snpcount_low_df)/2) + +snpcount_low_low_df = snpcount_low_df[1:low_half_snpcount,] +snpcount_low_high_df = snpcount_low_df[low_half_snpcount:nrow(snpcount_low_df),] + +#### Put sorted-by-pos_count bp_stability_hmaps here #### +#...if they look better :-) + + +low=round(max(merged_df3$position)/3) +centre=low*2 + +# Generate heat bar for whole DF +heatbar_df = generate_distance_colour_map(merged_df3) + +heat_bar_legend = get_legend( + ggplot(merged_df3, aes(factor(ligand_distance), y=0)) + + theme(legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = ligand_distance, + + ) + ) + + scale_colour_gradientn( + name="Ligand\nDistance (Å)", + colours = magma(5, direction = -1)) +) +# Position Tile Legend +# STR #00ff00 +tile_colour_map = c("STR" = "green", + "SAM" = "#2f4f4f" + ,"RNA" = "#ff1493" + ,"AMP" = "#000080") + +tile_legend=get_legend( + + ggplot(tile_map, aes(factor(tile),y=0 + , colour=tile_colour + , fill=tile_colour))+ + geom_tile() + + theme(legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt")) + + scale_colour_manual(name=NULL + #, values = tile_map$tile_colour + , values=tile_colour_map) + + scale_fill_manual(name=NULL + #,values=tile_map$tile_colour + , values = tile_colour_map) +) + +# SAV Tile Red/Blue Legend + +temp_snp_df=merged_df3%>%dplyr::add_count(position) + +snp_tile_legend = get_legend( + ggplot(temp_snp_df, aes(factor(avg_stability_scaled), y=0)) + + theme(#legend.title = element_blank(), + legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = avg_stability_scaled, + + ) + ) + + scale_colour_gradientn( + name="Protein Stability\n(ΔΔG Kcal/mol)", colours = rev( + c("#007ba4", # rev() because FUCK + "#00bfc4", + "#ffffff", + "#ffd9d0", + "#f8766d", + "#ae3121") + ) + ) +) +# Logop Chemistry Colours +temp_logo_p_snp=cbind(merged_df3) +tab_mt = table(temp_logo_p_snp[["mutant_type"]], temp_logo_p_snp[["position"]]) +tab_mt = unclass(tab_mt) +logoplot_chemistry_legend=get_legend( + ggplot() + + geom_logo( + tab_mt + , method = 'custom' + , col_scheme = "chemistry" + , seq_type = 'aa') + theme( + legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt") + ) +) + +stab=function(...){bp_stability_hmap( + small_df3, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "", + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4], + y_max_override=max(merged_df3$pos_count), + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + ... +) +} + + + +# LOW +small_df3 = merged_df3[merged_df3['position'] < low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting LOW Rows: ", nrow(small_df3))) + + +logo_p_snps_low = LogoPlotSnps( + mutable_df3, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +# needs to be: active_aa_pos AND odds ratio >10 +logo_or_low = LogoPlotCustomH( + small_df3, + y_axis_increment=10, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + rm_empty_y=T, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + + + +stability_low = stab() + + +# MID +small_df3 = merged_df3[merged_df3['position'] <= centre & merged_df3['position'] >= low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting MID Rows: ", nrow(small_df3))) +logo_p_snps_centre = LogoPlotSnps( + mutable_df3, + y_lab = "", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +logo_or_centre = LogoPlotCustomH( + small_df3, y_axis_increment=10, + rm_empty_y=T, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_centre = stab() + +# HIGH +small_df3 = merged_df3[merged_df3['position'] > centre,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting HIGH Rows: ", nrow(small_df3))) + + +logo_p_snps_high = LogoPlotSnps( + mutable_df3, + y_lab="" , + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +logo_or_high = LogoPlotCustomH( + small_df3, + y_axis_increment=10, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + rm_empty_y=T, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + + +stability_high = stab() +#### "sorted" stability plots #### + +# do stability plots for later assembly +stability_high_reorder = bp_stability_hmap( + snpcount_high_df, + reorder_position = T, + p_title = NULL, + my_ylab = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4], + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + y_max_override=max(merged_df3$pos_count) +) +stability_centre_reorder = bp_stability_hmap( + snpcount_low_high_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "", + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4], + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + y_max_override=max(merged_df3$pos_count) + +) + +stability_low_reorder = bp_stability_hmap( + snpcount_low_low_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = NULL, + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4], + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + y_max_override=max(merged_df3$pos_count) + +) + +png(paste0(outdir,"gid_average_stability-sorted.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high_reorder, stability_centre_reorder,stability_low_reorder, ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"gid_average_stability-unsorted.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_low, stability_centre, stability_high, ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"gid_snp.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_p_snps_low, logo_p_snps_centre, logo_p_snps_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"gid_or.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_or_low, logo_or_centre, logo_or_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() diff --git a/scripts/plotting/plotting_thesis/thesis-figs-katg.R b/scripts/plotting/plotting_thesis/thesis-figs-katg.R new file mode 100644 index 0000000..1525f04 --- /dev/null +++ b/scripts/plotting/plotting_thesis/thesis-figs-katg.R @@ -0,0 +1,383 @@ +# +# small_df3 = merged_df3[merged_df3['position'] <= 670 & merged_df3['position'] >= 300,] +# +# small_df3 = merged_df3[merged_df3['position'] <= 451 & merged_df3['position'] >= 300,] +# +# +# mutable_df3 = cbind(small_df3) +# plot_df = cbind(small_df3) + +source("~/git/LSHTM_analysis/config/katg.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +outdir="/home/pub/Work/LSHTM/Thesis_Plots/" + +# add "pos_count" position count column +merged_df3=merged_df3 %>% dplyr::add_count(position) +merged_df3$pos_count=merged_df3$n +merged_df3$n=NULL + +# re-sort the dataframe according to position count +sorted_df = cbind(merged_df3) +sorted_df = sorted_df %>% arrange(pos_count) + +# determine some split points + +# further split up our sorted dataframe into low/mid/high for stability plotting +#snpcount_high_df = sorted_df[mid_snpcount:nrow(sorted_df),] +#snpcount_mid_df = sorted_df[low_snpcount:mid_snpcount,] +#snpcount_low_df = sorted_df[1:low_snpcount,] + +snpcount_high_df = sorted_df[sorted_df$pos_count>1,] +snpcount_low_df = sorted_df[sorted_df$pos_count<2,] + +low_half_snpcount = round(nrow(snpcount_low_df)/2) + +snpcount_low_low_df = snpcount_low_df[1:low_half_snpcount,] +snpcount_low_high_df = snpcount_low_df[low_half_snpcount:nrow(snpcount_low_df),] + +#### Put sorted-by-pos_count bp_stability_hmaps here #### +#...if they look better :-) + +length(merged_df3$position) + +# manually picked numbers due to irregular position distribution +# low = 170 +# centre=300 +low=round(max(merged_df3$position)/3) +centre=low*2 + +# Generate heat bar for whole DF +heatbar_df = generate_distance_colour_map(merged_df3) + +heat_bar_legend = get_legend( + ggplot(merged_df3, aes(factor(ligand_distance), y=0)) + + theme(legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = ligand_distance, + + ) + ) + + scale_colour_gradientn( + name="Ligand\nDistance (Å)", + colours = magma(5, direction = -1)) +) +# Position Tile Legend +# STR #00ff00 +tile_colour_map = c("INH" = "#00ff00" + ,"HEME" = "darkslategrey") + +tile_legend=get_legend( + + ggplot(tile_map, aes(factor(tile),y=0 + , colour=tile_colour + , fill=tile_colour))+ + geom_tile() + + theme(legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt")) + + scale_colour_manual(name=NULL + #, values = tile_map$tile_colour + , values=tile_colour_map) + + scale_fill_manual(name=NULL + #,values=tile_map$tile_colour + , values = tile_colour_map) +) + +# SAV Tile Red/Blue Legend + +temp_snp_df=merged_df3%>%dplyr::add_count(position) + +snp_tile_legend = get_legend( + ggplot(temp_snp_df, aes(factor(avg_stability_scaled), y=0)) + + theme(#legend.title = element_blank(), + legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = avg_stability_scaled, + + ) + ) + + scale_colour_gradientn( + name="Protein Stability\n(ΔΔG Kcal/mol)", colours = rev( + c("#007ba4", # rev() because FUCK + "#00bfc4", + "#ffffff", + "#ffd9d0", + "#f8766d", + "#ae3121") + ) + ) +) +# Logop Chemistry Colours +temp_logo_p_snp=cbind(merged_df3) +tab_mt = table(temp_logo_p_snp[["mutant_type"]], temp_logo_p_snp[["position"]]) +tab_mt = unclass(tab_mt) +logoplot_chemistry_legend=get_legend( + ggplot() + + geom_logo( + tab_mt + , method = 'custom' + , col_scheme = "chemistry" + , seq_type = 'aa') + theme( + legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt") + ) +) + +stab=function(...){bp_stability_hmap( + small_df3, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "", + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4], + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + y_max_override=max(merged_df3$pos_count), + ... +) +} + + + +# LOW +small_df3 = merged_df3[merged_df3['position'] < low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting LOW Rows: ", nrow(small_df3))) + + +logo_p_snps_low = LogoPlotSnps( + mutable_df3, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +# needs to be: active_aa_pos AND odds ratio >10 +logo_or_low = LogoPlotCustomH( + small_df3, + y_lab = "Odds Ratio", + # y_axis_increment=1, + # y_axis_log = TRUE, + rm_empty_y=TRUE, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_low = stab() + + +# MID +small_df3 = merged_df3[merged_df3['position'] <= centre & merged_df3['position'] >= low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting MID Rows: ", nrow(small_df3))) + +logo_p_snps_centre = LogoPlotSnps( + mutable_df3, + y_lab = "", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +logo_or_centre = LogoPlotCustomH( + small_df3, + y_axis_increment=2, + rm_empty_y=TRUE, + y_lab = "Log(10) Odds Ratio", + y_axis_log = TRUE, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_centre = stab() + +# HIGH +small_df3 = merged_df3[merged_df3['position'] > centre,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting HIGH Rows: ", nrow(small_df3))) + + +logo_p_snps_high = LogoPlotSnps( + mutable_df3, + y_lab="" , + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +logo_or_high = LogoPlotCustomH( + small_df3, + y_axis_increment=10, + rm_empty_y=TRUE, + # y_axis_log = TRUE, + y_lab = "Odds Ratio", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + + +stability_high = stab() +#### "sorted" stability plots #### + +# do stability plots for later assembly +stability_high_reorder = bp_stability_hmap( + snpcount_high_df, + reorder_position = T, + p_title = NULL, + my_ylab = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_xaxls = 4, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] + +) +stability_centre_reorder = bp_stability_hmap( + snpcount_low_high_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "", + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_low_reorder = bp_stability_hmap( + snpcount_low_low_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = NULL, + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] + +) + +png(paste0(outdir,"katg_average_stability-sorted.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high_reorder, stability_centre_reorder, stability_low_reorder,ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"katg_average_stability-unsorted.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high, stability_centre, stability_low,ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"katg_snp.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_p_snps_low, logo_p_snps_centre, logo_p_snps_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"katg_or.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_or_low, logo_or_centre, logo_or_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() diff --git a/scripts/plotting/plotting_thesis/thesis-figs-pnca.R b/scripts/plotting/plotting_thesis/thesis-figs-pnca.R new file mode 100644 index 0000000..5055785 --- /dev/null +++ b/scripts/plotting/plotting_thesis/thesis-figs-pnca.R @@ -0,0 +1,271 @@ +# +# small_df3 = merged_df3[merged_df3['position'] <= 670 & merged_df3['position'] >= 300,] +# +# small_df3 = merged_df3[merged_df3['position'] <= 451 & merged_df3['position'] >= 300,] +# +# +# mutable_df3 = cbind(small_df3) +# plot_df = cbind(small_df3) + +source("~/git/LSHTM_analysis/config/pnca.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +outdir="/home/pub/Work/LSHTM/Thesis_Plots/" + +# add "pos_count" position count column +merged_df3=merged_df3 %>% dplyr::add_count(position) +merged_df3$pos_count=merged_df3$n +merged_df3$n=NULL + +# re-sort the dataframe according to position count +sorted_df = cbind(merged_df3) +sorted_df = sorted_df %>% arrange(pos_count) + +# determine some split points + +# further split up our sorted dataframe into low/mid/high for stability plotting +# snpcount_high_df = sorted_df[mid_snpcount:nrow(sorted_df),] +# snpcount_mid_df = sorted_df[low_snpcount:mid_snpcount,] +# snpcount_low_df = sorted_df[1:low_snpcount,] + +snpcount_high_df = sorted_df[sorted_df$pos_count>1,] +snpcount_low_df = sorted_df[sorted_df$pos_count<2,] + +low_half_snpcount = round(nrow(snpcount_low_df)/2) + +snpcount_low_low_df = snpcount_low_df[1:low_half_snpcount,] +snpcount_low_high_df = snpcount_low_df[low_half_snpcount:nrow(snpcount_low_df),] + + +# do stability plots for later assembly +stability_high = bp_stability_hmap(snpcount_high_df, + reorder_position = T, + p_title = NULL, + my_ylab = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, +) +stability_centre = bp_stability_hmap(snpcount_low_high_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "SAV Count", + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, +) +stability_low = bp_stability_hmap(snpcount_low_low_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = NULL, + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 6, # x-axis label size + my_yaxls = 8, # y-axis label size + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, +) + +low=round(max(merged_df3$position)/3) +centre=low*2 + +# Generate heat bar for whole DF +heatbar_df = generate_distance_colour_map(merged_df3) + +heat_bar_legend = get_legend( + ggplot(merged_df3, aes(factor(ligand_distance), y=0)) + + theme(legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = ligand_distance, + + ) + ) + + scale_colour_gradientn(name="Ligand\nDistance (Å)", colours = magma(5, direction = -1)) +) +# Position Tile Legend + +tile_colour_map = c("PZA" = "green" + ,"DPA" = "darkslategrey" + ,"CDL" = "navyblue" + ,"Ca" = "purple") + +tile_legend=get_legend( + + ggplot(tile_map, aes(factor(tile),y=0 + , colour=tile_colour + , fill=tile_colour))+ + geom_tile() + + theme(legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt")) + + scale_colour_manual(name=NULL + #, values = tile_map$tile_colour + , values=tile_colour_map) + + scale_fill_manual(name=NULL + #,values=tile_map$tile_colour + , values = tile_colour_map) +) + +# SAV Tile Red/Blue Legend + +temp_snp_df=merged_df3%>%dplyr::add_count(position) + +snp_tile_legend = get_legend( + ggplot(temp_snp_df, aes(factor(avg_stability_scaled), y=0)) + + theme(#legend.title = element_blank(), + legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = avg_stability_scaled, + + ) + ) + + scale_colour_gradientn(name="Protein Stabilityu\n(ΔΔG Kcal/mol)", colours = rev(c("#007ba4", # rev() because FUCK + "#00bfc4", + "#ffffff", + "#ffd9d0", + "#f8766d", + "#ae3121"))) +) +# Logop Chemistry Colours +temp_logo_p_snp=cbind(merged_df3) +tab_mt = table(temp_logo_p_snp[["mutant_type"]], temp_logo_p_snp[["position"]]) +tab_mt = unclass(tab_mt) +logoplot_chemistry_legend=get_legend( + ggplot() + + geom_logo(tab_mt + , method = 'custom' + , col_scheme = "chemistry" + , seq_type = 'aa') + theme(legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt")) +) + +# LOW +small_df3 = merged_df3[merged_df3['position'] < low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting LOW Rows: ", nrow(small_df3))) + + +logo_p_snps_low = LogoPlotSnps(mutable_df3, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + y_lab="" +) + +# needs to be: active_aa_pos AND odds ratio >10 +logo_or_low = LogoPlotCustomH(small_df3, + y_axis_increment=10, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + rm_empty_y=T) + + +# MID +small_df3 = merged_df3[merged_df3['position'] <= centre & merged_df3['position'] >= low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting MID Rows: ", nrow(small_df3))) +logo_p_snps_centre = LogoPlotSnps(mutable_df3, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + y_lab = "SAV Count" +) + +logo_or_centre = LogoPlotCustomH(small_df3, + y_axis_increment=10, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + rm_empty_y=T +) + +# HIGH +small_df3 = merged_df3[merged_df3['position'] > centre,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting HIGH Rows: ", nrow(small_df3))) + + +logo_p_snps_high = LogoPlotSnps(mutable_df3, + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + y_lab="" +) + +logo_or_high = LogoPlotCustomH(small_df3, + y_axis_increment=10, + y_lab="", + aa_pos_drug=aa_pos_drug, + active_aa_pos=active_aa_pos, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + rm_empty_y=T +) + +png(paste0(outdir,"pnca_average_stability.png"), height = 8.25, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high, stability_centre, stability_low,ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"pnca_snp.png"), height = 8.25, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_p_snps_low, logo_p_snps_centre, logo_p_snps_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"pnca_or.png"), height = 8.25, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_or_low, logo_or_centre, logo_or_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() diff --git a/scripts/plotting/plotting_thesis/thesis-figs-rpob.R b/scripts/plotting/plotting_thesis/thesis-figs-rpob.R new file mode 100644 index 0000000..36202ea --- /dev/null +++ b/scripts/plotting/plotting_thesis/thesis-figs-rpob.R @@ -0,0 +1,386 @@ +# +# small_df3 = merged_df3[merged_df3['position'] <= 670 & merged_df3['position'] >= 300,] +# +# small_df3 = merged_df3[merged_df3['position'] <= 451 & merged_df3['position'] >= 300,] +# +# +# mutable_df3 = cbind(small_df3) +# plot_df = cbind(small_df3) + +source("~/git/LSHTM_analysis/config/rpob.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") +outdir="/home/pub/Work/LSHTM/Thesis_Plots/" + +# add "pos_count" position count column +merged_df3=merged_df3 %>% dplyr::add_count(position) +merged_df3$pos_count=merged_df3$n +merged_df3$n=NULL + +# re-sort the dataframe according to position count +sorted_df = cbind(merged_df3) +sorted_df = sorted_df %>% arrange(pos_count) + +# determine some split points + +# further split up our sorted dataframe into low/mid/high for stability plotting +#snpcount_high_df = sorted_df[mid_snpcount:nrow(sorted_df),] +#snpcount_mid_df = sorted_df[low_snpcount:mid_snpcount,] +#snpcount_low_df = sorted_df[1:low_snpcount,] + +snpcount_high_df = sorted_df[sorted_df$pos_count>1,] +snpcount_low_df = sorted_df[sorted_df$pos_count<2,] + +low_half_snpcount = round(nrow(snpcount_low_df)/2) + +snpcount_low_low_df = snpcount_low_df[1:low_half_snpcount,] +snpcount_low_high_df = snpcount_low_df[low_half_snpcount:nrow(snpcount_low_df),] + +#### Put sorted-by-pos_count bp_stability_hmaps here #### +#...if they look better :-) + +length(merged_df3$position) + +# manually picked numbers due to irregular position distribution +# low = 170 +# centre=300 +low=round(max(merged_df3$position)/3) +centre=low*2 + +# Generate heat bar for whole DF +heatbar_df = generate_distance_colour_map(merged_df3) + +heat_bar_legend = get_legend( + ggplot(merged_df3, aes(factor(ligand_distance), y=0)) + + theme(legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = ligand_distance, + + ) + ) + + scale_colour_gradientn( + name="Ligand\nDistance (Å)", + colours = magma(5, direction = -1)) +) +# Position Tile Legend +# STR #00ff00 +tile_colour_map = c("RFP" = "#00ff00") + +tile_legend=get_legend( + + ggplot(tile_map, aes(factor(tile),y=0 + , colour=tile_colour + , fill=tile_colour))+ + geom_tile() + + theme(legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt")) + + scale_colour_manual(name=NULL + #, values = tile_map$tile_colour + , values=tile_colour_map) + + scale_fill_manual(name=NULL + #,values=tile_map$tile_colour + , values = tile_colour_map) +) + +# SAV Tile Red/Blue Legend + +temp_snp_df=merged_df3%>%dplyr::add_count(position) + +snp_tile_legend = get_legend( + ggplot(temp_snp_df, aes(factor(avg_stability_scaled), y=0)) + + theme(#legend.title = element_blank(), + legend.direction="horizontal", + legend.title = element_text(size=10), + #legend.text = element_text(size=10), + #legend.key.size = unit(10, "pt") + )+ + geom_point(aes( + colour = avg_stability_scaled, + + ) + ) + + scale_colour_gradientn( + name="Protein Stability\n(ΔΔG Kcal/mol)", colours = rev( + c("#007ba4", # rev() because FUCK + "#00bfc4", + "#ffffff", + "#ffd9d0", + "#f8766d", + "#ae3121") + ) + ) +) +# Logop Chemistry Colours +temp_logo_p_snp=cbind(merged_df3) +tab_mt = table(temp_logo_p_snp[["mutant_type"]], temp_logo_p_snp[["position"]]) +tab_mt = unclass(tab_mt) +logoplot_chemistry_legend=get_legend( + ggplot() + + geom_logo( + tab_mt + , method = 'custom' + , col_scheme = "chemistry" + , seq_type = 'aa') + theme( + legend.title = element_blank(), + legend.direction="horizontal", + legend.text = element_text(size=10), + legend.key.size = unit(10, "pt") + ) +) + +stab=function(...){bp_stability_hmap( + small_df3, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = "", + #my_xaxls = 4, # x-axis label size + my_yaxls = 8, # y-axis label size + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4], + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + y_max_override=max(merged_df3$pos_count), + ... +) +} + + + +# LOW +small_df3 = merged_df3[merged_df3['position'] < low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting LOW Rows: ", nrow(small_df3))) + + +logo_p_snps_low = LogoPlotSnps( + mutable_df3, + y_lab="", + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +# needs to be: active_aa_pos AND odds ratio >10 +logo_or_low = LogoPlotCustomH( + small_df3, + y_lab = " ", + y_axis_increment=50, + #y_axis_log = TRUE, + rm_empty_y=TRUE, + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_high = stab(my_xaxls = 5) + + +# MID +small_df3 = merged_df3[merged_df3['position'] <= centre & merged_df3['position'] >= low,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting MID Rows: ", nrow(small_df3))) + +logo_p_snps_centre = LogoPlotSnps( + mutable_df3, + y_lab = "", + x_ats = 4.5, + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +logo_or_centre = LogoPlotCustomH( + small_df3, + x_ats = 5, + rm_empty_y=TRUE, + y_lab = "Odds Ratio", + y_axis_increment=250, + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + + #y_axis_log = TRUE, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_centre = stab(my_xaxls = 4) + +# HIGH +small_df3 = merged_df3[merged_df3['position'] > centre,] +mutable_df3 = cbind(small_df3) +print(paste0("Plotting HIGH Rows: ", nrow(small_df3))) + + +logo_p_snps_high = LogoPlotSnps( + mutable_df3, + y_lab=" " , + x_ats = 5, + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +logo_or_high = LogoPlotCustomH( + small_df3, + rm_empty_y=TRUE, + y_lab = "", + y_axis_increment=50, + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + #y_axis_log = TRUE, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + + +stability_low = stab(my_xaxls = 4) +#### "sorted" stability plots #### + +# do stability plots for later assembly +stability_high_reorder = bp_stability_hmap( + snpcount_high_df, + reorder_position = T, + p_title = NULL, + my_ylab = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_xaxls = 4, # x-axis label size + my_yaxls = 8, # y-axis label size + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] + +) +stability_centre_reorder = bp_stability_hmap( + snpcount_low_high_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = NULL, + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 5, # x-axis label size + my_yaxls = 8, # y-axis label size + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] +) + +stability_low_reorder = bp_stability_hmap( + snpcount_low_low_df, + reorder_position = T, + p_title = NULL, + yvar_colname = 'avg_stability_scaled', + stability_colname = "avg_stability_scaled", + stability_outcome_colname = "avg_stability_outcome", + my_ylab = NULL, + y_max_override = max(snpcount_high_df$pos_count), + my_xaxls = 5, # x-axis label size + my_yaxls = 8, # y-axis label size + active_aa_pos=active_aa_pos, + aa_pos_drug=aa_pos_drug, + aa_pos_lig1=aa_pos_lig1, + aa_pos_lig2=aa_pos_lig2, + aa_pos_lig3=aa_pos_lig3, + drug_colour = tile_map$tile_colour[1], + lig1_colour = tile_map$tile_colour[2], + lig2_colour = tile_map$tile_colour[3], + lig3_colour = tile_map$tile_colour[4] + +) + +png(paste0(outdir,"rpob_average_stability-sorted.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high_reorder, stability_centre_reorder, stability_low_reorder,ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"rpob_average_stability-unsorted.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, snp_tile_legend, nrow=1), + cowplot::plot_grid(stability_high, stability_centre, stability_low,ncol=1), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"rpob_snp.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_p_snps_low, logo_p_snps_centre, logo_p_snps_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off() + +png(paste0(outdir,"rpob_or.png"), height = 6, width=11.5, unit="in",res=300) +cowplot::plot_grid( + cowplot::plot_grid(heat_bar_legend, tile_legend, logoplot_chemistry_legend, nrow=1), + cowplot::plot_grid(logo_or_low, logo_or_centre, logo_or_high, ncol=1, align='hv'), + ncol=1, + rel_heights = c(1,12) +) +dev.off()