From 06e53631125f378bfcdbab38c14108e6346fece8 Mon Sep 17 00:00:00 2001
From: Tanushree Tunstall <tanu@tunstall.in>
Date: Sun, 31 Jul 2022 16:29:47 +0100
Subject: [PATCH] added scriptaa mcsm_mean_stability_ensemble.R to get ensemble
 of averages across predictors for stability and affinity

---
 .../plotting/mcsm_mean_stability_ensemble.R   | 280 ++++++++++++++++++
 1 file changed, 280 insertions(+)
 create mode 100644 scripts/plotting/mcsm_mean_stability_ensemble.R

diff --git a/scripts/plotting/mcsm_mean_stability_ensemble.R b/scripts/plotting/mcsm_mean_stability_ensemble.R
new file mode 100644
index 0000000..9bc1be7
--- /dev/null
+++ b/scripts/plotting/mcsm_mean_stability_ensemble.R
@@ -0,0 +1,280 @@
+source("~/git/LSHTM_analysis/config/pnca.R")
+
+source("/home/tanu/git/LSHTM_analysis/my_header.R")
+#########################################################
+# TASK: Generate averaged stability values 
+# across all stability tools
+# for a given structure
+#########################################################
+
+#=======
+# output
+#=======
+outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
+
+outfile_mean_ens_st_aff = paste0(outdir_images, "/", tolower(gene)
+                                , "_mean_ens_stab_aff.csv")
+
+print(paste0("Output file:", outfile_mean_ens_st_aff))
+
+#%%===============================================================
+
+#=============
+# Input
+#=============
+df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv")
+df3 = read.csv(df3_filename)
+
+# mut_info checks
+table(df3$mutation_info)
+table(df3$mutation_info_orig)
+table(df3$mutation_info_labels_orig)
+
+# used in plots and analyses
+table(df3$mutation_info_labels) # different, and matches dst_mode
+table(df3$dst_mode)
+
+# create column based on dst mode with different colname
+table(is.na(df3$dst))
+table(is.na(df3$dst_mode))
+
+#===============
+# Create column: sensitivity mapped to dst_mode
+#===============
+df3$sensitivity = ifelse(df3$dst_mode == 1, "R", "S")
+table(df3$sensitivity)
+
+length(unique((df3$mutationinformation)))
+all_colnames = as.data.frame(colnames(df3))
+common_cols  = c("mutationinformation"
+                 , "position"
+                 , "dst_mode"
+                 #, "mutation_info_labels"
+                 , "sensitivity"
+                 , "ligand_distance")
+
+# ADD the ones for mcsm_na etc
+#optional_cols = c() 
+
+all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)]
+#TODO: affinity_cols
+scaled_cols = c("duet_scaled"          , "duet_stability_change"
+                ,"deepddg_scaled"      , "deepddg"   
+                ,"ddg_dynamut2_scaled" , "ddg_dynamut2"
+                ,"foldx_scaled"        , "ddg_foldx"
+                
+                #,"consurf_scaled"      , "consurf_score"
+                #,"snap2_scaled"        , "snap2_score"
+                #,"provean_scaled"      , "provean_score"
+                #,"affinity_scaled"     , "ligand_affinity_change"
+                #,"mmcsm_lig_scaled"    , "mmcsm_lig"
+                )
+all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)]
+outcome_cols = c("duet_outcome"
+                 , "deepddg_outcome"
+                 , "ddg_dynamut2_outcome"
+                 , "foldx_outcome"
+                 #, "ddg_foldx", "foldx_scaled"
+                 
+                 # consurf outcome doesn't exist
+                 #,"provean_outcome"
+                 #,"snap2_outcome"
+                 #,"ligand_outcome"
+                 #,"mmcsm_lig_outcome"
+                 )
+outcome_cols_affinity =  c("ligand_outcome"
+                           ,"mmcsm_lig_outcome")
+
+cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols, scaled_cols, outcome_cols, outcome_cols_affinity)]
+cols_to_extract = cols_to_consider[cols_to_consider%in%c(common_cols, outcome_cols)]
+
+foo = df3[, cols_to_consider]
+df3_plot_orig = df3[, cols_to_extract]
+
+##############################################################
+#####################
+# Ensemble stability
+#####################
+# extract outcome cols and map numeric values to the categories
+# Destabilising == 1, and stabilising == 0
+df3_plot = df3[, cols_to_extract]
+
+df3_plot[, outcome_cols] <- sapply(df3_plot[, outcome_cols]
+                             , function(x){ifelse(x == "Destabilising", 1, 0)})
+
+#=====================================
+# Stability (4 cols): average the scores
+# across predictors ==> average by
+# position ==> scale b/w -1 and 1
+
+# column to average: ens_stability
+#=====================================
+cols_to_average = which(colnames(df3_plot)%in%outcome_cols)
+
+# ensemble average across predictors
+df3_plot$ens_stability = rowMeans(df3_plot[,cols_to_average])
+
+head(df3_plot$position); head(df3_plot$mutationinformation)
+head(df3_plot$ens_stability)
+table(df3_plot$ens_stability)
+
+# ensemble average of predictors by position
+mean_ens_stability_by_position <- df3_plot %>%
+  dplyr::group_by(position) %>%
+  dplyr::summarize(avg_ens_stability = mean(ens_stability))
+
+# REscale b/w -1 and 1
+#en_stab_min = min(mean_ens_stability_by_position['avg_ens_stability'])
+#en_stab_max = max(mean_ens_stability_by_position['avg_ens_stability']) 
+
+# scale the average stability value between -1 and 1
+# mean_ens_by_position['averaged_stability3_scaled'] = lapply(mean_ens_by_position['averaged_stability3']
+#                                                        , function(x) ifelse(x < 0, x/abs(en3_min), x/en3_max))
+
+mean_ens_stability_by_position['avg_ens_stability_scaled'] = lapply(mean_ens_stability_by_position['avg_ens_stability']
+                                                                  , function(x) {
+                                                                    scales::rescale(x, to  = c(-1,1)
+                                                                                    #, from = c(en_stab_min,en_stab_max))
+                                                                                    , from = c(0,1))
+                                                                  })
+cat(paste0('Average stability scores:\n'
+           , head(mean_ens_stability_by_position['avg_ens_stability'])
+           , '\n---------------------------------------------------------------'
+           , '\nAverage stability scaled scores:\n'
+           , head(mean_ens_stability_by_position['avg_ens_stability_scaled'])))
+
+# convert to a data frame
+mean_ens_stability_by_position = as.data.frame(mean_ens_stability_by_position)
+
+#FIXME: sanity checks
+# TODO: predetermine the bounds
+# l_bound_ens = min(mean_ens_stability_by_position['avg_ens_stability_scaled'])
+# u_bound_ens = max(mean_ens_stability_by_position['avg_ens_stability_scaled'])
+# 
+# if ( (l_bound_ens == -1) && (u_bound_ens == 1) ){
+#   cat(paste0("PASS: ensemble stability scores averaged by position and then scaled"
+#         , "\nmin ensemble averaged stability: ", l_bound_ens
+#         , "\nmax ensemble averaged stability: ", u_bound_ens))
+# }else{
+#   cat(paste0("FAIL: avergaed duet scores could not be scaled b/w -1 and 1"
+#         , "\nmin ensemble averaged stability: ", l_bound_ens
+#         , "\nmax ensemble averaged stability: ", u_bound_ens))
+#   quit()
+# } 
+##################################################################
+############################
+# Ensemble affinity: ligand
+############################
+# extract ligand affinity outcome cols and map numeric values to the categories
+# Destabilising == 1, and stabilising == 0
+cols_to_extract_affinity = cols_to_consider[cols_to_consider%in%c(common_cols
+                                                                  , outcome_cols_affinity)]
+
+
+df3_plot_affinity = df3[, cols_to_extract_affinity]
+names(df3_plot_affinity)
+
+df3_plot_affinity[, outcome_cols_affinity] <- sapply(df3_plot_affinity[, outcome_cols_affinity]
+                                                     , function(x){ifelse(x == "Destabilising", 1, 0)})
+
+#=====================================
+# Affintiy (2 cols): average the scores
+# across predictors ==> average by
+# position ==> scale b/w -1 and 1
+
+# column to average: ens_affinity
+#=====================================
+cols_to_average_affinity = which(colnames(df3_plot_affinity)%in%outcome_cols_affinity)
+cols_to_average_affinity
+
+# ensemble average across predictors
+df3_plot_affinity$ens_affinity = rowMeans(df3_plot_affinity[,cols_to_average_affinity])
+
+head(df3_plot_affinity$position); head(df3_plot_affinity$mutationinformation)
+head(df3_plot_affinity$ens_affinity)
+table(df3_plot_affinity$ens_affinity)
+
+# ensemble average of predictors by position
+mean_ens_affinity_by_position <- df3_plot_affinity %>%
+  dplyr::group_by(position) %>%
+  dplyr::summarize(avg_ens_affinity = mean(ens_affinity))
+
+# REscale b/w -1 and 1
+#en_aff_min = min(mean_ens_affinity_by_position['ens_affinity'])
+#en_aff_max = max(mean_ens_affinity_by_position['ens_affinity']) 
+
+# scale the average affintiy value between -1 and 1
+# mean_ens_affinity_by_position['avg_ens_affinity_scaled'] = lapply(mean_ens_affinity_by_position['avg_ens_affinity']
+#                                                        , function(x) ifelse(x < 0, x/abs(en_aff_min), x/en_aff_max))
+
+
+mean_ens_affinity_by_position['avg_ens_affinity_scaled'] = lapply(mean_ens_affinity_by_position['avg_ens_affinity']
+                                                                  , function(x) {
+                                                                    scales::rescale(x, to  = c(-1,1)
+                                                                                    #, from = c(en_aff_min,en_aff_max))
+                                                                                    , from = c(0,1))
+                                                                  })
+cat(paste0('Average affintiy scores:\n'
+           , head(mean_ens_affinity_by_position['avg_ens_affinity'])
+           , '\n---------------------------------------------------------------'
+           , '\nAverage affintiy scaled scores:\n'
+           , head(mean_ens_affinity_by_position['avg_ens_affinity_scaled'])))
+
+#convert to a df
+mean_ens_affinity_by_position = as.data.frame(mean_ens_affinity_by_position)
+
+
+#FIXME: sanity checks
+# TODO: predetermine the bounds
+# l_bound_ens_aff = min(mean_ens_affintiy_by_position['avg_ens_affinity_scaled'])
+# u_bound_ens_aff = max(mean_ens_affintiy_by_position['avg_ens_affinity_scaled'])
+# 
+# if ( (l_bound_ens_aff == -1) && (u_bound_ens_aff == 1) ){
+#   cat(paste0("PASS: ensemble affinity scores averaged by position and then scaled"
+#         , "\nmin ensemble averaged affinity: ", l_bound_ens_aff
+#         , "\nmax ensemble averaged affinity: ", u_bound_ens_aff))
+# }else{
+#   cat(paste0("FAIL: ensemble affinity scores could not be scaled b/w -1 and 1"
+#         , "\nmin ensemble averaged affinity: ", l_bound_ens_aff
+#         , "\nmax ensemble averaged affinity: ", u_bound_ens_aff))
+#   quit()
+# } 
+
+
+######################################################################
+##################
+# merge: mean ensemble stability and affinity by_position
+####################
+# if ( class(mean_ens_stability_by_position) && class(mean_ens_affinity_by_position) != "data.frame"){
+#   cat("Y")
+# }
+
+common_cols = intersect(colnames(mean_ens_stability_by_position), colnames(mean_ens_affinity_by_position))
+
+if (dim(mean_ens_stability_by_position) && dim(mean_ens_affinity_by_position)){
+  print(paste0("PASS: dim's match, mering dfs by column :", common_cols))
+  #combined = as.data.frame(cbind(mean_duet_by_position, mean_affinity_by_position ))
+  combined_df = as.data.frame(merge(mean_ens_stability_by_position
+                                    , mean_ens_affinity_by_position
+                                    , by = common_cols
+                                    , all = T))
+  
+  cat(paste0("\nnrows combined_df:", nrow(combined_df)
+               , "\nnrows combined_df:", ncol(combined_df)))
+}else{
+    cat(paste0("FAIL: dim's mismatch, aborting cbind!"
+          , "\nnrows df1:", nrow(mean_duet_by_position)
+          , "\nnrows df2:", nrow(mean_affinity_by_position)))
+    quit()      
+}
+#%%============================================================
+# output
+write.csv(combined_df, outfile_mean_ens_st_aff
+          , row.names = F)
+cat("Finished writing file:\n"
+    , outfile_mean_ens_st_aff
+    , "\nNo. of rows:", nrow(combined_df)
+    , "\nNo. of cols:", ncol(combined_df))
+
+# end of script
+#===============================================================