From e4e8bd72786feacdae9d204c6c93dc47789d79bd Mon Sep 17 00:00:00 2001
From: Tanushree Tunstall <tanu@tunstall.in>
Date: Thu, 25 Aug 2022 17:50:49 +0100
Subject: [PATCH] added colors for tiles for gid plots and appendix tables

---
 config/gid.R                                  |  14 +-
 .../plotting_thesis/ORandSNP_results.R        |  11 +-
 .../plotting_thesis/appendix_tables.R         | 290 ++++++++++++++++--
 3 files changed, 273 insertions(+), 42 deletions(-)

diff --git a/config/gid.R b/config/gid.R
index 763a429..a44c1e1 100644
--- a/config/gid.R
+++ b/config/gid.R
@@ -129,15 +129,15 @@ cat("\n==================================================="
 
 ##############################################################
 # var for position customisation for plots
-aa_pos_drug
-aa_pos_lig1 = aa_pos_sam       
-aa_pos_lig2 = aa_pos_rna
-aa_pos_lig3 = aa_pos_amp 
+aa_pos_drug = c(0)  #00ff00 # green # as STR doesn't bind
+aa_pos_lig1 = aa_pos_sam #2f4f4f # darkslategrey
+aa_pos_lig2 = aa_pos_rna #ff1493 #deeppink
+aa_pos_lig3 = aa_pos_amp #000080 #navyblue
 
-tile_map=data.frame(tile=c("SRY","SAM","RNA","AMP"),
-                    tile_colour=c("green","darkslategrey","darkred","navyblue"))
+tile_map=data.frame(tile=c("STR","SAM","RNA","AMP"),
+                    tile_colour=c("#00ff00","#2f4f4f","#ff1493","#000080"))
 
 # green: #00ff00
 # darkslategrey : #2f4f4f
-# darkred : #8b0000
+# deeppink : #ff1493
 # navyblue :#000080
\ No newline at end of file
diff --git a/scripts/plotting/plotting_thesis/ORandSNP_results.R b/scripts/plotting/plotting_thesis/ORandSNP_results.R
index 127c441..8bbef3e 100644
--- a/scripts/plotting/plotting_thesis/ORandSNP_results.R
+++ b/scripts/plotting/plotting_thesis/ORandSNP_results.R
@@ -1,9 +1,10 @@
-#!/usr/bin/env Rscript       
-#source("~/git/LSHTM_analysis/config/alr.R")
-source("~/git/LSHTM_analysis/config/embb.R")
-#source("~/git/LSHTM_analysis/config/katg.R")
-#source("~/git/LSHTM_analysis/config/gid.R")
+#!/usr/bin/env Rscript
 #source("~/git/LSHTM_analysis/config/pnca.R")
+
+#source("~/git/LSHTM_analysis/config/embb.R")
+source("~/git/LSHTM_analysis/config/gid.R")
+#source("~/git/LSHTM_analysis/config/alr.R")
+#source("~/git/LSHTM_analysis/config/katg.R")
 #source("~/git/LSHTM_analysis/config/rpob.R")
 
 # get plottting dfs 
diff --git a/scripts/plotting/plotting_thesis/appendix_tables.R b/scripts/plotting/plotting_thesis/appendix_tables.R
index 1afba7a..7e9a43a 100644
--- a/scripts/plotting/plotting_thesis/appendix_tables.R
+++ b/scripts/plotting/plotting_thesis/appendix_tables.R
@@ -1,27 +1,29 @@
 #!/usr/bin/env Rscript       
-#source("~/git/LSHTM_analysis/config/alr.R")
 source("~/git/LSHTM_analysis/config/embb.R")
-#source("~/git/LSHTM_analysis/config/katg.R")
-#source("~/git/LSHTM_analysis/config/gid.R")
-#source("~/git/LSHTM_analysis/config/pnca.R")
-#source("~/git/LSHTM_analysis/config/rpob.R")
-
-# get plottting dfs 
 source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
 
-source("~/git/LSHTM_analysis/scripts/plotting/plotting_colnames.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)
 ###################################################################
-# FIXME: ADD distance to NA when SP replies
-
 geneL_normal  = c("pnca")
-geneL_na      = c("gid", "rpob")
+#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
@@ -31,22 +33,28 @@ 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("mutationinformation"
-                   , "position"
+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]
-                   , "mcsm_ppi2_affinity"
-                   , "mcsm_ppi2_outcome"
+                   # #, affinity_dist_colnames[3]
+                   # , "mcsm_na_affinity"
+                   # , "mcsm_na_outcome"
+                   # #, "mcsm_ppi2_affinity"
+                   # #, "mcsm_ppi2_outcome"
+                   , gene_colnames
                    , "maf"
                    , "or_mychisq"
                    , "pval_fisher")
@@ -54,7 +62,7 @@ cols_to_output = c("mutationinformation"
 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
 #==============================================
@@ -91,7 +99,6 @@ head(df3_output$af); head(df3_output$maf);head(df3_output$maf_percent)
 #----------
 df3_output$pval_fisher = round(df3_output$pval_fisher,2)
 
-
 class(df3_output)
 head(df3_output)
 
@@ -100,8 +107,9 @@ head(df3_output)
 ####################################
 df_lig = df3_output[df3_output[[LigDist_colname]]<DistCutOff,]
 
-cols_to_output_lig = c("mutationinformation"
-                       , "position"
+cols_to_output_lig = c("position"
+                       , "sensitivity"
+                       , "mutationinformation"
                        , LigDist_colname
                        , "ligand_affinity_change"
                        , "ligand_outcome"
@@ -122,8 +130,9 @@ 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("Mutation"
-                            , "position"
+colsNames_to_output_lig = c("position"
+                            , "sensitivity"
+                            , "Mutation"
                             , paste0("Lig-Dist (", angstroms_symbol, ")")
                             , "mCSM-ligand affinity"
                             , "mCSM ligand_outcome"
@@ -150,18 +159,20 @@ write.csv(Out_df_ligS, Out_ligT, row.names = FALSE)
 
 ########################################################################
 ####################################
-# Appendix: PPi2 affinity
+# Appendix: ppi2 affinity
 ####################################
 
 # Filtered data
 df_ppi2 = df3_output[df3_output[[ppi2Dist_colname]]<DistCutOff,]
 
 # select cols
-cols_to_output_ppi2 = c("mutationinformation"
-                        , "position"
+cols_to_output_ppi2 = c("position"
+                        , "sensitivity"
+                        , "mutationinformation"
                         , ppi2Dist_colname
-                        , "mcsm_ppi2_affinity"
-                        , "mcsm_ppi2_outcome"
+                        #, "mcsm_ppi2_affinity"
+                        #, "mcsm_ppi2_outcome"
+                        , gene_colnames
                         , "maf_percent"
                         , "or_mychisq"
                         , "pval_fisher"
@@ -176,8 +187,9 @@ Out_df_ppi2 = df_ppi2[, cols_to_output_ppi2]
 Out_df_ppi2S = Out_df_ppi2[order(-Out_df_ppi2$or_mychisq, Out_df_ppi2$maf_percent), ]
 
 
-colsNames_to_output_ppi2 = c("Mutation"
-                             , "position"
+colsNames_to_output_ppi2 = c("position"
+                             , "sensitivity"
+                             , "Mutation"
                              , paste0("PPI2-Dist (", angstroms_symbol, ")")
                              , paste0("mCSM-PPI2 (", delta_symbol,delta_symbol,"G)")
                              , "mCSM-PPI2 outcome"
@@ -199,3 +211,221 @@ Out_ppi2T = paste0(outdir_stats
 cat("Output of PPI2 muts:", Out_ppi2T )
 write.csv(Out_df_ppi2S, Out_ppi2T, 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
+#-------------------
+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 = "Most 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_dd$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")
+}
+
+#--------------------
+# 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)
+