From d7f348318cc45a13ac054d18ace21a9d4efe9b48 Mon Sep 17 00:00:00 2001
From: Tanushree Tunstall <tanu@tunstall.in>
Date: Mon, 29 Aug 2022 23:28:53 +0100
Subject: [PATCH] removed template files from rpob

---
 .../rpob/TEMP_Ppi2_appendix_tables.R          | 460 ------------------
 .../rpob/TEMp_na_appendix_tables.R            | 441 -----------------
 2 files changed, 901 deletions(-)
 delete mode 100644 scripts/plotting/plotting_thesis/rpob/TEMP_Ppi2_appendix_tables.R
 delete mode 100644 scripts/plotting/plotting_thesis/rpob/TEMp_na_appendix_tables.R

diff --git a/scripts/plotting/plotting_thesis/rpob/TEMP_Ppi2_appendix_tables.R b/scripts/plotting/plotting_thesis/rpob/TEMP_Ppi2_appendix_tables.R
deleted file mode 100644
index bf1e0f6..0000000
--- a/scripts/plotting/plotting_thesis/rpob/TEMP_Ppi2_appendix_tables.R
+++ /dev/null
@@ -1,460 +0,0 @@
-#!/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)
diff --git a/scripts/plotting/plotting_thesis/rpob/TEMp_na_appendix_tables.R b/scripts/plotting/plotting_thesis/rpob/TEMp_na_appendix_tables.R
deleted file mode 100644
index fb89a4e..0000000
--- a/scripts/plotting/plotting_thesis/rpob/TEMp_na_appendix_tables.R
+++ /dev/null
@@ -1,441 +0,0 @@
-#!/usr/bin/env Rscript       
-
-source("~/git/LSHTM_analysis/config/gid.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_nca    = 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_nca) {
-  gene_colnames = c("mcsm_nca_affinity", "mcsm_nca_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)
-
-#--------------------
-# 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 affinity
-# naDist_colname
-####################################
-# Filtered data
-df_nca = df3_output[df3_output[[naDist_colname]]<DistCutOff,]
-
-# select cols
-cols_to_output_nca = c("position"
-                        , "sensitivity"
-                        , "mutationinformation"
-                        , naDist_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("MAF ","(%)")
-                             , "Odds Ratio"
-                             , "P-value"
-                             , "Adjusted P-value"
-                             , "Adjusted P-value significance")
-
-colnames(Out_df_ncaS) = colsNames_to_output_nca
-Out_df_ncaS
-#--------------------
-# write output file: KS test within grpup
-#----------------------
-Out_ncaT = paste0(outdir_stats
-                   , tolower(gene)
-                   , "_na_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_na_dd, mut_h_na_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$rna_site = ifelse(combined_table$position%in%aa_pos_rna, "rna", "no")
-table(combined_table$rna_site)
-
-combined_table$sam_site = ifelse(combined_table$position%in%aa_pos_sam, "sam", "no")
-table(combined_table$sam_site)
-
-combined_table$amp_site = ifelse(combined_table$position%in%aa_pos_amp, "amp", "no")
-table(combined_table$amp_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)