separating mcsm_mean_stability_ensemble from combined script

scripts/plotting/mcsm_mean_stability_ensemble.R

@@ -1,4 +1,9 @@
 source("~/git/LSHTM_analysis/config/pnca.R")
+#source("~/git/LSHTM_analysis/config/alr.R")
+#source("~/git/LSHTM_analysis/config/gid.R")
+#source("~/git/LSHTM_analysis/config/embb.R")
+#source("~/git/LSHTM_analysis/config/katg.R")
+#source("~/git/LSHTM_analysis/config/rpob.R")
 
 source("/home/tanu/git/LSHTM_analysis/my_header.R")
 #########################################################
@@ -11,10 +16,8 @@ source("/home/tanu/git/LSHTM_analysis/my_header.R")
 # 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")
+                                , "_mean_ens_stability.csv")
-
 print(paste0("Output file:", outfile_mean_ens_st_aff))
 
 #%%===============================================================
@@ -24,6 +27,7 @@ print(paste0("Output file:", outfile_mean_ens_st_aff))
 #=============
 df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv")
 df3 = read.csv(df3_filename)
+length(df3$mutationinformation)
 
 # mut_info checks
 table(df3$mutation_info)
@@ -57,12 +61,12 @@ common_cols  = c("mutationinformation"
 #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"
-                
+                , "mcsm_ppi2_scaled"   , "mcsm_ppi2_affinity"
+                , "mcsm_na_scaled"     , "mcsm_na_affinity"
                 #,"consurf_scaled"      , "consurf_score"
                 #,"snap2_scaled"        , "snap2_score"
                 #,"provean_scaled"      , "provean_score"
@@ -81,27 +85,26 @@ outcome_cols = c("duet_outcome"
                  #,"snap2_outcome"
                  #,"ligand_outcome"
                  #,"mmcsm_lig_outcome"
+                 #, "mcsm_ppi2_outcome"
+                 #, "mcsm_na_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_consider = colnames(df3)[colnames(df3)%in%c(common_cols, scaled_cols,outcome_cols)]
-cols_to_extract = cols_to_consider[cols_to_consider%in%c(common_cols, outcome_cols)]
+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
+# Destabilising == 0, and stabilising == 1, so rescaling can let -1 be destabilising
 df3_plot = df3[, cols_to_extract]
 
+# assign numeric values to outcome
 df3_plot[, outcome_cols] <- sapply(df3_plot[, outcome_cols]
-                             , function(x){ifelse(x == "Destabilising", 1, 0)})
+                             , function(x){ifelse(x == "Destabilising", 0, 1)})
-
+table(df3$duet_outcome)
+table(df3_plot$duet_outcome)
 #=====================================
 # Stability (4 cols): average the scores
 # across predictors ==> average by
@@ -162,119 +165,15 @@ mean_ens_stability_by_position = as.data.frame(mean_ens_stability_by_position)
 #   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
+#write.csv(combined_df, outfile_mean_ens_st_aff
+write.csv(mean_ens_stability_by_position
+          , 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 rows:", nrow(mean_ens_stability_by_position)
-    , "\nNo. of cols:", ncol(combined_df))
+    , "\nNo. of cols:", ncol(mean_ens_stability_by_position))
 
 # end of script
 #===============================================================

scripts/plotting/replaceBfactor_pdb.R

@@ -10,7 +10,7 @@
 # rendering on chimera
 
 # read mcsm mean stability value files
-# extract the respecitve mean values and assign to the
+# extract the respective mean values and assign to the
 # b-factor column within their respective pdbs
 
 # generate some distribution plots for inspection
@@ -52,7 +52,8 @@ cat(gene_match)
 datadir = paste0("~/git/Data")
 indir = paste0(datadir, "/", drug, "/input")
 outdir = paste0("~/git/Data", "/", drug, "/output")
-outdir_plots = paste0("~/git/Data", "/", drug, "/output/plots")
+#outdir_plots = paste0("~/git/Data", "/", drug, "/output/plots")
+outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
 
 #======
 # input
@@ -61,14 +62,19 @@ in_filename_pdb = paste0(tolower(gene), "_complex.pdb")
 infile_pdb = paste0(indir, "/", in_filename_pdb)
 cat(paste0("Input file:", infile_pdb) )
 
-in_filename_mean_stability = paste0(tolower(gene), "_mean_stability.csv") 
+#in_filename_mean_stability = paste0(tolower(gene), "_mean_stability.csv")
-infile_mean_stability = paste0(outdir, "/", in_filename_mean_stability)
+#infile_mean_stability = paste0(outdir, "/", in_filename_mean_stability)
+
+in_filename_mean_stability = paste0(tolower(gene), "_mean_ens_stab_aff.csv")
+infile_mean_stability = paste0(outdir_plots, "/", in_filename_mean_stability)
+
 cat(paste0("Input file:", infile_mean_stability) )
 
 #=======
 # output
 #=======
-out_filename_duet_mspdb = paste0(tolower(gene), "_complex_bduet_ms.pdb") 
+#out_filename_duet_mspdb = paste0(tolower(gene), "_complex_bduet_ms.pdb")
+out_filename_duet_mspdb = paste0(tolower(gene), "_complex_b_stab_ms.pdb") 
 outfile_duet_mspdb = paste0(outdir_plots, "/", out_filename_duet_mspdb)
 print(paste0("Output file:", outfile_duet_mspdb))
 
@@ -77,6 +83,8 @@ outfile_lig_mspdb = paste0(outdir_plots, "/", out_filename_lig_mspdb)
 print(paste0("Output file:", outfile_lig_mspdb))
 
 #%%===============================================================
+#NOTE: duet here refers to the ensemble stability values
+
 ###########################
 # Read file: average stability values
 # or mcsm_normalised file
@@ -133,17 +141,17 @@ par(oma = c(3,2,3,0)
     #, mfrow = c(3,2)
     , mfrow = c(3,4))
 
-#************
+#=============
 # Row 1 plots: original B-factors
 # duet and affinity
-#************
+#=============
 hist(df_duet$b
      , xlab = "" 
-     , main = "Bfactor duet")
+     , main = "Bfactor stability")
 
 plot(density(df_duet$b)
      , xlab = ""
-     , main = "Bfactor duet")
+     , main = "Bfactor stability")
 
 
 hist(df_lig$b
@@ -154,32 +162,36 @@ plot(density(df_lig$b)
      , xlab = ""
      , main = "Bfactor affinity")
 
-#************
+#=============
 # Row 2 plots: original mean stability values
 # duet and affinity
-#************
+#=============
-hist(my_df$averaged_duet
+
+#hist(my_df$averaged_duet
+hist(my_df$avg_ens_stability_scaled
      , xlab = "" 
-     , main = "mean duet values")
+     , main = "mean stability values")
 
-plot(density(my_df$averaged_duet)
+#plot(density(my_df$averaged_duet)
+plot(density(my_df$avg_ens_stability_scaled)
      , xlab = ""
-     , main = "mean duet values")
+     , main = "mean stability values")
 
-
+#hist(my_df$averaged_affinity
-hist(my_df$averaged_affinity
+hist(my_df$avg_ens_affinity_scaled
      , xlab = "" 
      , main = "mean affinity values")
 
-plot(density(my_df$averaged_affinity)
+#plot(density(my_df$averaged_affinity)
+plot(density(my_df$avg_ens_affinity_scaled)
      , xlab = ""
      , main = "mean affinity values")
 
-#************
+#==============
 # Row 3 plots: replaced B-factors with mean stability values
 # After actual replacement in the b factor column
-#*************
+#===============
-#=========================================================
+################################################################
 #=========
 # step 0_P1: DONT RUN once you have double checked the matched output
 #=========
@@ -192,49 +204,54 @@ plot(density(my_df$averaged_affinity)
 #=========
 # Be brave and replace in place now (don"t run sanity check)
 # this makes all the B-factor values in the non-matched positions as NA
-df_duet$b = my_df$averaged_duet_scaled[match(df_duet$resno, my_df$position)]
+
-df_lig$b = my_df$averaged_affinity_scaled[match(df_lig$resno, my_df$position)]
+#df_duet$b = my_df$averaged_duet_scaled[match(df_duet$resno, my_df$position)]
+#df_lig$b = my_df$averaged_affinity_scaled[match(df_lig$resno, my_df$position)]
+
+df_duet$b = my_df$avg_ens_stability_scaled[match(df_duet$resno, my_df$position)]
+df_lig$b  = my_df$avg_ens_affinity_scaled[match(df_lig$resno, my_df$position)]
 
 #=========
 # step 2_P1
 #=========
 # count NA in Bfactor
 b_na_duet = sum(is.na(df_duet$b)) ; b_na_duet
-b_na_lig = sum(is.na(df_lig$b)) ; b_na_lig 
+b_na_lig  = sum(is.na(df_lig$b)) ; b_na_lig 
 
 # count number of 0"s in Bactor
 sum(df_duet$b == 0)
-sum(df_lig$b == 0)
+sum(df_lig$b  == 0)
 
 # replace all NA in b factor with 0
-df_duet$b[is.na(df_duet$b)] = 0
+na_rep = 2
-df_lig$b[is.na(df_lig$b)] = 0
+df_duet$b[is.na(df_duet$b)] = na_rep
+df_lig$b[is.na(df_lig$b)] = na_rep
 
-# sanity check: should be 0 and True
+# # sanity check: should be 0 and True
-# duet and lig
+# # duet and lig
-if ( (sum(df_duet$b == 0) == b_na_duet) && (sum(df_lig$b == 0) == b_na_lig) ) {
+# if ( (sum(df_duet$b == na_rep) == b_na_duet) && (sum(df_lig$b == na_rep) == b_na_lig) ) {
-  print ("PASS: NA's replaced with 0s successfully in df_duet and df_lig")
+#   print ("PASS: NA's replaced with 0s successfully in df_duet and df_lig")
-} else {
+# } else {
-  print("FAIL: NA replacement in df_duet NOT successful")
+#   print("FAIL: NA replacement in df_duet NOT successful")
-  quit()
+#   quit()
-}
+# }
+# 
+# max(df_duet$b); min(df_duet$b)
+# 
+# # sanity checks: should be True
+# if( (max(df_duet$b) == max(my_df$avg_ens_stability_scaled)) & (min(df_duet$b) == min(my_df$avg_ens_stability_scaled)) ){
+#   print("PASS: B-factors replaced correctly in df_duet")
+# } else {
+#   print ("FAIL: To replace B-factors in df_duet")
+#   quit()
+# }
 
-max(df_duet$b); min(df_duet$b)
+# if( (max(df_lig$b) == max(my_df$avg_ens_affinity_scaled)) & (min(df_lig$b) == min(my_df$avg_ens_affinity_scaled)) ){
-
+#   print("PASS: B-factors replaced correctly in df_lig")
-# sanity checks: should be True
+# } else {
-if( (max(df_duet$b) == max(my_df$averaged_duet_scaled)) & (min(df_duet$b) == min(my_df$averaged_duet_scaled)) ){
+#   print ("FAIL: To replace B-factors in df_lig")
-  print("PASS: B-factors replaced correctly in df_duet")
+#   quit()
-} else {
+# }
-  print ("FAIL: To replace B-factors in df_duet")
-  quit()
-}
-
-if( (max(df_lig$b) == max(my_df$averaged_affinity_scaled)) & (min(df_lig$b) == min(my_df$averaged_affinity_scaled)) ){
-  print("PASS: B-factors replaced correctly in df_lig")
-} else {
-  print ("FAIL: To replace B-factors in df_lig")
-  quit()
-}
 
 #=========
 # step 3_P1
@@ -255,6 +272,8 @@ if ( (dim(df_duet)[1] == dim(d2_duet)[1]) & (dim(df_lig)[1] == dim(d2_lig)[1]) &
 # assign it back to the pdb file
 my_pdb_duet[['atom']] = df_duet
 max(df_duet$b); min(df_duet$b)
+table(df_duet$b)
+sum(is.na(df_duet$b))
 
 my_pdb_lig[['atom']] = df_lig
 max(df_lig$b); min(df_lig$b)
@@ -268,9 +287,9 @@ write.pdb(my_pdb_duet, outfile_duet_mspdb)
 cat(paste0("output file ligand mean stability pdb:", outfile_lig_mspdb))
 write.pdb(my_pdb_lig, outfile_lig_mspdb)
 
-#********************************
+#============================
 # Add the 3rd histogram and density plots for comparisons
-#********************************
+#============================
 # Plots continued...
 # Row 3 plots: hist and density of replaced B-factors with stability values
 hist(df_duet$b
@@ -296,16 +315,17 @@ mtext(text = "Frequency"
       , line = 0
       , outer = TRUE)
 
-mtext(text = "Stability Distribution"
+mtext(text = paste0(tolower(gene), ": Stability Distribution")
       , side = 3
       , line = 0
       , outer = TRUE)
-#********************************
+#============================================
 
 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 # NOTE: This replaced B-factor distribution has the same
 # x-axis as the PredAff normalised values, but the distribution
-# is affected since 0 is overinflated. This is because all the positions
+# is affected since 0 is overinflated/or hs an additional blip because
+# of the positions not associated with resistance. This is because all the positions
 # where there are no SNPs have been assigned 0???
 #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!