diff --git a/scripts/plotting/plotting_thesis/alr/DONOTUSEdm_om_plots_layout_alr.R b/scripts/plotting/plotting_thesis/alr/DONOTUSEdm_om_plots_layout_alr.R
new file mode 100644
index 0000000..cab13ae
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/alr/DONOTUSEdm_om_plots_layout_alr.R
@@ -0,0 +1,177 @@
+# # source dm_om_plots.R
+# source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/dm_om_plots.R")
+# 
+# ##### plots to combine ####
+# duetP
+# foldxP
+# deepddgP
+# dynamut2P
+# genomicsP
+# consurfP
+# proveanP
+# snap2P
+# mcsmligP
+# mcsmlig2P
+# mcsmppi2P
+# 
+# 
+# # Plot labels
+# tit1 = "Stability changes"
+# tit2 = "Genomic measure"
+# tit3 = "Distance to partners"
+# tit4 = "Evolutionary Conservation"
+# tit5 = "Affinity changes"
+# pt_size = 30
+# 
+# theme_georgia <- function(...) {
+#   theme_gray(base_family = "sans", ...) + 
+#     theme(plot.title = element_text(face = "bold"))
+# }
+# 
+# 
+# title_theme <- calc_element("plot.title", theme_georgia())
+# 
+# pt1 = ggdraw() + 
+#   draw_label(
+#     tit1,
+#     fontfamily = title_theme$family,
+#     fontface = title_theme$face,
+#     #size = title_theme$size
+#     size = pt_size
+#   )
+# 
+# pt2 = ggdraw() + 
+#   draw_label(
+#     tit2,
+#     fontfamily = title_theme$family,
+#     fontface = title_theme$face,
+#     size = pt_size
+#   )
+# 
+# pt3 = ggdraw() + 
+#   draw_label(
+#     tit3,
+#     fontfamily = title_theme$family,
+#     fontface = title_theme$face,
+#     size = pt_size
+#   )
+# 
+# pt4 = ggdraw() + 
+#   draw_label(
+#     tit4,
+#     fontfamily = title_theme$family,
+#     fontface = title_theme$face,
+#     size = pt_size
+#   )
+# 
+# 
+# pt5 = ggdraw() + 
+#   draw_label(
+#     tit5,
+#     fontfamily = title_theme$family,
+#     fontface = title_theme$face,
+#     size = pt_size
+#   )
+# 
+# #======================
+# # Output plot function
+# #======================
+# OutPlot_dm_om = function(x){
+#   
+#   # dist b/w plot title and plot
+#   relH_tp = c(0.08, 0.92)
+#   
+#   my_label_size = 25
+#   #----------------
+#   # Top panel
+#   #----------------
+#   top_panel = cowplot::plot_grid(
+#     cowplot::plot_grid(pt1, 
+#                        cowplot::plot_grid(duetP, foldxP, deepddgP, dynamut2P
+#                                           , nrow = 1
+#                                           , labels = c("A", "B", "C", "D")
+#                                           , label_size = my_label_size)
+#                        , ncol = 1
+#                        , rel_heights = relH_tp
+#     ),
+#     NULL,
+#     cowplot::plot_grid(pt2, 
+#                        cowplot::plot_grid(genomicsP
+#                                           , nrow = 1
+#                                           , labels = c("E")
+#                                           , label_size = my_label_size)
+#                        , ncol = 1
+#                        , rel_heights = relH_tp
+#     ),
+#     NULL,
+#     cowplot::plot_grid(pt3, 
+#                        cowplot::plot_grid( #distanceP
+#                          distanceP_lig
+#                          , distanceP_ppi2
+#                          , nrow = 1
+#                          , labels = c("F", "G")
+#                          , label_size = my_label_size)
+#                        , ncol = 1
+#                        , rel_heights = relH_tp
+#     ),
+#     nrow = 1,
+#     rel_widths = c(2/7, 0.1/7, 0.5/7, 0.1/7,  1/7)
+#   )
+#   
+#   #----------------
+#   # Bottom panel
+#   #----------------
+#   bottom_panel = cowplot::plot_grid(
+#     cowplot::plot_grid(pt4, 
+#                        cowplot::plot_grid(consurfP, proveanP, snap2P
+#                                           , nrow = 1
+#                                           , labels = c("H", "I", "J")
+#                                           , label_size = my_label_size)
+#                        , ncol = 1
+#                        , rel_heights =relH_tp
+#     ),NULL,
+#     cowplot::plot_grid(pt5, 
+#                        cowplot::plot_grid(mcsmligP
+#                                           , mcsmlig2P
+#                                           , mcsmppi2P
+#                                           , nrow = 1
+#                                           , labels = c("K", "L", "M")
+#                                           , label_size = my_label_size)
+#                        , ncol = 1
+#                        , rel_heights = relH_tp
+#     ),NULL,
+#     nrow = 1, 
+#     rel_widths = c(3/6,0.1/6,3/6, 0.1/6 )
+#   )
+#   
+#   #-------------------------------
+#   # combine: Top and Bottom panel
+#   #-------------------------------
+#   cowplot::plot_grid (top_panel, bottom_panel
+#                       , nrow =2
+#                       , rel_widths = c(1, 1)
+#                       , align = "hv")
+# }
+# 
+# #=====================
+# # OutPlot: svg and png
+# #======================
+# dm_om_combinedP =  paste0(outdir_images
+#                           ,tolower(gene)
+#                           ,"_dm_om_all.svg")
+# 
+# cat("DM OM plots with stats:", dm_om_combinedP)
+# svg(dm_om_combinedP, width = 32, height = 18)
+# 
+# OutPlot_dm_om()
+# dev.off()
+# 
+# 
+# dm_om_combinedP_png =  paste0(outdir_images
+#                               ,tolower(gene)
+#                               ,"_dm_om_all.png")
+# cat("DM OM plots with stats:", dm_om_combinedP_png)
+# png(dm_om_combinedP_png, width = 32, height = 18, units = "in", res = 300)
+# 
+# OutPlot_dm_om()
+# dev.off()
diff --git a/scripts/plotting/plotting_thesis/alr/alr_ORandSNP_results.R b/scripts/plotting/plotting_thesis/alr/alr_ORandSNP_results.R
new file mode 100644
index 0000000..5567926
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/alr/alr_ORandSNP_results.R
@@ -0,0 +1,198 @@
+#!/usr/bin/env Rscript
+#source("~/git/LSHTM_analysis/config/gid.R")
+#source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
+
+#=======
+# output
+#=======
+outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
+###################################################################
+# FIXME: ADD distance to NA when SP replies
+
+geneL_normal  = c("pnca")
+geneL_na      = c("gid", "rpob")
+geneL_ppi2    = c("alr", "embb", "katg", "rpob")
+
+# LigDist_colname # from globals used
+# ppi2Dist_colname #from globals used 
+# naDist_colname #from globals used
+
+common_cols  = c("mutationinformation"
+                 , "X5uhc_position"
+                 , "X5uhc_offset"
+                 , "position"
+                 , "dst_mode"
+                 , "mutation_info_labels"
+                 , "sensitivity", dist_columns )
+
+
+########################################
+categ_cols_to_factor = grep( "_outcome|_info", colnames(merged_df3) )
+fact_cols = colnames(merged_df3)[categ_cols_to_factor]
+
+if (any(lapply(merged_df3[, fact_cols], class) == "character")){
+  cat("\nChanging", length(categ_cols_to_factor), "cols to factor")
+  merged_df3[, fact_cols] <- lapply(merged_df3[, fact_cols], as.factor)
+  if (all(lapply(merged_df3[, fact_cols], class) == "factor")){
+    cat("\nSuccessful: cols changed to factor")
+  }
+}else{
+  cat("\nRequested cols aready factors")
+}
+
+cat("\ncols changed to factor are:\n", colnames(merged_df3)[categ_cols_to_factor] )
+
+####################################
+# merged_df3: NECESSARY pre-processing
+###################################
+#df3 = merged_df3
+plot_cols = c("mutationinformation", "mutation_info_labels", "position", "dst_mode"
+                   , all_cols)
+
+all_cols = c(common_cols
+             , all_stability_cols
+             , all_affinity_cols
+             , all_conserv_cols)
+
+
+# counting
+foo = merged_df3[, c("mutationinformation"
+                     , "wild_pos"
+                     , "position"
+                     , "sensitivity"
+                     , "avg_lig_affinity"
+                     , "avg_lig_affinity_scaled"
+                     , "avg_lig_affinity_outcome"
+                     , "ligand_distance"
+                     , "ligand_affinity_change"
+                     , "affinity_scaled"
+                     , "ligand_outcome"
+                     , "consurf_colour_rev"
+                     , "consurf_outcome")]
+
+table(foo$consurf_outcome)
+foo2 = foo[foo$ligand_distance<10,]
+
+table(foo2$ligand_outcome)
+
+#############################
+# wide plots SNP
+# DRUG
+length(aa_pos_drug); aa_pos_drug
+drug = foo[foo$position%in%aa_pos_drug,]
+drug$wild_pos
+#CA
+
+
+###############################################
+# OR plot
+
+bar = merged_df3[, c("mutationinformation"
+                     , "wild_pos"
+                     , "position"
+                     , "sensitivity"
+                     , affinity_dist_colnames
+                     , "or_mychisq"
+                     , "pval_fisher"
+                     #, "pval_chisq"
+                     , "neglog_pval_fisher"
+                     , "log10_or_mychisq")]
+
+# bar$p_adj_bonferroni = p.adjust(bar$pval_fisher, method = "bonferroni")
+# bar$signif_bon = bar$p_adj_bonferroni
+# bar = dplyr::mutate(bar
+#                     , 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'))
+
+bar$p_adj_fdr = p.adjust(bar$pval_fisher, method = "BH")
+bar$signif_fdr = bar$p_adj_fdr
+bar = dplyr::mutate(bar
+                    , 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'))
+
+# sort df
+bar = bar[order(bar$or_mychisq, decreasing = T), ]
+bar = bar[, c("mutationinformation"
+               , "wild_pos"
+               , "position"
+               , "sensitivity"
+               , affinity_dist_colnames
+               , "or_mychisq"
+               #, "pval_fisher"
+               #, "pval_chisq"
+               #, "neglog_pval_fisher"
+               #, "log10_or_mychisq"
+               #, "signif_bon"
+               , "p_adj_fdr"
+               , "signif_fdr")]
+
+table(bar$sensitivity)
+
+str(bar)
+sen = bar[bar$or_mychisq<1,]
+sen = na.omit(sen)
+
+res = bar[bar$or_mychisq>1,]
+res = na.omit(res)
+
+# comp 
+bar_or = bar[!is.na(bar$or_mychisq),]
+table(bar_or$sensitivity)
+      
+sen1 = bar_or[bar_or$or_mychisq<1,] # sen and res ~OR
+res1 = bar_or[bar_or$or_mychisq>1,] # sen and res ~OR
+
+# sanity check
+if (nrow(bar_or) == nrow(sen1) + nrow(res1) ){
+  cat("\nPASS: df with or successfully sourced"
+      , "\nCalculating % of muts with OR>1")
+}else{
+  stop("Abort: df with or numbers mimatch")
+}
+
+# percent for OR muts
+pc_orR = nrow(res1)/(nrow(sen1) + nrow(res1)); pc_orR
+
+cat("Number of R muts with OR>1:", nrow(res1)
+  , "\nPercentage of muts with OR>1 i.e resistant:"
+    , pc_orR *100 )
+
+# muts with highest OR
+head(bar_or$mutationinformation, 10)
+
+# sort df
+bar_or = bar_or[order(bar_or$or_mychisq
+                      , bar_or$ligand_distance
+                      #, bar_or$nca_dist
+                      , bar_or$interface_dist
+                      , decreasing = T), ]
+
+nrow(bar_or)
+bar_or$drug_site = ifelse(bar_or$position%in%aa_pos_drug, "drug", "no")
+table(bar_or$drug_site)
+
+bar_or$plp_site = ifelse(bar_or$position%in%aa_pos_plp, "plp", "no")
+table(bar_or$plp_site)
+top10_or = bar_or[1:10,]
+
+# are these active sites
+top10_or$position[top10_or$position%in%active_aa_pos]
+
+
+# clostest most sig
+bar_or_lig = bar_or[bar_or$ligand_distance<10,]
+bar_or_lig = bar_or_lig[order(bar_or_lig$ligand_distance, -bar_or_lig$or_mychisq), ]
+table(bar_or_lig$signif_fdr)
+
+
+bar_or_ppi = bar_or[bar_or$interface_dist<10,]
+bar_or_ppi = bar_or_ppi[order(bar_or_ppi$interface_dist, -bar_or_ppi$or_mychisq), ]
+table(bar_or_ppi$signif_fdr)
diff --git a/scripts/plotting/plotting_thesis/alr/alr_appendix_tables.R b/scripts/plotting/plotting_thesis/alr/alr_appendix_tables.R
new file mode 100644
index 0000000..948ab78
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/alr/alr_appendix_tables.R
@@ -0,0 +1,459 @@
+#!/usr/bin/env Rscript       
+#source("~/git/LSHTM_analysis/config/alr.R")
+#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/") 
+cat("\nOutput dir for stats:", outdir_stats)
+###################################################################
+geneL_normal  = c("pnca")
+#geneL_na      = c("gid", "rpob")
+geneL_na_v2   = c("gid")
+geneL_ppi2    = c("alr", "embb", "katg", "rpob")
+geneL_both    = c("rpob")
+
+
+if (tolower(gene)%in%geneL_na_v2) {
+  gene_colnames = c("mcsm_na_affinity", "mcsm_na_outcome")
+}
+
+if (tolower(gene)%in%geneL_ppi2) {
+  gene_colnames = c("mcsm_ppi2_affinity", "mcsm_ppi2_outcome")
+}
+
+#from plotting_globals()
+LigDist_colname
+ppi2Dist_colname 
+naDist_colname
+
+delta_symbol #delta_symbol = "\u0394"; delta_symbol
+angstroms_symbol
+
+cat("\nAffinity Distance colnames:", length(affinity_dist_colnames)
+    , "\nThese are:", affinity_dist_colnames)
+#===========
+# Data used
+#===========
+df3 = merged_df3
+
+cols_to_output = c("position"
+                   , "sensitivity"
+                   , "mutationinformation"
+                   , affinity_dist_colnames[1]
+                   , "ligand_affinity_change"
+                   , "ligand_outcome"
+                   , "mmcsm_lig"
+                   , "mmcsm_lig_outcome"
+                   , affinity_dist_colnames[2]
+                   # #, affinity_dist_colnames[3]
+                   # , "mcsm_na_affinity"
+                   # , "mcsm_na_outcome"
+                   # #, "mcsm_nca_affinity"
+                   # #, "mcsm_nca_outcome"
+                   , gene_colnames
+                   , "maf"
+                   , "or_mychisq"
+                   , "pval_fisher")
+
+cols_to_output
+df3_output = df3[, cols_to_output]
+colnames(df3_output)
+cat("\nSelecting columns:", length(colnames(df3_output)))
+#===============================================
+# Add COLS and rounding: adjusted P-values + MAF
+#==============================================
+#-----------------------------
+# adjusted P-values
+#-----------------------------
+# add cols: p_adj_fdr and signif_fdr
+df3_output$p_adj_fdr = p.adjust(df3_output$pval_fisher, method = "fdr")
+df3_output$signif_fdr = df3_output$p_adj_fdr
+df3_output = dplyr::mutate(df3_output
+                           , 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'))
+# rounding
+df3_output$or_mychisq = round(df3_output$or_mychisq,2)
+df3_output$p_adj_fdr  = round(df3_output$p_adj_fdr,2)
+head(df3_output)
+
+#----------
+# MAF (%)
+#----------
+# add col maf_percent
+df3_output$maf_percent = df3_output$maf*100
+
+# rounding
+df3_output$maf_percent = round(df3_output$maf_percent,2)
+head(df3_output$af); head(df3_output$maf);head(df3_output$maf_percent)
+
+#----------
+# P-value
+#----------
+df3_output$pval_fisher = round(df3_output$pval_fisher,2)
+
+class(df3_output)
+head(df3_output)
+
+####################################
+# Appendix: ligand affinity
+####################################
+df_lig = df3_output[df3_output[[LigDist_colname]]<DistCutOff,]
+
+cols_to_output_lig = c("position"
+                       , "sensitivity"
+                       , "mutationinformation"
+                       , LigDist_colname
+                       , "ligand_affinity_change"
+                       , "ligand_outcome"
+                       , "mmcsm_lig"
+                       , "mmcsm_lig_outcome"
+                       , "maf_percent"
+                       , "or_mychisq"
+                       , "pval_fisher"
+                       , "p_adj_fdr"
+                       , "signif_fdr")
+# select cols
+Out_df_lig = df_lig[, cols_to_output_lig]
+
+# sort df by OR and then MAF: highest OR and highest MAF
+#Out_df_ligS1 = Out_df_lig[order(Out_df_lig$or_mychisq, decreasing = T), ]
+Out_df_ligS = Out_df_lig[order(-Out_df_lig$or_mychisq, Out_df_lig$maf_percent), ]
+
+#head(Out_df_ligS1); tail(Out_df_ligS1)
+head(Out_df_ligS); tail(Out_df_ligS)
+
+colsNames_to_output_lig = c("position"
+                            , "sensitivity"
+                            , "Mutation"
+                            , paste0("Lig-Dist (", angstroms_symbol, ")")
+                            , "mCSM-ligand affinity"
+                            , "mCSM ligand_outcome"
+                            , "mmCSM-ligand affinity"
+                            , "mmCSM ligand_outcome"
+                            , paste0("MAF ","(%)")
+                            , "Odds Ratio"
+                            , "P-value"
+                            , "Adjusted P-value"
+                            , "Adjusted P-value significance")
+
+colnames(Out_df_ligS) = colsNames_to_output_lig
+head(Out_df_ligS)
+
+# ADD: active site annot
+nrow(Out_df_ligS)
+Out_df_ligS$drug_site = ifelse(Out_df_ligS$position%in%aa_pos_drug, "drug", "no")
+table(Out_df_ligS$drug_site)
+
+Out_df_ligS$plp_site = ifelse(Out_df_ligS$position%in%aa_pos_plp, "plp", "no")
+table(Out_df_ligS$plp_site)
+
+#--------------------
+# write output file: KS test within grpup
+#----------------------
+Out_ligT = paste0(outdir_stats
+                  , tolower(gene)
+                  , "_lig_muts.csv")
+
+cat("Output of Ligand muts:", Out_ligT )
+write.csv(Out_df_ligS, Out_ligT, row.names = FALSE)
+
+########################################################################
+####################################
+# Appendix: NA affinity or PPI2
+# naDist_colname
+####################################
+# Filtered data
+#df_nca = df3_output[df3_output[[naDist_colname]]<DistCutOff,]
+df_nca = df3_output[df3_output[[ppi2Dist_colname]]<DistCutOff,]
+
+# select cols
+cols_to_output_nca = c("position"
+                        , "sensitivity"
+                        , "mutationinformation"
+                        #, naDist_colname
+                        , ppi2Dist_colname
+                        , gene_colnames
+                        , "maf_percent"
+                        , "or_mychisq"
+                        , "pval_fisher"
+                        , "p_adj_fdr"
+                        , "signif_fdr")
+
+# extract output cols
+Out_df_nca = df_nca[, cols_to_output_nca]
+
+# sort df by OR and then MAF: Highest OR and Highest MAF
+#Out_df_ncaS = Out_df_nca[order(Out_df_nca$or_mychisq, decreasing = T), ]
+Out_df_ncaS = Out_df_nca[order(-Out_df_nca$or_mychisq, Out_df_nca$maf_percent), ]
+
+colsNames_to_output_nca = c("position"
+                             , "sensitivity"
+                             , "Mutation"
+                            
+                             #, paste0("NA-Dist (", angstroms_symbol, ")")
+                             #, paste0("mCSM-NA (", delta_symbol,delta_symbol,"G)")
+                             #, "mCSM-NA outcome"
+                             , paste0("PPI-Dist (", angstroms_symbol, ")")
+                             , paste0("mCSM-PPI (", delta_symbol,delta_symbol,"G)")
+                             , "mCSM-PPI outcome"
+                            
+                             , paste0("MAF ","(%)")
+                             , "Odds Ratio"
+                             , "P-value"
+                             , "Adjusted P-value"
+                             , "Adjusted P-value significance")
+
+colnames(Out_df_ncaS) = colsNames_to_output_nca
+Out_df_ncaS
+
+# ADD: active site annot
+nrow(Out_df_ncaS)
+Out_df_ncaS$drug_site = ifelse(Out_df_ncaS$position%in%aa_pos_drug, "drug", "no")
+table(Out_df_ncaS$drug_site)
+
+Out_df_ncaS$plp_site = ifelse(Out_df_ncaS$position%in%aa_pos_plp, "plp", "no")
+table(Out_df_ncaS$plp_site)
+
+#--------------------
+# write output file: KS test within grpup
+#----------------------
+Out_ncaT = paste0(outdir_stats
+                   , tolower(gene)
+                   , "_na_muts.csv")
+
+cat("Output of NA muts:", Out_ncaT )
+write.csv(Out_df_ncaS, Out_ncaT, row.names = FALSE)
+
+#################################################################################
+#################################################################################
+#################################################################################
+##########################################################
+# higest or/maf and stability effects
+###########################################################
+# convert to percet
+df3$maf_percent = df3$maf*100
+
+cols_to_output_effects = c("position"
+                   , "sensitivity"
+                   , "mutationinformation"
+                   , "avg_stability"
+                   , "avg_stability_outcome"
+                   , affinity_dist_colnames[1]
+                   , "avg_lig_affinity"
+                   , "avg_lig_affinity_outcome"
+                   , affinity_dist_colnames[2]
+                   , gene_colnames
+                   , "maf_percent"
+                   , "or_mychisq"
+                   , "pval_fisher")
+
+df3_effects = df3[, cols_to_output_effects]
+nrow(df3_effects); ncol(df3_effects)
+
+# add cols: p_adj_fdr and signif_fdr
+df3_effects$p_adj_fdr = p.adjust(df3_effects$pval_fisher, method = "fdr")
+df3_effects$signif_fdr = df3_effects$p_adj_fdr
+df3_effects = dplyr::mutate(df3_effects
+                           , 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'))
+# rounding
+df3_effects$or_mychisq = round(df3_effects$or_mychisq,2)
+df3_effects$p_adj_fdr  = round(df3_effects$p_adj_fdr,2)
+head(df3_effects)
+
+#-------------------
+# Highest OR and MAF
+#-------------------
+mut_hor  = df3_effects[df3_effects$or_mychisq%in%(max(df3_effects$or_mychisq, na.rm = T)),]
+mut_hmaf = df3_effects[df3_effects$maf_percent%in%(max(df3_effects$maf_percent, na.rm = T)),]
+
+if (identical(colnames(mut_hor), colnames(mut_hmaf)) ){
+  
+  # add cols
+  mut_hor$mutational_effect  = "Mutation with highest OR"
+  mut_hmaf$mutational_effect = "Most frequent mutation"
+  cat("\nPass: or and maf")
+}else{
+  quit("Abort: colnames or and maf mismatch")
+}
+
+#-------------------
+# Avg stability
+# most DD/SS: average stability
+#-------------------
+mut_h_avs_dd = df3_effects[df3_effects$avg_stability%in%(min(df3_effects$avg_stability, na.rm = T)),]
+mut_h_avs_ss = df3_effects[df3_effects$avg_stability%in%(max(df3_effects$avg_stability, na.rm = T)),]
+
+if (identical(colnames(mut_h_avs_dd), colnames(mut_h_avs_ss)) ){
+  
+  # add cols
+  mut_h_avs_dd$mutational_effect = "Most Destabilising for protomer"
+  mut_h_avs_ss$mutational_effect = "Most Stabilising for protomer"
+  
+  cat("\nPass : avg stability")
+}else{
+  quit("Abort: colnames stability mismatch")
+}
+
+#-------------------
+# Filtered columns
+# most DD/SS: ligand
+# FIXME DUBIOUS as min and max can be both negative
+#-------------------
+df3_effects_lig = df3_effects[df3_effects[[LigDist_colname]]<DistCutOff,]
+nrow(df3_effects_lig)
+
+mut_h_lig_dd = df3_effects_lig[df3_effects_lig$avg_lig_affinity%in%(min(df3_effects_lig$avg_lig_affinity, na.rm = T)),]
+mut_h_lig_ss = df3_effects_lig[df3_effects_lig$avg_lig_affinity%in%(max(df3_effects_lig$avg_lig_affinity, na.rm = T)),]
+
+if (identical(colnames(mut_h_lig_dd), colnames(mut_h_lig_ss)) ){
+  
+  # add cols
+  mut_h_lig_dd$mutational_effect = "Most Destabilising for Ligand affinity"
+  mut_h_lig_ss$mutational_effect = "CAUTION: Most DE/Stabilising for Ligand affinity"
+  
+  cat("\nPass: avg ligand affinity")
+}else{
+  quit("Abort: colnames lig mismatch")
+}
+
+#-------------------
+# Filtered columns
+# most DD/SS: NA
+#-------------------
+if (tolower(gene)%in%geneL_na_v2 ){
+  
+  df3_effects_na = df3_effects[df3_effects[[naDist_colname]]<DistCutOff,]
+  nrow(df3_effects_na)
+  mut_h_na_dd = df3_effects_na[df3_effects_na$mcsm_na_affinity%in%(min(df3_effects_na$mcsm_na_affinity, na.rm = T)),]
+  mut_h_na_ss = df3_effects_na[df3_effects_na$mcsm_na_affinity%in%(max(df3_effects_na$mcsm_na_affinity, na.rm = T)),]
+  
+  # add cols
+  mut_h_na_dd$mutational_effect = "Most Destabilising for NA affinity"
+  mut_h_na_ss$mutational_effect = "Most Stabilising for NA affinity"
+  
+  if (identical(colnames(mut_h_na_dd), colnames(mut_h_na_ss)) ){
+    cat("\nPass 1: NCA")
+  }else{
+    quit("Abort: colnames nca mismatch")
+  }
+  
+  if (identical(colnames(mut_h_na_dd), colnames(mut_h_lig_dd)) ){
+    cat("\nPass 2: NCA")
+  }else{
+    quit("Abort: colnames ppi2 mismatch")
+  }
+  #combine
+  gene_aff_combined = rbind(mut_h_na_dd, mut_h_na_ss)
+  
+}
+
+#-------------------
+# Filtered columns
+# most DD/SS: ppi2
+#-------------------
+if (tolower(gene)%in%geneL_ppi2 ){
+
+  df3_effects_ppi2 = df3_effects[df3_effects[[ppi2Dist_colname]]<DistCutOff,]
+  nrow(df3_effects_ppi2)
+
+  mut_h_ppi2_dd = df3_effects_ppi2[df3_effects_ppi2$mcsm_ppi2_affinity%in%(min(df3_effects_ppi2$mcsm_ppi2_affinity, na.rm = T)),]
+  mut_h_ppi2_ss = df3_effects_ppi2[df3_effects_ppi2$mcsm_ppi2_affinity%in%(max(df3_effects_ppi2$mcsm_ppi2_affinity, na.rm = T)),]
+  
+  # add cols
+  mut_h_ppi2_dd$mutational_effect = "Most Destabilising for PPI affinity"
+  mut_h_ppi2_ss$mutational_effect = "Most Stabilising for PPI affinity"
+  
+  if (identical(colnames(mut_h_ppi2_dd), colnames(mut_h_ppi2_ss)) ){
+    cat("\nPass 1: ppi2")
+  }else{
+    quit("Abort 1: colnames ppi2 mismatch")
+  }
+  
+  if (identical(colnames(mut_h_ppi2_dd), colnames(mut_h_lig_dd)) ){
+    cat("\nPass 2 : ppi2")
+  }else{
+    quit("Abort 2: colnames ppi2 mismatch")
+  }
+  #combine
+  gene_aff_combined = rbind(mut_h_ppi2_dd, mut_h_ppi2_ss)
+  
+}
+
+#============
+# final combine
+#============
+if ( identical(colnames(mut_hor), colnames(mut_h_lig_dd)) ){
+  cat("PASS: all")
+  combined_table = rbind(mut_hor, mut_hmaf,
+                         mut_h_avs_dd, mut_h_avs_ss,
+                         mut_h_lig_dd, mut_h_lig_ss,
+                         gene_aff_combined)
+  cat("\nCombined table dim:", "\nnrow:", nrow(combined_table), "\nncol:", ncol(combined_table))
+}else{
+  quit("Abort: colnames ppi2 mismatch")
+}
+
+# Assign pretty colnames
+colnames(combined_table)
+colsNames_combined_table = c("position"
+                            , "sensitivity"
+                            , "Mutation"
+                            , paste0("Avg stability (", delta_symbol,delta_symbol,"G)")
+                            , "avg stability outcome"
+                            
+                            , paste0("Lig-Dist (", angstroms_symbol, ")")
+                            , "Avg ligand affinity"
+                            , "Ligand affinity outcome"
+                            
+                            # , paste0("NA-Dist (", angstroms_symbol, ")")
+                            # , paste0("mCSM-NA (", delta_symbol,delta_symbol,"G)")
+                            # , "mCSM-NA outcome"
+                            # 
+                            
+                            , paste0("PPI-Dist (", angstroms_symbol, ")")
+                            , paste0("mCSM-PPI (", delta_symbol,delta_symbol,"G)")
+                            , "mCSM-PPI outcome"
+                            
+                            , paste0("MAF ","(%)")
+                            , "Odds Ratio"
+                            , "P-value"
+                            , "Adjusted P-value"
+                            , "Adjusted P-value significance"
+                            , "Mutational effect")
+
+if ( length(colnames(combined_table)) == length(colsNames_combined_table) ) {
+  cat("Assiging pretty colnames for output")
+  colnames(combined_table) <- colsNames_combined_table
+  #colnames(combined_table)
+}else{
+  stop("\nAbort: No. of cols mismatch. Cannot assign pretty colnames for output")
+}
+
+nrow(combined_table)
+combined_table$drug_site = ifelse(combined_table$position%in%aa_pos_drug, "drug", "no")
+table(combined_table$drug_site)
+
+combined_table$plp_site = ifelse(combined_table$position%in%aa_pos_plp, "plp", "no")
+table(combined_table$plp_site)
+
+#--------------------
+# write output file: KS test within grpup
+#----------------------
+Out_combined_effectsT = paste0(outdir_stats
+                  , tolower(gene)
+                  , "_mut_effects.csv")
+
+cat("Output of effects:", Out_combined_effectsT )
+write.csv(combined_table, Out_combined_effectsT, row.names = FALSE)
diff --git a/scripts/plotting/plotting_thesis/alr/alr_ks_test_all_PS.R b/scripts/plotting/plotting_thesis/alr/alr_ks_test_all_PS.R
new file mode 100755
index 0000000..0f5f961
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/alr/alr_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/alr.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)
diff --git a/scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist.R b/scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist.R
new file mode 100644
index 0000000..19a412f
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist.R
@@ -0,0 +1,165 @@
+#!/usr/bin/env Rscript  
+
+#########################################################
+# TASK: Lineage plots [merged_df2]
+# Count
+# Diversity
+# Average stability dist
+# Avergae affinity dist: optional
+#########################################################
+#=======
+# output
+#=======
+# outdir_images = paste0("~/git/Writing/thesis/images/results/"
+#                        , tolower(gene), "/")
+# cat("plots will output to:", outdir_images)
+#########################################################
+
+#===============
+#Quick numbers checks
+#===============
+nsample_lin = merged_df2[merged_df2$lineage%in%c("L1", "L2", "L3", "L4"),]
+
+if ( all(table(nsample_lin$sensitivity)== table(nsample_lin$mutation_info_labels)) ){
+  cat("\nTotal no. of samples belonging to L1-l4 for", gene,":", nrow(nsample_lin)
+      , "\nCounting R and S samples")
+  if( sum(table(nsample_lin$sensitivity)) ==  nrow(nsample_lin) ){
+    cat("\nPASSNumbers cross checked:")
+    print(table(nsample_lin$sensitivity))
+  }
+}else{
+  stop("Abort: Numbers mismatch. Please check")
+}
+########################################################################
+###################################################
+#                  Lineage barplots               #
+###################################################
+my_xats = 13 # x axis text size # were  25
+my_yats = 13# y axis text sized_lab_size
+my_xals = 13 # x axis label size
+my_yals = 13 # y axis label size
+my_lls  = 13 # legend label size
+d_lab_size = 5
+#===============================
+# lineage sample and SNP count
+#===============================
+lin_countP = lin_count_bp(lf_data = lineage_dfL[['lin_lf']]
+             , all_lineages = F
+             , x_categ = "sel_lineages"
+             , y_count = "p_count"
+             , use_lineages = c("L1", "L2", "L3", "L4")
+             , bar_fill_categ = "count_categ"
+             , display_label_col = "p_count"
+             , bar_stat_stype = "identity"
+             , d_lab_size = d_lab_size
+             , d_lab_col = "black"
+             , my_xats = my_xats # x axis text size
+             , my_yats = my_yats # y axis text sized_lab_size
+             , my_xals = my_xals # x axis label size
+             , my_yals = my_yals # y axis label size
+             , my_lls  = my_lls # legend label size
+             , bar_col_labels =  c("nsSNPs", "Total Samples")
+             , bar_col_values = c("grey50", "gray75")
+             , bar_leg_name = ""
+             , leg_location = "top"
+             , y_log10 = F
+             , y_scale_percent = FALSE
+             , y_label = c("Count")
+             )
+lin_countP
+#===============================
+# lineage SNP diversity count
+#===============================
+lin_diversityP = lin_count_bp_diversity(lf_data = lineage_dfL[['lin_wf']]
+                              , x_categ = "sel_lineages"
+                              , y_count = "snp_diversity"
+                              #, all_lineages = F
+                              , use_lineages = c("L1", "L2", "L3", "L4")
+                              , display_label_col = "snp_diversity_f"
+                              , bar_stat_stype = "identity"
+                              , x_lab_angle = 90
+                              , d_lab_size = d_lab_size
+                              , my_xats = my_xats # x axis text size
+                              , my_yats = my_yats # y axis text sized_lab_size
+                              , my_xals = my_xals # x axis label size
+                              , my_yals = my_yals # y axis label size
+                              , my_lls  = my_lls # legend label size
+                              , y_log10 = F
+                              , y_scale_percent = F
+                              , leg_location = "top"
+                              , y_label = "Percent" #"SNP diversity"
+                              , bp_plot_title = "nsSNP diversity"
+                              , title_colour = "black" #"chocolate4"
+                              , subtitle_text = NULL
+                              , sts = 10
+                              , subtitle_colour = "#350E20FF")
+lin_diversityP
+
+
+###################################################
+#                   Stability dist                #
+###################################################
+# scaled_cols_stability = c("duet_scaled"       
+#                           , "deepddg_scaled"   
+#                           , "ddg_dynamut2_scaled"
+#                           , "foldx_scaled"
+#                           , "avg_stability_scaled")
+
+my_ats = 13 # x axis text size # were  25
+my_als = 13# y axis text sized_lab_size
+my_leg_ts = 13 # x axis label size
+my_leg_title =11 # y axis label size
+my_strip_ts  = 13 #
+
+
+my_xlabel = paste0("Average stability ", "(", stability_suffix, ")"); my_xlabel
+#plotdf = merged_df2[merged_df2$lineage%in%c("L1", "L2", "L3", "L4"),]
+
+linP_dm_om = lineage_distP(merged_df2
+                           , with_facet = F
+                           , x_axis = "avg_stability_scaled"
+                           , y_axis = "lineage_labels"
+                           , x_lab = my_xlabel
+                           , use_lineages = c("L1", "L2", "L3", "L4")
+                           #, fill_categ = "mutation_info_orig", fill_categ_cols = c("#E69F00", "#999999")
+                           , fill_categ = "sensitivity"
+                           , fill_categ_cols = c("red", "blue")
+                           , label_categories = c("Resistant", "Sensitive")
+                           , leg_label = "Mutation group"
+                           , my_ats = my_ats # axis text size
+                           , my_als = my_als # axis label size
+                           , my_leg_ts = my_leg_ts
+                           , my_leg_title = my_leg_title
+                           , my_strip_ts = my_strip_ts
+                           , alpha = 0.56
+)
+
+linP_dm_om
+
+###################################################
+#                 Affinity dist [OPTIONAL]        #
+###################################################
+# scaled_cols_affinity = c("affinity_scaled" 
+#                          , "mmcsm_lig_scaled" 
+#                          , "mcsm_ppi2_scaled" 
+#                          , "mcsm_na_scaled"
+#                          , "avg_lig_affinity_scaled")
+
+# lineage_distP(merged_df2
+#               , with_facet = F
+#               , x_axis = "avg_lig_affinity_scaled"
+#               , y_axis = "lineage_labels"
+#               , x_lab = my_xlabel
+#               , use_lineages = c("L1", "L2", "L3", "L4")
+#               #, fill_categ = "mutation_info_orig", fill_categ_cols = c("#E69F00", "#999999")
+#               , fill_categ = "sensitivity"
+#               , fill_categ_cols = c("red", "blue")
+#               , label_categories = c("Resistant", "Sensitive")
+#               , leg_label = "Mutation group"
+#               , my_ats = 22 # axis text size
+#               , my_als = 22 # axis label size
+#               , my_leg_ts = 22
+#               , my_leg_title = 22
+#               , my_strip_ts = 22
+#               , alpha = 0.56
+# )
diff --git a/scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist_layout.R b/scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist_layout.R
new file mode 100644
index 0000000..7144093
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist_layout.R
@@ -0,0 +1,75 @@
+#!/usr/bin/env Rscript  
+source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/alr/alr_lineage_bp_dist.R")
+
+###########################################
+# TASK: generate plots for lineage
+# Individual plots in 
+#lineage_bp_both.R
+#linage_dist_ens_stability.R
+###########################################
+#### just A+B ####
+# png
+my_label_size = 12
+linPlots_combined = paste0(outdir_images
+                           , tolower(gene)
+                           ,"_linP_AB_combined.png")
+
+
+cat("\nOutput plot:", linPlots_combined)
+png(linPlots_combined, width = 9, height = 6, units = "in", res = 300)
+
+
+cowplot::plot_grid(lin_countP, lin_diversityP
+                     , ncol = 2
+                     , labels = "AUTO"
+                     , label_size = my_label_size)
+dev.off()
+
+# svg
+# linPlots_combined = paste0(outdir_images
+#                         , tolower(gene)
+#                         ,"_linP_combined.svg")
+# 
+# cat("\nOutput plot:", linPlots_combined)
+# svg(linPlots_combined, width = 18, height = 12)
+# 
+# cowplot::plot_grid(
+#   cowplot::plot_grid(lin_countP, lin_diversityP
+#                      , nrow = 2
+#                      , rel_heights = c(1.2,1)
+#                      , labels = "AUTO"
+#                      , label_size = my_label_size),
+#   NULL,
+#   linP_dm_om,
+#   nrow = 1,
+#   labels = c("", "", "C"),
+#   label_size = my_label_size,
+#   rel_widths = c(35, 3, 52)
+# )
+# dev.off()
+
+# png
+# my_label_size = 12
+# linPlots_combined = paste0(outdir_images
+#                            , tolower(gene)
+#                            ,"_linP_combined.png")
+# 
+# cat("\nOutput plot:", linPlots_combined)
+# png(linPlots_combined, width = 9, height = 6, units = "in" ,res = 300)
+# 
+# cowplot::plot_grid(
+#   cowplot::plot_grid(lin_countP, lin_diversityP
+#                      , nrow = 2
+#                      , rel_heights = c(1.2,1)
+#                      , labels = "AUTO"
+#                      , label_size = my_label_size),
+#   NULL,
+#   linP_dm_om,
+#   nrow = 1,
+#   labels = c("", "", "C"),
+#   label_size = my_label_size,
+#   rel_widths = c(35, 3, 52)
+# )
+# dev.off()
+
+
diff --git a/scripts/plotting/plotting_thesis/lineage_bp_dist.R b/scripts/plotting/plotting_thesis/lineage_bp_dist.R
new file mode 100644
index 0000000..d863e8f
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/lineage_bp_dist.R
@@ -0,0 +1,165 @@
+#!/usr/bin/env Rscript  
+
+#########################################################
+# TASK: Lineage plots [merged_df2]
+# Count
+# Diversity
+# Average stability dist
+# Avergae affinity dist: optional
+#########################################################
+#=======
+# output
+#=======
+# outdir_images = paste0("~/git/Writing/thesis/images/results/"
+#                        , tolower(gene), "/")
+# cat("plots will output to:", outdir_images)
+#########################################################
+
+#===============
+#Quick numbers checks
+#===============
+nsample_lin = merged_df2[merged_df2$lineage%in%c("L1", "L2", "L3", "L4"),]
+
+if ( all(table(nsample_lin$sensitivity)== table(nsample_lin$mutation_info_labels)) ){
+  cat("\nTotal no. of samples belonging to L1-l4 for", gene,":", nrow(nsample_lin)
+      , "\nCounting R and S samples")
+  if( sum(table(nsample_lin$sensitivity)) ==  nrow(nsample_lin) ){
+    cat("\nPASSNumbers cross checked:")
+    print(table(nsample_lin$sensitivity))
+  }
+}else{
+  stop("Abort: Numbers mismatch. Please check")
+}
+########################################################################
+###################################################
+#                  Lineage barplots               #
+###################################################
+my_xats = 8 # x axis text size # were  25
+my_yats = 8# y axis text sized_lab_size
+my_xals = 8 # x axis label size
+my_yals = 8 # y axis label size
+my_lls  = 8 # legend label size
+d_lab_size = 2.3
+#===============================
+# lineage sample and SNP count
+#===============================
+lin_countP = lin_count_bp(lf_data = lineage_dfL[['lin_lf']]
+             , all_lineages = F
+             , x_categ = "sel_lineages"
+             , y_count = "p_count"
+             , use_lineages = c("L1", "L2", "L3", "L4")
+             , bar_fill_categ = "count_categ"
+             , display_label_col = "p_count"
+             , bar_stat_stype = "identity"
+             , d_lab_size = d_lab_size
+             , d_lab_col = "black"
+             , my_xats = my_xats # x axis text size
+             , my_yats = my_yats # y axis text sized_lab_size
+             , my_xals = my_xals # x axis label size
+             , my_yals = my_yals # y axis label size
+             , my_lls  = my_lls # legend label size
+             , bar_col_labels =  c("nsSNPs", "Total Samples")
+             , bar_col_values = c("grey50", "gray75")
+             , bar_leg_name = ""
+             , leg_location = "top"
+             , y_log10 = F
+             , y_scale_percent = FALSE
+             , y_label = c("Count")
+             )
+lin_countP
+#===============================
+# lineage SNP diversity count
+#===============================
+lin_diversityP = lin_count_bp_diversity(lf_data = lineage_dfL[['lin_wf']]
+                              , x_categ = "sel_lineages"
+                              , y_count = "snp_diversity"
+                              #, all_lineages = F
+                              , use_lineages = c("L1", "L2", "L3", "L4")
+                              , display_label_col = "snp_diversity_f"
+                              , bar_stat_stype = "identity"
+                              , x_lab_angle = 90
+                              , d_lab_size = d_lab_size
+                              , my_xats = my_xats # x axis text size
+                              , my_yats = my_yats # y axis text sized_lab_size
+                              , my_xals = my_xals # x axis label size
+                              , my_yals = my_yals # y axis label size
+                              , my_lls  = my_lls # legend label size
+                              , y_log10 = F
+                              , y_scale_percent = F
+                              , leg_location = "top"
+                              , y_label = "Percent" #"SNP diversity"
+                              , bp_plot_title = "nsSNP diversity"
+                              , title_colour = "black" #"chocolate4"
+                              , subtitle_text = NULL
+                              , sts = 10
+                              , subtitle_colour = "#350E20FF")
+lin_diversityP
+
+
+###################################################
+#                   Stability dist                #
+###################################################
+# scaled_cols_stability = c("duet_scaled"       
+#                           , "deepddg_scaled"   
+#                           , "ddg_dynamut2_scaled"
+#                           , "foldx_scaled"
+#                           , "avg_stability_scaled")
+
+my_ats = 8 # x axis text size # were  25
+my_als = 8# y axis text sized_lab_size
+my_leg_ts = 8 # x axis label size
+my_leg_title = 8 # y axis label size
+my_strip_ts  = 8 #
+
+
+my_xlabel = paste0("Average stability ", "(", stability_suffix, ")"); my_xlabel
+#plotdf = merged_df2[merged_df2$lineage%in%c("L1", "L2", "L3", "L4"),]
+
+linP_dm_om = lineage_distP(merged_df2
+                           , with_facet = F
+                           , x_axis = "avg_stability_scaled"
+                           , y_axis = "lineage_labels"
+                           , x_lab = my_xlabel
+                           , use_lineages = c("L1", "L2", "L3", "L4")
+                           #, fill_categ = "mutation_info_orig", fill_categ_cols = c("#E69F00", "#999999")
+                           , fill_categ = "sensitivity"
+                           , fill_categ_cols = c("red", "blue")
+                           , label_categories = c("Resistant", "Sensitive")
+                           , leg_label = "Mutation group"
+                           , my_ats = my_ats # axis text size
+                           , my_als = my_als # axis label size
+                           , my_leg_ts = my_leg_ts
+                           , my_leg_title = my_leg_title
+                           , my_strip_ts = my_strip_ts
+                           , alpha = 0.56
+)
+
+linP_dm_om
+
+###################################################
+#                 Affinity dist [OPTIONAL]        #
+###################################################
+# scaled_cols_affinity = c("affinity_scaled" 
+#                          , "mmcsm_lig_scaled" 
+#                          , "mcsm_ppi2_scaled" 
+#                          , "mcsm_na_scaled"
+#                          , "avg_lig_affinity_scaled")
+
+# lineage_distP(merged_df2
+#               , with_facet = F
+#               , x_axis = "avg_lig_affinity_scaled"
+#               , y_axis = "lineage_labels"
+#               , x_lab = my_xlabel
+#               , use_lineages = c("L1", "L2", "L3", "L4")
+#               #, fill_categ = "mutation_info_orig", fill_categ_cols = c("#E69F00", "#999999")
+#               , fill_categ = "sensitivity"
+#               , fill_categ_cols = c("red", "blue")
+#               , label_categories = c("Resistant", "Sensitive")
+#               , leg_label = "Mutation group"
+#               , my_ats = 22 # axis text size
+#               , my_als = 22 # axis label size
+#               , my_leg_ts = 22
+#               , my_leg_title = 22
+#               , my_strip_ts = 22
+#               , alpha = 0.56
+# )
diff --git a/scripts/plotting/plotting_thesis/lineage_bp_dist_layout.R b/scripts/plotting/plotting_thesis/lineage_bp_dist_layout.R
new file mode 100644
index 0000000..1f229df
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/lineage_bp_dist_layout.R
@@ -0,0 +1,56 @@
+#!/usr/bin/env Rscript  
+source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/linage_bp_dist.R")
+###########################################
+# TASK: generate plots for lineage
+# Individual plots in 
+#lineage_bp_both.R
+#linage_dist_ens_stability.R
+###########################################
+# svg
+# linPlots_combined = paste0(outdir_images
+#                         , tolower(gene)
+#                         ,"_linP_combined.svg")
+# 
+# cat("\nOutput plot:", linPlots_combined)
+# svg(linPlots_combined, width = 18, height = 12)
+# 
+# cowplot::plot_grid(
+#   cowplot::plot_grid(lin_countP, lin_diversityP
+#                      , nrow = 2
+#                      , rel_heights = c(1.2,1)
+#                      , labels = "AUTO"
+#                      , label_size = my_label_size),
+#   NULL,
+#   linP_dm_om,
+#   nrow = 1,
+#   labels = c("", "", "C"),
+#   label_size = my_label_size,
+#   rel_widths = c(35, 3, 52)
+# )
+# dev.off()
+
+# png
+my_label_size = 12
+linPlots_combined = paste0(outdir_images
+                            , tolower(gene)
+                            ,"_linP_combined.png")
+
+cat("\nOutput plot:", linPlots_combined)
+png(linPlots_combined, width = 9, height = 6, units = "in" ,res = 300)
+ 
+ cowplot::plot_grid(
+   cowplot::plot_grid(lin_countP, lin_diversityP
+#                      , nrow = 2
+#                      , rel_heights = c(1.2,1)
+#                      , labels = "AUTO"
+#                      , label_size = my_label_size),
+#   NULL,
+#   linP_dm_om,
+#   nrow = 1,
+#   labels = c("", "", "C"),
+#   label_size = my_label_size,
+#   rel_widths = c(35, 3, 52)
+# )
+# dev.off()
+
+