From 1bf66b145c5ff40022214fc34c270a81f0cdc05d Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Sun, 31 Jul 2022 19:24:35 +0100 Subject: [PATCH] separating mcsm_mean_stability_ensemble from combined script --- .../plotting/mcsm_mean_stability_ensemble.R | 147 +++--------------- scripts/plotting/replaceBfactor_pdb.R | 132 +++++++++------- 2 files changed, 99 insertions(+), 180 deletions(-) diff --git a/scripts/plotting/mcsm_mean_stability_ensemble.R b/scripts/plotting/mcsm_mean_stability_ensemble.R index 9bc1be7..2f5abb2 100644 --- a/scripts/plotting/mcsm_mean_stability_ensemble.R +++ b/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_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] +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)] ############################################################## ##################### # 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 cols:", ncol(combined_df)) + , "\nNo. of rows:", nrow(mean_ens_stability_by_position) + , "\nNo. of cols:", ncol(mean_ens_stability_by_position)) # end of script #=============================================================== diff --git a/scripts/plotting/replaceBfactor_pdb.R b/scripts/plotting/replaceBfactor_pdb.R index 0d5cd1a..a840c19 100755 --- a/scripts/plotting/replaceBfactor_pdb.R +++ b/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") -infile_mean_stability = paste0(outdir, "/", in_filename_mean_stability) +#in_filename_mean_stability = paste0(tolower(gene), "_mean_stability.csv") +#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 -df_lig$b[is.na(df_lig$b)] = 0 +na_rep = 2 +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 -# duet and lig -if ( (sum(df_duet$b == 0) == b_na_duet) && (sum(df_lig$b == 0) == b_na_lig) ) { - print ("PASS: NA's replaced with 0s successfully in df_duet and df_lig") -} else { - print("FAIL: NA replacement in df_duet NOT successful") - quit() -} +# # sanity check: should be 0 and True +# # duet and 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") +# } else { +# print("FAIL: NA replacement in df_duet NOT successful") +# 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) - -# sanity checks: should be True -if( (max(df_duet$b) == max(my_df$averaged_duet_scaled)) & (min(df_duet$b) == min(my_df$averaged_duet_scaled)) ){ - print("PASS: B-factors replaced correctly in df_duet") -} 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() -} +# 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") +# } 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??? #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!