diff --git a/scripts/plotting/plotting_thesis/embb/embb_ORandSNP_results.R b/scripts/plotting/plotting_thesis/embb/embb_ORandSNP_results.R
new file mode 100644
index 0000000..4716963
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/embb/embb_ORandSNP_results.R
@@ -0,0 +1,223 @@
+#!/usr/bin/env Rscript
+#source("~/git/LSHTM_analysis/config/katg.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/") 
+
+###################################################################
+# 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
+###############################################
+# OR plot
+df3_or = merged_df3
+df3_or$maf_percent = df3_or$maf*100
+
+bar = df3_or[, c("mutationinformation"
+                     , "wild_pos"
+                     , "position"
+                     , "sensitivity"
+                     , drug
+                     , affinity_dist_colnames
+                     , "maf_percent"
+                     , "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"
+               , drug
+               , affinity_dist_colnames
+               , "maf_percent"
+               , "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("\nNo.of DST muts:", nrow(bar_or)
+    , "\nNo of DST (R):", table(bar_or$sensitivity)[[1]]
+    , "\nNo of DST (S):", table(bar_or$sensitivity)[[2]]
+    , "\nNumber of R muts with OR >1 (n = ", nrow(res1),")"
+    , "\nPercentage of muts with OR>1 i.e resistant:" , pc_orR *100 )
+
+table(bar_or$sensitivity)
+
+# 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$heme_site = ifelse(bar_or$position%in%aa_pos_hem, "heme", "no")
+table(bar_or$heme_site)
+
+top10_or = bar_or[1:10,]
+top10_or
+
+write.csv(bar_or, paste0(outdir_stats, "katg_OR_10.csv"))
+
+# are these active sites
+top10_or$position[top10_or$position%in%active_aa_pos]
+
+
+# maf
+bar_maf = bar_or[order(bar_or$maf_percent
+                      , bar_or$ligand_distance
+                      # bar_or$nca_dist
+                      , bar_or$interface_dist
+                      , decreasing = T), ]
+
+head(bar_maf)
+#########################################################
+# closest 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/embb/embb_appendix_tables.R b/scripts/plotting/plotting_thesis/embb/embb_appendix_tables.R
new file mode 100644
index 0000000..bf1e0f6
--- /dev/null
+++ b/scripts/plotting/plotting_thesis/embb/embb_appendix_tables.R
@@ -0,0 +1,460 @@
+#!/usr/bin/env Rscript       
+#source("~/git/LSHTM_analysis/config/katg.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$heme_site = ifelse(Out_df_ligS$position%in%aa_pos_hem, "heme", "no")
+table(Out_df_ligS$heme_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/PPI2 affinity
+# naDist_colname
+# ppi2Dist_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$heme_site = ifelse(Out_df_ncaS$position%in%aa_pos_hem, "heme", "no")
+table(Out_df_ncaS$heme_site)
+
+#--------------------
+# write output file: KS test within grpup
+#----------------------
+Out_ncaT = paste0(outdir_stats
+                   , tolower(gene)
+                   , "_ppi2_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: colnames ppi2 mismatch")
+  }
+  
+  if (identical(colnames(mut_h_ppi2_dd), colnames(mut_h_lig_dd)) ){
+    cat("\nPass 2 : ppi2")
+  }else{
+    quit("Abort: 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$heme_site = ifelse(combined_table$position%in%aa_pos_hem, "heme", "no")
+table(combined_table$heme_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)