diff --git a/scripts/plotting/plotting_thesis/dm_om_plots.R b/scripts/plotting/plotting_thesis/dm_om_plots.R
new file mode 100644
index 0000000..937aadf
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/dm_om_plots.R
@@ -0,0 +1,222 @@
+geneL_normal  = c("pnca")
+geneL_na      = c("gid", "rpob")
+geneL_ppi2    = c("alr", "embb", "katg", "rpob")
+
+
+if (tolower(gene)%in%geneL_na){
+  lf_mcsm_na
+  }
+
+if (tolower(gene)%in%geneL_ppi2){
+  lf_mcsm_ppi2
+}
+
+colnames(lf_duet)
+table(lf_duet$param_type)
+
+static_colsP = c("Lig Dist(Å)","ASA", "RSA","RD","KD","Log10(MAF)")
+
+stability_suffix <- paste0(delta_symbol, delta_symbol, "G")
+ 
+lf_commonP   = lf_duet[!lf_duet$param_type%in%c("DUET ΔΔG"),]
+lf_commonP$param_type = levels(droplevels(lf_commonP$param_type))
+table(lf_commonP$param_type); colnames(lf_commonP)
+lf_commonP$outcome = lf_commonP$duet_outcome
+lf_commonP$duet_outcome = NULL
+
+
+lf_duet$outcome = lf_duet$duet_outcome
+lf_duet$duet_outcome = NULL
+lf_duetP                = lf_duet[!lf_duet$param_type%in%c(static_colsP, "outcome"),]
+lf_duetP$param_type     = levels(droplevels(lf_duetP$param_type))
+table(lf_duetP$param_type); colnames(lf_duetP)
+colnames(lf_duetP)
+
+lf_foldx$outcome = lf_foldx$foldx_outcome
+lf_foldx$foldx_outcome = NULL
+lf_foldxP               = lf_foldx[!lf_foldx$param_type%in%c(static_colsP,"outcome"),]
+lf_foldxP$param_type    = levels(droplevels(lf_foldxP$param_type))
+table(lf_foldxP$param_type) 
+colnames(lf_foldxP)
+
+
+lf_deepddg$outcome = lf_deepddg$deepddg_outcome
+lf_deepddg$deepddg_outcome = NULL
+lf_deepddgP             = lf_deepddg[!lf_deepddg$param_type%in%c(static_colsP, "outcome"),]
+lf_deepddgP$param_type  = levels(droplevels(lf_deepddgP$param_type))
+table(lf_deepddgP$param_type)
+colnames(lf_deepddgP)
+
+
+lf_dynamut2$outcome = lf_dynamut2$ddg_dynamut2_outcome
+lf_dynamut2$ddg_dynamut2_outcome = NULL
+lf_dynamut2P            = lf_dynamut2[!lf_dynamut2$param_type%in%c(static_colsP, "outcome"),]
+lf_dynamut2P$param_type = levels(droplevels(lf_dynamut2P$param_type))
+table(lf_dynamut2P$param_type)
+colnames(lf_dynamut2P)
+
+
+lf_consurf$outcome = lf_consurf$consurf_outcome
+lf_consurf$consurf_outcome = NULL
+lf_consurfP             = lf_consurf[!lf_consurf$param_type%in%c(static_colsP),]
+lf_consurfP$param_type  = levels(droplevels(lf_consurfP$param_type))
+table(lf_consurfP$param_type)
+colnames(lf_consurfP)
+
+
+lf_snap2$outcome = lf_snap2$snap2_outcome
+lf_snap2$snap2_outcome = NULL
+lf_snap2P               = lf_snap2[!lf_snap2$param_type%in%c(static_colsP),]
+lf_snap2P$param_type    = levels(droplevels(lf_snap2P$param_type))
+table(lf_snap2P$param_type)
+colnames(lf_snap2P)
+
+
+lf_provean$outcome = lf_provean$provean_outcome
+lf_provean$provean_outcome = NULL
+lf_proveanP             = lf_provean[!lf_provean$param_type%in%c(static_colsP),]
+lf_proveanP$param_type  = levels(droplevels(lf_proveanP$param_type))
+table(lf_proveanP$param_type)
+colnames(lf_proveanP)
+
+bar = rbind(colnames(lf_duetP)
+               , colnames(lf_foldxP)
+               , colnames(lf_deepddgP)
+               , colnames(lf_dynamut2P)
+               , colnames(lf_consurfP)
+               , colnames(lf_snap2P)
+               , colnames(lf_proveanP)
+)
+bar
+
+lf_df_stabP = rbind((lf_duetP)
+                      , (lf_foldxP)
+                      , (lf_deepddgP)
+                      , (lf_dynamut2P))
+
+lf_df_consP = rbind((lf_consurfP)
+                      , (lf_snap2P)
+                      , (lf_proveanP))
+
+table(lf_df_stabP$param_type)
+# VERY USEFUL for seeing numbers for param types
+table(lf_df_stabP$param_type,lf_df_stabP$outcome)
+table(lf_df_consP$param_type,lf_df_consP$outcome)
+
+#==============
+# Plot:BP
+#==============
+stability_suffix <- paste0(delta_symbol, delta_symbol, "G")
+
+# lf_bp(lf_df_stabP, p_title = paste0("Stability",stability_suffix)
+#       , violin_quantiles = c(0.5))
+# lf_bp(lf_df_consP, p_title = "Evolutionary Conservation"
+# , violin_quantiles = c(0.5))
+
+lf_bp2(lf_df_stabP, p_title = paste0("Stability",stability_suffix)
+       , violin_quantiles = c(0.5), monochrome = F)
+
+lf_bp2(lf_df_consP, p_title = "Evolutionary Conservation"
+      , violin_quantiles = c(0.5), monochrome = F)
+
+#HMMM: Bollocks!
+lf_bp2(lf_commonP,  p_title = paste0("Residue level properties")
+       , violin_quantiles = c(0.5)
+      , monochrome = T) # doesn't plot stat bars
+
+lf_bp(lf_commonP,  p_title = paste0("Residue level properties")
+       , violin_quantiles = c(0.5)) #plots stat bars but incorrect result
+
+lf_unpaired_stats(lf_duet)
+wilcox.test(wf_duet$`Lig Dist(Å)`[wf_duet$mutation_info_labels=="R"]
+            , wf_duet$`Lig Dist(Å)`[wf_duet$mutation_info_labels=="S"])
+
+wilcox.test(wf_duet$ASA[wf_duet$mutation_info_labels=="R"]
+            , wf_duet$ASA[wf_duet$mutation_info_labels=="S"])
+
+
+# 1: variable
+# 9: conserved
+
+# CHECK THESE
+foo = merged_df3[c("dst_mode", "mutation_info_labels", "consurf_colour_rev"
+                   , "consurf_scaled"
+                   , "consurf_score"
+                   , "consurf_outcome"
+                   , "snap2_score"
+                   , "snap2_scaled"
+                   , "snap2_outcome"
+                   , "provean_score"
+                   , "provean_scaled"
+                   , "provean_outcome")]
+
+################
+# Affinity
+################
+# ligand
+lf_mcsm_lig$outcome = lf_mcsm_lig$ligand_outcome
+lf_mcsm_lig$ligand_outcome = NULL
+colnames(lf_mcsm_lig)
+table(lf_mcsm_lig$param_type)
+
+
+lf_mcsm_lig$outcome = lf_mcsm_lig$ligand_outcome
+lf_mcsm_ligP            = lf_mcsm_lig[!lf_mcsm_lig$param_type%in%c(static_colsP, "outcome"),]
+#lf_mcsm_ligP$param_type = levels(droplevels(lf_mcsm_ligP))
+table(lf_mcsm_ligP$param_type)
+lf_mcsm_ligP$ligand_outcome = NULL
+colnames(lf_mcsm_ligP)
+
+if (tolower(gene)%in%geneL_na){
+  lf_mcsm_na$outcome = lf_mcsm_na$mcsm_na_outcome
+  #lf_mcsm_na$mcsm_na_outcome = NULL
+  lf_mcsm_naP            = lf_mcsm_na[!lf_mcsm_na$param_type%in%c(static_colsP, "outcome"),]
+  #lf_mcsm_naP$param_type = levels(droplevels(lf_mcsm_naP))
+  table(lf_mcsm_naP$param_type)
+  lf_mcsm_naP$mcsm_na_outcome = NULL
+  colnames(lf_mcsm_naP)
+}
+
+if (tolower(gene)%in%geneL_ppi2){
+  
+  lf_mcsm_ppi2$outcome = lf_mcsm_ppi2$mcsm_ppi2_outcome
+  colnames(lf_mcsm_ppi2)
+  #lf_mcsm_ppi2$mcsm_ppi2_outcome = NULL
+  lf_mcsm_ppi2P           = lf_mcsm_ppi2[!lf_mcsm_ppi2$param_type%in%c(static_colsP, "outcome"),]
+  #lf_mcsm_ppi2P$param_type = levels(droplevels(lf_mcsm_ppi2P))
+  table(lf_mcsm_ppi2P$param_type)
+  lf_mcsm_ppi2P$mcsm_ppi2_outcome = NULL
+  colnames(lf_mcsm_ppi2P)
+  
+}
+
+
+bar = rbind(colnames(lf_mcsm_ligP)
+            #, colnames(lf_mcsm_naP)
+            , colnames(lf_mcsm_ppi2P))
+bar
+
+lf_df_affP = rbind((lf_mcsm_ligP)
+                    , (lf_mcsm_ppi2P))
+
+lf_bp(lf_df_affP,  p_title = paste0("Affinity changes")
+       , violin_quantiles = c(0.5))
+       #, monochrome = T) # doesn't plot stat bars
+
+wilcox.test(wf_mcsm_lig$ASA[wf_mcsm_lig$mutation_info_labels=="R"]
+            , wf_mcsm_lig$ASA[wf_mcsm_lig$mutation_info_labels=="S"])
+
+wilcox.test(wf_mcsm_ppi2$ASA[wf_mcsm_ppi2$mutation_info_labels=="R"]
+            , wf_mcsm_ppi2$ASA[wf_mcsm_ppi2$mutation_info_labels=="S"])
+
+#===============================
+p1 =   lf_bp2(lf_df_stabP, p_title = paste0("Stability",stability_suffix)
+       , violin_quantiles = c(0.5), monochrome = F)
+
+p2 = lf_bp2(lf_df_consP, p_title = "Evolutionary Conservation"
+       , violin_quantiles = c(0.5), monochrome = F)
+
+p3 = lf_bp2(lf_df_affP,  p_title = paste0("Affinity changes")
+      , violin_quantiles = c(0.5), monochrome = F)
+
+cowplot::plot_grid(p1,cowplot::plot_grid(p2, p3), nrow=2)