diff --git a/scripts/functions/bp_subcolours.R b/scripts/functions/bp_subcolours.R
index a3cc403..3db4079 100755
--- a/scripts/functions/bp_subcolours.R
+++ b/scripts/functions/bp_subcolours.R
@@ -3,7 +3,7 @@
 # LINK: https://stackoverflow.com/questions/49818271/stacked-barplot-with-colour-gradients-for-each-bar
 #########################################################
 
-ColourPalleteMulti <- function(df, group, subgroup){
+ColourPalleteMulti = function(df, group, subgroup){
   
   # Find how many colour categories to create and the number of colours in each
   categories <- aggregate(as.formula(paste(subgroup, group, sep="~" ))
@@ -24,4 +24,81 @@ ColourPalleteMulti <- function(df, group, subgroup){
                                                          , category.end[i]))(categories[i,2])}))
   return(colours)
 }
-#########################################################
\ No newline at end of file
+#########################################################################
+
+bp_stability_hmap <- function(plotdf = merged_df3
+                              , xvar_colname = "position"
+                              #, bar_col_colname = "group"
+                              , stability_colname = "duet_scaled"
+                              , stability_outcome_colname = "duet_outcome"
+                              , p_title = ""  # "Protein stability (DUET)"
+                              , my_xaxls = 12 # x-axis label size
+                              , my_yaxls = 20 # y-axis label size
+                              , my_xaxts = 18 # x-axis text size
+                              , my_yaxts = 20 # y-axis text size
+                              , my_pts  = 20  # plot-title size
+                              , my_xlab = "Position"
+                              , my_ylab = "No. of nsSNPs"
+)
+{
+  
+  # order the df by position and ensure it is a factor
+  plotdf = plotdf[order(plotdf[[xvar_colname]]), ]
+  plotdf[[xvar_colname]] = factor(plotdf[[xvar_colname]])
+  
+  #cat("\nSneak peak:\n")
+  head(data.frame( plotdf[[xvar_colname]], plotdf[[stability_colname]] ) )
+  
+  # stability values isolated to help with generating column called: 'group'
+  my_grp = plotdf[[stability_colname]]
+  cat( "\nLength of nsSNPs:", length(my_grp)
+       , "\nLength of unique values for nsSNPs:", length(unique(my_grp)) )
+  
+  # Add col: 'group'
+  plotdf$group = paste0(plotdf[[stability_outcome_colname]], "_", my_grp, sep = "")
+  
+  # check unique values in normalised data
+  cat("\nNo. of unique values in", stability_colname, "no rounding:"
+      , length(unique(plotdf[[stability_colname]])))
+  
+  # Call the function to create the palette based on the group defined above
+  #subcols_ps
+  subcols_bp_hmap = ColourPalleteMulti(plotdf, stability_outcome_colname, stability_colname)
+  
+  cat("\nNo. of sub colours generated:", length(subcols_bp_hmap))
+  
+  #-------------------------------
+  # Generate the subcols barplot
+  #-------------------------------
+  
+  #g = ggplot(plotdf, aes(x = factor(position, ordered = T)))
+  g = ggplot(plotdf, aes_string(x = xvar_colname
+                                # , ordered = T)
+  ))
+  
+  
+  OutWidePlot = g + geom_bar(aes(fill = group)
+                             , colour = "grey") +
+    
+    scale_fill_manual( values = subcols_bp_hmap
+                       , guide = "none") +
+    
+    theme( axis.text.x = element_text(size = my_xaxls
+                                      , angle = 90
+                                      , hjust = 1
+                                      , vjust = 0.4)
+           , axis.text.y = element_text(size = my_yaxls
+                                        , angle = 0
+                                        , hjust = 1
+                                        , vjust = 0)
+           , axis.title.x = element_text(size = my_xaxts)
+           , axis.title.y = element_text(size = my_yaxts ) 
+           , plot.title = element_text(size = my_pts
+                                       , hjust = 0.5)) +
+    
+    labs(title = p_title
+         , x = my_xlab
+         , y = my_ylab)
+  
+  return(OutWidePlot)
+}
diff --git a/scripts/plotting/get_plotting_dfs.R b/scripts/plotting/get_plotting_dfs.R
index c1ce5b2..a50b0a9 100755
--- a/scripts/plotting/get_plotting_dfs.R
+++ b/scripts/plotting/get_plotting_dfs.R
@@ -129,7 +129,7 @@ merged_df3_comp = all_plot_dfs[[4]]
 #                        Data for subcols barplot (~heatmap)
 ####################################################################
 
-source("coloured_bp_data.R")
+#source("coloured_bp_data.R")
 
 ####################################################################
 #                        Data for logoplots