diff --git a/scripts/plotting/mcsm_affinity_data_only.R b/scripts/plotting/mcsm_affinity_data_only.R
new file mode 100644
index 0000000..cb20c00
--- /dev/null
+++ b/scripts/plotting/mcsm_affinity_data_only.R
@@ -0,0 +1,243 @@
+#source("~/git/LSHTM_analysis/config/pnca.R")
+#source("~/git/LSHTM_analysis/config/alr.R")
+#source("~/git/LSHTM_analysis/config/gid.R")
+#source("~/git/LSHTM_analysis/config/embb.R")
+#source("~/git/LSHTM_analysis/config/katg.R")
+source("~/git/LSHTM_analysis/config/rpob.R")
+
+source("/home/tanu/git/LSHTM_analysis/my_header.R")
+#########################################################
+# TASK: Generate averaged affinity values 
+# across all affinity tools for a given structure
+# as applicable...
+#########################################################
+
+#=======
+# output
+#=======
+outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
+
+#OutFile1
+outfile_mean_aff = paste0(outdir_images, "/", tolower(gene)
+                                , "_mean_affinity_all.csv")
+print(paste0("Output file:", outfile_mean_aff))
+
+#OutFile2
+outfile_mean_aff_priorty = paste0(outdir_images, "/", tolower(gene)
+                          , "_mean_affinity_priority.csv")
+print(paste0("Output file:", outfile_mean_aff_priorty))
+
+#%%===============================================================
+
+#=============
+# Input
+#=============
+df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv")
+df3 = read.csv(df3_filename)
+length(df3$mutationinformation)
+
+# mut_info checks
+table(df3$mutation_info)
+table(df3$mutation_info_orig)
+table(df3$mutation_info_labels_orig)
+
+# used in plots and analyses
+table(df3$mutation_info_labels) # different, and matches dst_mode
+table(df3$dst_mode)
+
+# create column based on dst mode with different colname
+table(is.na(df3$dst))
+table(is.na(df3$dst_mode))
+
+#===============
+# Create column: sensitivity mapped to dst_mode
+#===============
+df3$sensitivity = ifelse(df3$dst_mode == 1, "R", "S")
+table(df3$sensitivity)
+
+length(unique((df3$mutationinformation)))
+all_colnames = as.data.frame(colnames(df3))
+common_cols  = c("mutationinformation"
+                 , "X5uhc_position"
+                 , "X5uhc_offset"
+                 , "position"
+                 , "dst_mode"
+                 , "mutation_info_labels"
+                 , "sensitivity"
+                 , "ligand_distance"
+                 , "interface_dist")
+
+all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)]
+all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)]
+
+#===================
+# stability cols
+#===================
+# raw_cols_stability =  c("duet_stability_change"
+#                         , "deepddg"
+#                         , "ddg_dynamut2"
+#                         , "ddg_foldx")
+# 
+# scaled_cols_stability = c("duet_scaled"       
+#                           , "deepddg_scaled"   
+#                           , "ddg_dynamut2_scaled"
+#                           , "foldx_scaled")
+# 
+# outcome_cols_stability = c("duet_outcome"
+#                            , "deepddg_outcome"
+#                            , "ddg_dynamut2_outcome"
+#                            , "foldx_outcome")
+
+#===================
+# affinity cols
+#===================
+raw_cols_affinity =  c("ligand_affinity_change"
+                       , "mmcsm_lig"
+                       , "mcsm_ppi2_affinity"
+                       , "mcsm_na_affinity")
+
+scaled_cols_affinity = c("affinity_scaled" 
+                         , "mmcsm_lig_scaled" 
+                         , "mcsm_ppi2_scaled" 
+                         , "mcsm_na_scaled" )
+
+outcome_cols_affinity  = c( "ligand_outcome"
+                            , "mmcsm_lig_outcome"
+                            , "mcsm_ppi2_outcome"
+                            , "mcsm_na_outcome")
+
+#===================
+# conservation cols
+#===================
+# raw_cols_conservation =  c("consurf_score"
+#                            , "snap2_score"
+#                            , "provean_score")
+# 
+# scaled_cols_conservation = c("consurf_scaled"
+#                              , "snap2_scaled"
+#                              , "provean_scaled")
+# 
+# # CANNOT strictly be used, as categories are not identical with conssurf missing altogether
+# outcome_cols_conservation = c("provean_outcome"
+#                               , "snap2_outcome"
+#                               #consurf outcome doesn't exist
+# )
+
+######################################################################
+cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols
+                                                    , raw_cols_affinity
+                                                    , scaled_cols_affinity
+                                                    , outcome_cols_affinity
+                                                    # , raw_cols_stability
+                                                    # , scaled_cols_stability
+                                                    # , outcome_cols_stability             
+                                                    )]
+
+cols_to_extract  = cols_to_consider[cols_to_consider%in%c(common_cols
+                                                          , raw_cols_affinity
+                                                          , scaled_cols_affinity)]
+
+df3_plot = df3[, cols_to_extract]
+
+DistCutOff_colnames = c("ligand_distance", "interface_dist")
+DistCutOff = 10
+
+df3_affinity_filtered = df3_plot[ (df3_plot$ligand_distance<10 | df3_plot$interface_dist <10),]
+c0u = unique(df3_affinity_filtered$position)
+length(c0u)
+
+foo = df3_affinity_filtered[df3_affinity_filtered$ligand_distance<10,]
+bar = df3_affinity_filtered[df3_affinity_filtered$interface_dist<10,]
+
+wilcox.test(foo$mmcsm_lig_scaled~foo$sensitivity)
+wilcox.test(foo$mmcsm_lig~foo$sensitivity)
+
+wilcox.test(foo$affinity_scaled~foo$sensitivity)
+wilcox.test(foo$ligand_affinity_change~foo$sensitivity)
+
+wilcox.test(bar$mcsm_na_scaled~bar$sensitivity)
+wilcox.test(bar$mcsm_na_affinity~bar$sensitivity)
+
+wilcox.test(bar$mcsm_ppi2_scaled~bar$sensitivity)
+wilcox.test(bar$mcsm_ppi2_affinity~bar$sensitivity)
+
+##############################################################
+df = df3_affinity_filtered
+sum(is.na(df))
+df2 = na.omit(df)                           # Apply na.omit function
+
+a = df2[df2$position==37,]
+sel_cols = c("mutationinformation", "position", scaled_cols_affinity)
+a = a[, sel_cols]
+
+##############################################################
+# FIXME: ADD distance to NA when SP replies
+
+#####################
+# Ensemble affinity: affinity_cols
+# mcsm_lig, mmcsm_lig and mcsm_na
+#####################
+# extract outcome cols and map numeric values to the categories
+# Destabilising == 0, and stabilising == 1 so rescaling can let -1 be destabilising
+#########################################
+#=====================================
+# Affintiy (2 cols): average the scores
+# across predictors ==> average by
+# position ==> scale b/w -1 and 1
+
+# column to average: ens_affinity
+#=====================================
+cols_mcsm_lig  = c("mutationinformation"
+                   , "position"
+                   , "sensitivity"
+                   , "X5uhc_position"
+                   , "X5uhc_offset"
+                   , "ligand_distance"
+                   , "ligand_outcome"
+                   , "mmcsm_lig_outcome")
+
+
+
+
+
+
+
+
+
+######################################################################
+##################
+# merge: mean ensemble stability and affinity by_position
+####################
+# if ( class(mean_ens_stability_by_position) && class(mean_ens_affinity_by_position) != "data.frame"){
+#   cat("Y")
+# }
+
+common_cols = intersect(colnames(mean_ens_stability_by_position), colnames(mean_ens_affinity_by_position))
+
+if (dim(mean_ens_stability_by_position) && dim(mean_ens_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_ens_stability_by_position
+                                    , mean_ens_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_ens_st_aff
+          , row.names = F)
+cat("Finished writing file:\n"
+    , outfile_mean_ens_st_aff
+    , "\nNo. of rows:", nrow(combined_df)
+    , "\nNo. of cols:", ncol(combined_df))
+
+# end of script
+#===============================================================
diff --git a/scripts/plotting/mcsm_mean_stability_ensemble_5uhc_rpob.R b/scripts/plotting/mcsm_mean_stability_ensemble_5uhc_rpob.R
new file mode 100644
index 0000000..c75a6ba
--- /dev/null
+++ b/scripts/plotting/mcsm_mean_stability_ensemble_5uhc_rpob.R
@@ -0,0 +1,176 @@
+#source("~/git/LSHTM_analysis/config/pnca.R")
+#source("~/git/LSHTM_analysis/config/alr.R")
+#source("~/git/LSHTM_analysis/config/gid.R")
+#source("~/git/LSHTM_analysis/config/embb.R")
+#source("~/git/LSHTM_analysis/config/katg.R")
+#source("~/git/LSHTM_analysis/config/rpob.R")
+
+source("/home/tanu/git/LSHTM_analysis/my_header.R")
+#########################################################
+# TASK: Generate averaged stability values 
+# across all stability tools
+# for a given structure
+#########################################################
+
+#=======
+# output
+#=======
+outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
+outfile_mean_ens_st_aff = paste0(outdir_images, "/5uhc_", tolower(gene)
+                                , "_mean_ens_stability.csv")
+print(paste0("Output file:", outfile_mean_ens_st_aff))
+
+#%%===============================================================
+
+#=============
+# Input
+#=============
+df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv")
+df3 = read.csv(df3_filename)
+length(df3$mutationinformation)
+
+# mut_info checks
+table(df3$mutation_info)
+table(df3$mutation_info_orig)
+table(df3$mutation_info_labels_orig)
+
+# used in plots and analyses
+table(df3$mutation_info_labels) # different, and matches dst_mode
+table(df3$dst_mode)
+
+# create column based on dst mode with different colname
+table(is.na(df3$dst))
+table(is.na(df3$dst_mode))
+
+#===============
+# Create column: sensitivity mapped to dst_mode
+#===============
+df3$sensitivity = ifelse(df3$dst_mode == 1, "R", "S")
+table(df3$sensitivity)
+
+length(unique((df3$mutationinformation)))
+all_colnames = as.data.frame(colnames(df3))
+common_cols  = c("mutationinformation"
+                 , "X5uhc_position"
+                 , "dst_mode"
+                 , "mutation_info_labels"
+                 , "sensitivity"
+                 , "X5uhc_position"
+                 , "X5uhc_offset"
+                 , "ligand_distance"
+                 , "interface_dist")
+
+all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)]
+all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)]
+
+#===================
+# stability cols
+#===================
+raw_cols_stability =  c("duet_stability_change"
+                        , "deepddg"
+                        , "ddg_dynamut2"
+                        , "ddg_foldx")
+
+scaled_cols_stability = c("duet_scaled"       
+                          , "deepddg_scaled"   
+                          , "ddg_dynamut2_scaled"
+                          , "foldx_scaled")
+
+outcome_cols_stability = c("duet_outcome"
+                           , "deepddg_outcome"
+                           , "ddg_dynamut2_outcome"
+                           , "foldx_outcome")
+
+###########################################################
+cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols
+                                                    , raw_cols_stability
+                                                    , scaled_cols_stability
+                                                    , outcome_cols_stability)]
+
+cols_to_extract  = cols_to_consider[cols_to_consider%in%c(common_cols
+                                                          , outcome_cols_stability)]
+##############################################################
+#####################
+# Ensemble stability: outcome_cols_stability
+#####################
+# extract outcome cols and map numeric values to the categories
+# Destabilising == 0, and stabilising == 1, so rescaling can let -1 be destabilising
+df3_plot = df3[, cols_to_extract]
+
+# assign numeric values to outcome
+df3_plot[, outcome_cols_stability] <- sapply(df3_plot[, outcome_cols_stability]
+                             , function(x){ifelse(x == "Destabilising", 0, 1)})
+table(df3$duet_outcome)
+table(df3_plot$duet_outcome)
+#=====================================
+# Stability (4 cols): average the scores
+# across predictors ==> average by
+# X5uhc_position ==> scale b/w -1 and 1
+
+# column to average: ens_stability
+#=====================================
+cols_to_average = which(colnames(df3_plot)%in%outcome_cols_stability)
+
+# ensemble average across predictors
+df3_plot$ens_stability = rowMeans(df3_plot[,cols_to_average])
+
+head(df3_plot$X5uhc_position); head(df3_plot$mutationinformation)
+head(df3_plot$ens_stability)
+table(df3_plot$ens_stability)
+
+# ensemble average of predictors by X5uhc_position
+mean_ens_stability_by_position <- df3_plot %>%
+  dplyr::group_by(X5uhc_position) %>%
+  dplyr::summarize(avg_ens_stability = mean(ens_stability))
+
+# REscale b/w -1 and 1
+#en_stab_min = min(mean_ens_stability_by_position['avg_ens_stability'])
+#en_stab_max = max(mean_ens_stability_by_position['avg_ens_stability']) 
+
+# scale the average stability value between -1 and 1
+# mean_ens_by_position['averaged_stability3_scaled'] = lapply(mean_ens_by_position['averaged_stability3']
+#                                                        , function(x) ifelse(x < 0, x/abs(en3_min), x/en3_max))
+
+mean_ens_stability_by_position['avg_ens_stability_scaled'] = lapply(mean_ens_stability_by_position['avg_ens_stability']
+                                                                  , function(x) {
+                                                                    scales::rescale(x, to  = c(-1,1)
+                                                                                    #, from = c(en_stab_min,en_stab_max))
+                                                                                    , from = c(0,1))
+                                                                  })
+cat(paste0('Average stability scores:\n'
+           , head(mean_ens_stability_by_position['avg_ens_stability'])
+           , '\n---------------------------------------------------------------'
+           , '\nAverage stability scaled scores:\n'
+           , head(mean_ens_stability_by_position['avg_ens_stability_scaled'])))
+
+# convert to a data frame
+mean_ens_stability_by_position = as.data.frame(mean_ens_stability_by_position)
+
+#FIXME: sanity checks
+# TODO: predetermine the bounds
+# l_bound_ens = min(mean_ens_stability_by_position['avg_ens_stability_scaled'])
+# u_bound_ens = max(mean_ens_stability_by_position['avg_ens_stability_scaled'])
+# 
+# if ( (l_bound_ens == -1) && (u_bound_ens == 1) ){
+#   cat(paste0("PASS: ensemble stability scores averaged by X5uhc_position and then scaled"
+#         , "\nmin ensemble averaged stability: ", l_bound_ens
+#         , "\nmax ensemble averaged stability: ", u_bound_ens))
+# }else{
+#   cat(paste0("FAIL: avergaed duet scores could not be scaled b/w -1 and 1"
+#         , "\nmin ensemble averaged stability: ", l_bound_ens
+#         , "\nmax ensemble averaged stability: ", u_bound_ens))
+#   quit()
+# } 
+##################################################################
+# output
+#write.csv(combined_df, outfile_mean_ens_st_aff
+write.csv(mean_ens_stability_by_position
+          , outfile_mean_ens_st_aff
+          , row.names = F)
+cat("Finished writing file:\n"
+    , outfile_mean_ens_st_aff
+    , "\nNo. of rows:", nrow(mean_ens_stability_by_position)
+    , "\nNo. of cols:", ncol(mean_ens_stability_by_position))
+
+# end of script
+#===============================================================
diff --git a/scripts/plotting/replaceBfactor_pdb_stability_5uhc_rpob.R b/scripts/plotting/replaceBfactor_pdb_stability_5uhc_rpob.R
new file mode 100644
index 0000000..25cb705
--- /dev/null
+++ b/scripts/plotting/replaceBfactor_pdb_stability_5uhc_rpob.R
@@ -0,0 +1,283 @@
+#!/usr/bin/env Rscript                                                  
+
+#########################################################
+# TASK: Replace B-factors in the pdb file with the mean
+# normalised stability values.
+
+# read pdb file
+
+# 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/5uhc_rpob_", tolower(gene))
+
+#======
+# input
+#======
+#in_filename_pdb = paste0(tolower(gene), "_complex.pdb") 
+in_filename_pdb = "5uhc.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_stability.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))
+
+#%%===============================================================
+#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
+
+#=========================================================
+# 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']]
+
+# make a copy: required for downstream sanity checks
+d2_duet = df_duet
+
+# sanity checks: B factor
+max(df_duet$b); min(df_duet$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))
+    , mfrow  = c(3,2))
+
+
+#=============
+# 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")
+
+#=============
+# 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")
+
+#==============
+# 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_duet$b = my_df$avg_ens_stability_scaled[match(df_duet$resno, my_df$position)]
+
+#=========
+# step 2_P1
+#=========
+# count NA in Bfactor
+b_na_duet = sum(is.na(df_duet$b)) ; b_na_duet
+
+# count number of 0"s in Bactor
+sum(df_duet$b == 0)
+
+# replace all NA in b factor with 0
+na_rep = 2
+df_duet$b[is.na(df_duet$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_duet)[2] == dim(d2_duet)[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))
+
+#=========
+# step 5_P1
+#=========
+cat(paste0("output file duet mean stability pdb:", outfile_duet_mspdb))
+write.pdb(my_pdb_duet, outfile_duet_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")
+
+# 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???
+#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
\ No newline at end of file