diff --git a/scripts/functions/lineage_dist.R b/scripts/functions/lineage_dist.R
index 5f5db56..4e3689a 100644
--- a/scripts/functions/lineage_dist.R
+++ b/scripts/functions/lineage_dist.R
@@ -19,6 +19,7 @@ lineage_distP <- function(plotdf
                           , facet_wrap_var = "" # FIXME: document what this is for
                           , fill_categ = "mutation_info_labels"
                           , fill_categ_cols = c("#E69F00", "#999999")
+                          , label_categories = c("LABEL1", "LABEL2")
                           , my_ats = 15 # axis text size
                           , my_als = 20 # axis label size
                           , my_leg_ts = 16
@@ -27,7 +28,8 @@ lineage_distP <- function(plotdf
                           , leg_pos = c(0.8, 0.9)
                           , leg_pos_wf = c("top", "left", "bottom", "right")
                           , leg_dir_wf = c("horizontal", "vertical")
-                          , leg_label = "")
+                          , leg_label = ""
+                          , alpha = 0.7)
 
 {
   
@@ -41,10 +43,11 @@ lineage_distP <- function(plotdf
     geom_density_ridges(aes_string(fill = fill_categ)
                         , scale = 3
                         , size = 0.3
-                        , alpha = 0.8) +
+                        , alpha = alpha) +
     scale_x_continuous(expand = c(0.01, 0.01)) +
     #coord_cartesian( xlim = c(-1, 1)) +
-    scale_fill_manual(values = fill_categ_cols) +
+    scale_fill_manual(values = fill_categ_cols
+                      , labels = label_categories) +
     theme(axis.text.x = element_text(size = my_ats
                                      , angle = 90
                                      , hjust = 1
diff --git a/scripts/functions/stability_count_bp.R b/scripts/functions/stability_count_bp.R
index 190a71b..a40e1bf 100644
--- a/scripts/functions/stability_count_bp.R
+++ b/scripts/functions/stability_count_bp.R
@@ -62,10 +62,13 @@ stability_count_bp <- function(plotdf
          , subtitle = subtitle_text
          , y        = yaxis_title) + 
     
-    scale_fill_discrete(name = leg_title
-                        , labels = label_categories) +
+    # scale_fill_discrete(name = leg_title
+    #                     , labels = label_categories) +
 
-    scale_fill_manual("", values=bar_fill_values)
+    scale_fill_manual(name = ""
+                      # name = leg_title
+                      , values = bar_fill_values
+                      , labels = label_categories)
   
   
   return(OutPlot_count)
diff --git a/scripts/plotting/structure_figures/mcsm_mean_stability.R b/scripts/plotting/structure_figures/mcsm_mean_stability.R
deleted file mode 100755
index 7e35e79..0000000
--- a/scripts/plotting/structure_figures/mcsm_mean_stability.R
+++ /dev/null
@@ -1,163 +0,0 @@
-getwd()
-setwd("~/git/LSHTM_analysis/scripts/plotting")
-getwd()
-
-#########################################################
-# TASK:
-
-#########################################################
-#source("~/git/LSHTM_analysis/scripts/Header_TT.R")
-#require(data.table)
-#require(dplyr)
-
-source("plotting_data.R")
-# should return
-#my_df
-#my_df_u
-#dup_muts
-
-# cmd parse arguments
-#require('getopt', quietly = TRUE)
-#========================================================
-
-
-#========================================================
-#	Read file: call script for combining df for PS
-
-#source("../combining_two_df.R")
-
-#========================================================
-
-# plotting_data.R imports all the dir names, etc
-
-#=======
-# output
-#=======
-out_filename_mean_stability = paste0(tolower(gene), "_mean_stability.csv") 
-outfile_mean_stability = paste0(outdir, "/", out_filename_mean_stability)
-print(paste0("Output file:", outfile_mean_stability))
-
-#%%===============================================================
-
-#================
-# Data for plots
-#================
-# REASSIGNMENT as necessary
-df  = my_df_u
-rm(my_df)
-
-###########################
-# Data for bfactor figure
-# PS (duet) average 
-# Ligand affinity average
-###########################
-head(df$position); head(df$mutationinformation)
-head(df$duet_stability_change)
-
-# order data frame 
-#df = df[order(df$position),] #already done
-#head(df$position); head(df$mutationinformation)
-#head(df$duet_stability_change)
-
-#***********
-# PS(duet): average by position and then scale b/w -1 and 1
-# column to average: duet_stability_change (NOT scaled!)
-#***********
-mean_duet_by_position <- df %>%
-  group_by(position) %>%
-  summarize(averaged_duet = mean(duet_stability_change))
-
-# scale b/w -1 and 1
-duet_min = min(mean_duet_by_position['averaged_duet'])
-duet_max = max(mean_duet_by_position['averaged_duet']) 
-
-# scale the averaged_duet values
-mean_duet_by_position['averaged_duet_scaled'] = lapply(mean_duet_by_position['averaged_duet']
-                                                       , function(x) ifelse(x < 0, x/abs(duet_min), x/duet_max))
-
-cat(paste0('Average duet scores:\n', head(mean_duet_by_position['averaged_duet'])
-           , '\n---------------------------------------------------------------'
-           , '\nScaled duet scores:\n', head(mean_duet_by_position['averaged_duet_scaled'])))
-
-# sanity checks
-l_bound_duet = min(mean_duet_by_position['averaged_duet_scaled'])
-u_bound_duet = max(mean_duet_by_position['averaged_duet_scaled'])
-
-if ( (l_bound_duet == -1) && (u_bound_duet == 1) ){
-  cat(paste0("PASS: duet scores averaged by position and then scaled"
-        , "\nmin averaged duet: ", l_bound_duet
-        , "\nmax averaged duet: ", u_bound_duet))
-}else{
-  cat(paste0("FAIL: avergaed duet scores could not be scaled b/w -1 and 1"
-        , "\nmin averaged duet: ", l_bound_duet
-        , "\nmax averaged duet: ", u_bound_duet))
-  quit()
-} 
-
-#***********
-# Lig: average by position and then scale b/w -1 and 1
-# column: ligand_affinity_change (NOT scaled!)
-#***********
-mean_affinity_by_position <- df %>%
-  group_by(position) %>%
-  summarize(averaged_affinity = mean(ligand_affinity_change))
-
-# scale b/w -1 and 1
-affinity_min = min(mean_affinity_by_position['averaged_affinity'])
-affinity_max = max(mean_affinity_by_position['averaged_affinity']) 
-
-# scale the averaged_affinity values
-mean_affinity_by_position['averaged_affinity_scaled'] = lapply(mean_affinity_by_position['averaged_affinity']
-                                                               , function(x) ifelse(x < 0, x/abs(affinity_min), x/affinity_max))
-
-cat(paste0('Average affinity scores:\n', head(mean_affinity_by_position['averaged_affinity'])
-           , '\n---------------------------------------------------------------'
-           , '\nScaled affinity scores:\n', head(mean_affinity_by_position['averaged_affinity_scaled'])))
-
-# sanity checks
-l_bound_affinity = min(mean_affinity_by_position['averaged_affinity_scaled'])
-u_bound_affinity = max(mean_affinity_by_position['averaged_affinity_scaled'])
-
-if ( (l_bound_affinity == -1) && (u_bound_affinity == 1) ){
-  cat(paste0("PASS: affinity scores averaged by position and then scaled"
-             , "\nmin averaged affintiy: ", l_bound_affinity
-             , "\nmax averaged affintiy: ", u_bound_affinity))
-}else{
-  cat(paste0("FAIL: avergaed affinity scores could not be scaled b/w -1 and 1"
-             , "\nmin averaged affintiy: ", l_bound_affinity
-             , "\nmax averaged affintiy: ", u_bound_affinity))
-  quit()
-} 
-
-#***********
-# merge: mean_duet_by_position and mean_affinity_by_position
-#***********
-common_cols = intersect(colnames(mean_duet_by_position), colnames(mean_affinity_by_position))
-
-if (dim(mean_duet_by_position) && dim(mean_affinity_by_position)){
-  print(paste0("PASS: dim's match, mering dfs by column :", common_cols))
-  #combined = as.data.frame(cbind(mean_duet_by_position, mean_affinity_by_position ))
-  combined_df = as.data.frame(merge(mean_duet_by_position
-                                    , mean_affinity_by_position
-                                    , by = common_cols
-                                    , all = T))
-  
-  cat(paste0("\nnrows combined_df:", nrow(combined_df)
-               , "\nnrows combined_df:", ncol(combined_df)))
-}else{
-    cat(paste0("FAIL: dim's mismatch, aborting cbind!"
-          , "\nnrows df1:", nrow(mean_duet_by_position)
-          , "\nnrows df2:", nrow(mean_affinity_by_position)))
-    quit()      
-}
-#%%============================================================
-# output
-write.csv(combined_df, outfile_mean_stability
-          , row.names = F)
-cat("Finished writing file:\n"
-    , outfile_mean_stability
-    , "\nNo. of rows:", nrow(combined_df)
-    , "\nNo. of cols:", ncol(combined_df))
-
-# end of script
-#===============================================================
diff --git a/scripts/plotting/structure_figures/replaceBfactor_pdb.R b/scripts/plotting/structure_figures/replaceBfactor_pdb.R
deleted file mode 100755
index a840c19..0000000
--- a/scripts/plotting/structure_figures/replaceBfactor_pdb.R
+++ /dev/null
@@ -1,332 +0,0 @@
-#!/usr/bin/env Rscript                                                  
-
-#########################################################
-# TASK: Replace B-factors in the pdb file with the mean
-# normalised stability values.
-
-# read pdb file
-# make two copies so you can replace B factors for 1)duet
-# 2)affinity values and output 2 separate pdbs for
-# rendering on chimera
-
-# read mcsm mean stability value files
-# extract the respective mean values and assign to the
-# b-factor column within their respective pdbs
-
-# generate some distribution plots for inspection
-
-#########################################################
-# working dir and loading libraries
-getwd()
-setwd("~/git/LSHTM_analysis/scripts/plotting")
-cat(c(getwd(),"\n"))
-
-#source("~/git/LSHTM_analysis/scripts/Header_TT.R")
-library(bio3d)
-require("getopt", quietly = TRUE) # cmd parse arguments
-#========================================================
-#drug = "pyrazinamide"
-#gene = "pncA"
-
-# command line args
-spec = matrix(c(
-  "drug"   , "d", 1, "character",
-  "gene"   , "g", 1, "character"
-), byrow = TRUE, ncol = 4)
-
-opt = getopt(spec)
-
-drug = opt$drug
-gene = opt$gene
-
-if(is.null(drug)|is.null(gene)) {
-  stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)")
-}
-#========================================================
-gene_match = paste0(gene,"_p.")
-cat(gene_match)
-
-#=============
-# directories
-#=============
-datadir = paste0("~/git/Data")
-indir = paste0(datadir, "/", drug, "/input")
-outdir = paste0("~/git/Data", "/", drug, "/output")
-#outdir_plots = paste0("~/git/Data", "/", drug, "/output/plots")
-outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
-
-#======
-# input
-#======
-in_filename_pdb = paste0(tolower(gene), "_complex.pdb") 
-infile_pdb = paste0(indir, "/", in_filename_pdb)
-cat(paste0("Input file:", infile_pdb) )
-
-#in_filename_mean_stability = paste0(tolower(gene), "_mean_stability.csv")
-#infile_mean_stability = paste0(outdir, "/", in_filename_mean_stability)
-
-in_filename_mean_stability = paste0(tolower(gene), "_mean_ens_stab_aff.csv")
-infile_mean_stability = paste0(outdir_plots, "/", in_filename_mean_stability)
-
-cat(paste0("Input file:", infile_mean_stability) )
-
-#=======
-# output
-#=======
-#out_filename_duet_mspdb = paste0(tolower(gene), "_complex_bduet_ms.pdb")
-out_filename_duet_mspdb = paste0(tolower(gene), "_complex_b_stab_ms.pdb") 
-outfile_duet_mspdb = paste0(outdir_plots, "/", out_filename_duet_mspdb)
-print(paste0("Output file:", outfile_duet_mspdb))
-
-out_filename_lig_mspdb  = paste0(tolower(gene), "_complex_blig_ms.pdb") 
-outfile_lig_mspdb = paste0(outdir_plots, "/", out_filename_lig_mspdb)
-print(paste0("Output file:", outfile_lig_mspdb))
-
-#%%===============================================================
-#NOTE: duet here refers to the ensemble stability values
-
-###########################
-# Read file: average stability values
-# or mcsm_normalised file
-###########################
-my_df <- read.csv(infile_mean_stability, header = T)
-str(my_df)
-
-#############
-# Read pdb
-#############
-# list of 8
-my_pdb = read.pdb(infile_pdb
-                  , maxlines = -1
-                  , multi = FALSE 
-                  , rm.insert = FALSE
-                  , rm.alt = TRUE
-                  , ATOM.only = FALSE 
-                  , hex = FALSE
-                  , verbose = TRUE)
-
-rm(in_filename_mean_stability, in_filename_pdb)
-
-# assign separately for duet and ligand 
-my_pdb_duet = my_pdb
-my_pdb_lig = my_pdb
-
-#=========================================================
-# Replacing B factor with mean stability scores
-# within the respective dfs
-#==========================================================
-# extract atom list into a variable
-# since in the list this corresponds to data frame, variable will be a df
-#df_duet = my_pdb_duet[[1]]
-df_duet= my_pdb_duet[['atom']]
-df_lig = my_pdb_lig[['atom']]
-
-# make a copy: required for downstream sanity checks
-d2_duet = df_duet
-d2_lig = df_lig
-
-# sanity checks: B factor
-max(df_duet$b); min(df_duet$b)
-max(df_lig$b); min(df_lig$b)
-
-#*******************************************
-# histograms and density plots for inspection
-# 1: original B-factors
-# 2: original mean stability values
-# 3: replaced B-factors with mean stability values
-#*********************************************
-# Set the margin on all sides
-par(oma = c(3,2,3,0)
-    , mar = c(1,3,5,2)
-    #, mfrow = c(3,2)
-    , mfrow = c(3,4))
-
-#=============
-# Row 1 plots: original B-factors
-# duet and affinity
-#=============
-hist(df_duet$b
-     , xlab = "" 
-     , main = "Bfactor stability")
-
-plot(density(df_duet$b)
-     , xlab = ""
-     , main = "Bfactor stability")
-
-
-hist(df_lig$b
-     , xlab = "" 
-     , main = "Bfactor affinity")
-
-plot(density(df_lig$b)
-     , xlab = ""
-     , main = "Bfactor affinity")
-
-#=============
-# Row 2 plots: original mean stability values
-# duet and affinity
-#=============
-
-#hist(my_df$averaged_duet
-hist(my_df$avg_ens_stability_scaled
-     , xlab = "" 
-     , main = "mean stability values")
-
-#plot(density(my_df$averaged_duet)
-plot(density(my_df$avg_ens_stability_scaled)
-     , xlab = ""
-     , main = "mean stability values")
-
-#hist(my_df$averaged_affinity
-hist(my_df$avg_ens_affinity_scaled
-     , xlab = "" 
-     , main = "mean affinity values")
-
-#plot(density(my_df$averaged_affinity)
-plot(density(my_df$avg_ens_affinity_scaled)
-     , xlab = ""
-     , main = "mean affinity values")
-
-#==============
-# Row 3 plots: replaced B-factors with mean stability values
-# After actual replacement in the b factor column
-#===============
-################################################################
-#=========
-# step 0_P1: DONT RUN once you have double checked the matched output
-#=========
-# sanity check:  match and assign to a separate column to double check
-# colnames(my_df)
-# df_duet$duet_scaled = my_df$averge_duet_scaled[match(df_duet$resno, my_df$position)]
-
-#=========
-# step 1_P1
-#=========
-# Be brave and replace in place now (don"t run sanity check)
-# this makes all the B-factor values in the non-matched positions as NA
-
-#df_duet$b = my_df$averaged_duet_scaled[match(df_duet$resno, my_df$position)]
-#df_lig$b = my_df$averaged_affinity_scaled[match(df_lig$resno, my_df$position)]
-
-df_duet$b = my_df$avg_ens_stability_scaled[match(df_duet$resno, my_df$position)]
-df_lig$b  = my_df$avg_ens_affinity_scaled[match(df_lig$resno, my_df$position)]
-
-#=========
-# step 2_P1
-#=========
-# count NA in Bfactor
-b_na_duet = sum(is.na(df_duet$b)) ; b_na_duet
-b_na_lig  = sum(is.na(df_lig$b)) ; b_na_lig 
-
-# count number of 0"s in Bactor
-sum(df_duet$b == 0)
-sum(df_lig$b  == 0)
-
-# replace all NA in b factor with 0
-na_rep = 2
-df_duet$b[is.na(df_duet$b)] = na_rep
-df_lig$b[is.na(df_lig$b)] = na_rep
-
-# # sanity check: should be 0 and True
-# # duet and lig
-# if ( (sum(df_duet$b == na_rep) == b_na_duet) && (sum(df_lig$b == na_rep) == b_na_lig) ) {
-#   print ("PASS: NA's replaced with 0s successfully in df_duet and df_lig")
-# } else {
-#   print("FAIL: NA replacement in df_duet NOT successful")
-#   quit()
-# }
-# 
-# max(df_duet$b); min(df_duet$b)
-# 
-# # sanity checks: should be True
-# if( (max(df_duet$b) == max(my_df$avg_ens_stability_scaled)) & (min(df_duet$b) == min(my_df$avg_ens_stability_scaled)) ){
-#   print("PASS: B-factors replaced correctly in df_duet")
-# } else {
-#   print ("FAIL: To replace B-factors in df_duet")
-#   quit()
-# }
-
-# if( (max(df_lig$b) == max(my_df$avg_ens_affinity_scaled)) & (min(df_lig$b) == min(my_df$avg_ens_affinity_scaled)) ){
-#   print("PASS: B-factors replaced correctly in df_lig")
-# } else {
-#   print ("FAIL: To replace B-factors in df_lig")
-#   quit()
-# }
-
-#=========
-# step 3_P1
-#=========
-# sanity check: dim should be same before reassignment
-if ( (dim(df_duet)[1] == dim(d2_duet)[1]) & (dim(df_lig)[1] == dim(d2_lig)[1]) &
-     (dim(df_duet)[2] == dim(d2_duet)[2]) & (dim(df_lig)[2] == dim(d2_lig)[2])
-    ){
-  print("PASS: Dims of both dfs as expected")
-} else {
-  print ("FAIL: Dims mismatch")
-  quit()}
-
-#=========
-# step 4_P1:
-# VERY important
-#=========
-# assign it back to the pdb file
-my_pdb_duet[['atom']] = df_duet
-max(df_duet$b); min(df_duet$b)
-table(df_duet$b)
-sum(is.na(df_duet$b))
-
-my_pdb_lig[['atom']] = df_lig
-max(df_lig$b); min(df_lig$b)
-
-#=========
-# step 5_P1
-#=========
-cat(paste0("output file duet mean stability pdb:", outfile_duet_mspdb))
-write.pdb(my_pdb_duet, outfile_duet_mspdb)
-
-cat(paste0("output file ligand mean stability pdb:", outfile_lig_mspdb))
-write.pdb(my_pdb_lig, outfile_lig_mspdb)
-
-#============================
-# Add the 3rd histogram and density plots for comparisons
-#============================
-# Plots continued...
-# Row 3 plots: hist and density of replaced B-factors with stability values
-hist(df_duet$b
-     , xlab = ""
-     , main = "repalcedB duet")
-
-plot(density(df_duet$b)
-     , xlab = ""
-     , main = "replacedB duet")
-
-
-hist(df_lig$b
-     , xlab = ""
-     , main = "repalcedB affinity")
-
-plot(density(df_lig$b)
-     , xlab = ""
-     , main = "replacedB affinity")
-
-# graph titles
-mtext(text = "Frequency"
-      , side = 2
-      , line = 0
-      , outer = TRUE)
-
-mtext(text = paste0(tolower(gene), ": Stability Distribution")
-      , side = 3
-      , line = 0
-      , outer = TRUE)
-#============================================
-
-#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-# NOTE: This replaced B-factor distribution has the same
-# x-axis as the PredAff normalised values, but the distribution
-# is affected since 0 is overinflated/or hs an additional blip because
-# of the positions not associated with resistance. This is because all the positions
-# where there are no SNPs have been assigned 0???
-#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-