From aabe46659935974cf0352b06ccdc15c1254e3065 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Wed, 3 Aug 2022 21:32:47 +0100 Subject: [PATCH] added plots for thesis --- scripts/Header_TT.R | 13 + scripts/functions/lineage_plot_data.R | 2 + scripts/functions/position_count_bp.R | 5 +- scripts/functions/stability_count_bp.R | 9 +- scripts/plotting/LINEAGE2.R | 3 - scripts/plotting/mcsm_affinity_data_only.R | 241 ------------- .../plotting/mcsm_mean_affinity_ensemble.R | 316 ----------------- scripts/plotting/mcsm_mean_stability.R | 163 --------- .../plotting/mcsm_mean_stability_ensemble.R | 212 ----------- .../mcsm_mean_stability_ensemble_5uhc_rpob.R | 176 ---------- scripts/plotting/mut_landscape.R | 155 -------- scripts/plotting/mut_landscape_5uhc_rpob.R | 191 ---------- scripts/plotting/replaceBfactor_pdb.R | 332 ------------------ .../plotting/replaceBfactor_pdb_affinity.R | 281 --------------- scripts/plotting/replaceBfactor_pdb_ppi2.R | 277 --------------- .../plotting/replaceBfactor_pdb_stability.R | 281 --------------- .../replaceBfactor_pdb_stability_5uhc_rpob.R | 280 --------------- 17 files changed, 24 insertions(+), 2913 deletions(-) delete mode 100644 scripts/plotting/mcsm_affinity_data_only.R delete mode 100644 scripts/plotting/mcsm_mean_affinity_ensemble.R delete mode 100755 scripts/plotting/mcsm_mean_stability.R delete mode 100644 scripts/plotting/mcsm_mean_stability_ensemble.R delete mode 100644 scripts/plotting/mcsm_mean_stability_ensemble_5uhc_rpob.R delete mode 100644 scripts/plotting/mut_landscape.R delete mode 100644 scripts/plotting/mut_landscape_5uhc_rpob.R delete mode 100755 scripts/plotting/replaceBfactor_pdb.R delete mode 100644 scripts/plotting/replaceBfactor_pdb_affinity.R delete mode 100644 scripts/plotting/replaceBfactor_pdb_ppi2.R delete mode 100644 scripts/plotting/replaceBfactor_pdb_stability.R delete mode 100644 scripts/plotting/replaceBfactor_pdb_stability_5uhc_rpob.R diff --git a/scripts/Header_TT.R b/scripts/Header_TT.R index e336ed2..e500edb 100755 --- a/scripts/Header_TT.R +++ b/scripts/Header_TT.R @@ -222,6 +222,19 @@ consurf_palette2 = c("0" = "yellow2" , "8" = "orchid4" , "9" = "darkorchid4") + +consurf_colours = c("0" = rgb(1.00,1.00,0.59) + , "1" = rgb(0.63,0.16,0.37) + , "2" = rgb(0.94,0.49,0.67) + , "3" = rgb(0.98,0.78,0.86) + , "4" = rgb(0.98,0.92,0.96) + , "5" = rgb(1.00,1.00,1.00) + , "6" = rgb(0.84,0.94,0.94) + , "7" = rgb(0.65,0.86,0.90) + , "8" = rgb(0.29,0.69,0.75) + , "9" = rgb(0.04,0.49,0.51) + ) + ################################################## # Function name clashes with plyr and dplyr diff --git a/scripts/functions/lineage_plot_data.R b/scripts/functions/lineage_plot_data.R index 116f156..3aa75a5 100644 --- a/scripts/functions/lineage_plot_data.R +++ b/scripts/functions/lineage_plot_data.R @@ -34,6 +34,8 @@ lineage_plot_data <- function(df ################################################################ # Get WF and LF data with lineage count, and snp diversity ################################################################ + + df[lineage_column_name] = # Initialise output list lineage_dataL = list( lin_wf = data.frame() diff --git a/scripts/functions/position_count_bp.R b/scripts/functions/position_count_bp.R index c08f2d5..876a0fc 100755 --- a/scripts/functions/position_count_bp.R +++ b/scripts/functions/position_count_bp.R @@ -96,7 +96,7 @@ site_snp_count_bp <- function (plotdf # but atm being using as plot title #my_leg_title bp_plot_title = paste0("Total nsSNPs: ", tot_muts - , ", Total no. of nsSNPs sites: ", tot_sites) + , "\nTotal sites: ", tot_sites) #------------- # start plot 2 @@ -123,7 +123,8 @@ site_snp_count_bp <- function (plotdf #, legend.text = element_text(size = leg_text_size) #, legend.title = element_text(size = axis_label_size) , plot.title = element_text(size = leg_text_size - , colour = title_colour) + , colour = title_colour + , hjust = 0.5) , plot.subtitle = element_text(size = subtitle_size , hjust = 0.5 , colour = subtitle_colour)) + diff --git a/scripts/functions/stability_count_bp.R b/scripts/functions/stability_count_bp.R index 1a38059..136cd67 100644 --- a/scripts/functions/stability_count_bp.R +++ b/scripts/functions/stability_count_bp.R @@ -30,7 +30,8 @@ stability_count_bp <- function(plotdf , sts = 20 , subtitle_colour = "pink" #, leg_position = c(0.73,0.8) # within plot area - , leg_position = "top"){ + , leg_position = "top" + , bar_fill_values = c("#F8766D", "#00BFC4")){ # OutPlot_count = ggplot(plotdf, aes(x = eval(parse(text = df_colname)))) + OutPlot_count = ggplot(plotdf, aes_string(x = df_colname)) + @@ -57,8 +58,10 @@ stability_count_bp <- function(plotdf , subtitle = subtitle_text , y = yaxis_title) + scale_fill_discrete(name = leg_title - #, labels = c("Destabilising", "Stabilising") - , labels = label_categories) + , labels = label_categories) + + + scale_fill_manual("", values=bar_fill_values) + return(OutPlot_count) } diff --git a/scripts/plotting/LINEAGE2.R b/scripts/plotting/LINEAGE2.R index e680a2f..716f2fc 100644 --- a/scripts/plotting/LINEAGE2.R +++ b/scripts/plotting/LINEAGE2.R @@ -128,9 +128,6 @@ plot_df$pvalRF = ifelse(plot_df$pvalR > 0.05, paste0("p=",plot_df$pvalR), plot_d # plot_df$pvalF = ifelse(plot_df$pval < 0.05, paste0(round(plot_df$pval, 3), "* "), plot_df$pval ) # plot_df$pvalF = ifelse(plot_df$pval == 0.05, paste0(round(plot_df$pval, 3), ". "), plot_df$pval ) -round(plot_df$pvalF, 3) - - #================================================ # Plot attempt 1 [no stats]: WORKS beeautifully #================================================ diff --git a/scripts/plotting/mcsm_affinity_data_only.R b/scripts/plotting/mcsm_affinity_data_only.R deleted file mode 100644 index f962935..0000000 --- a/scripts/plotting/mcsm_affinity_data_only.R +++ /dev/null @@ -1,241 +0,0 @@ -#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") -######################################################### -# TASK: Generate averaged affinity values -# across all affinity tools for a given structure -# as applicable... -######################################################### - -#======= -# output -#======= -outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) - -#OutFile1 -outfile_mean_aff = paste0(outdir_images, "/", tolower(gene) - , "_mean_affinity_all.csv") -print(paste0("Output file:", outfile_mean_aff)) - -#OutFile2 -outfile_mean_aff_priorty = paste0(outdir_images, "/", tolower(gene) - , "_mean_affinity_priority.csv") -print(paste0("Output file:", outfile_mean_aff_priorty)) - -#%%=============================================================== - -#============= -# Input -#============= -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) -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)) - -# FIXME: ADD distance to NA when SP replies -dist_columns = c("ligand_distance", "interface_dist") -DistCutOff = 10 -common_cols = c("mutationinformation" - , "X5uhc_position" - , "X5uhc_offset" - , "position" - , "dst_mode" - , "mutation_info_labels" - , "sensitivity", dist_columns ) - -all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)] -all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)] - -#=================== -# stability cols -#=================== -raw_cols_stability = c("duet_stability_change" - , "deepddg" - , "ddg_dynamut2" - , "ddg_foldx") - -scaled_cols_stability = c("duet_scaled" - , "deepddg_scaled" - , "ddg_dynamut2_scaled" - , "foldx_scaled") - -outcome_cols_stability = c("duet_outcome" - , "deepddg_outcome" - , "ddg_dynamut2_outcome" - , "foldx_outcome") - -#=================== -# affinity cols -#=================== -raw_cols_affinity = c("ligand_affinity_change" - , "mmcsm_lig" - , "mcsm_ppi2_affinity" - , "mcsm_na_affinity") - -scaled_cols_affinity = c("affinity_scaled" - , "mmcsm_lig_scaled" - , "mcsm_ppi2_scaled" - , "mcsm_na_scaled" ) - -outcome_cols_affinity = c( "ligand_outcome" - , "mmcsm_lig_outcome" - , "mcsm_ppi2_outcome" - , "mcsm_na_outcome") - -#=================== -# conservation cols -#=================== -# raw_cols_conservation = c("consurf_score" -# , "snap2_score" -# , "provean_score") -# -# scaled_cols_conservation = c("consurf_scaled" -# , "snap2_scaled" -# , "provean_scaled") -# -# # CANNOT strictly be used, as categories are not identical with conssurf missing altogether -# outcome_cols_conservation = c("provean_outcome" -# , "snap2_outcome" -# #consurf outcome doesn't exist -# ) - -gene_aff_cols = colnames(df3)[colnames(df3)%in%scaled_cols_affinity] -gene_stab_cols = colnames(df3)[colnames(df3)%in%scaled_cols_stability] -gene_common_cols = colnames(df3)[colnames(df3)%in%common_cols] - -sel_cols = c(gene_common_cols - , gene_stab_cols - , gene_aff_cols) - -######################################### -#df3_plot = df3[, cols_to_extract] -df3_plot = df3[, sel_cols] - -###################### -#FILTERING HMMMM! -#all dist <10 -#for embb this results in 2 muts -###################### -df3_affinity_filtered = df3_plot[ (df3_plot$ligand_distance<10 | df3_plot$interface_dist <10),] -df3_affinity_filtered = df3_plot[ (df3_plot$ligand_distance<10 & df3_plot$interface_dist <10),] - -c0u = unique(df3_affinity_filtered$position) -length(c0u) - -#df = df3_affinity_filtered -########################################## -#NO FILTERING: annotate the whole df -df = df3_plot -sum(is.na(df)) -df2 = na.omit(df) -c0u = unique(df2$position) -length(c0u) - -# reassign orig -my_df_raw = df3 - -# now subset -df3 = df2 -####################################################### -#================= -# affinity effect -#================= -give_col=function(x,y,df=df3){ - df[df$position==x,y] -} - -for (i in unique(df3$position) ){ - #print(i) - biggest = max(abs(give_col(i,gene_aff_cols))) - - df3[df3$position==i,'abs_max_effect'] = biggest - df3[df3$position==i,'effect_type']= names( - give_col(i,gene_aff_cols)[which( - abs( - give_col(i,gene_aff_cols) - ) == biggest, arr.ind=T - )[, "col"]]) - - # effect_name = unique(df3[df3$position==i,'effect_type']) - effect_name = df3[df3$position==i,'effect_type'][1] # pick first one in case we have multiple exact values - - ind = rownames(which(abs(df3[df3$position==i,c('position',effect_name)][effect_name])== biggest, arr.ind=T)) - df3[df3$position==i,'effect_sign'] = sign(df3[effect_name][ind,]) -} - -df3$effect_type = sub("\\.[0-9]+", "", df3$effect_type) # cull duplicate effect types that happen when there are exact duplicate values -df3U = df3[!duplicated(df3[c('position')]), ] -table(df3U$effect_type) -######################################################### -#%% consider stability as well -df4 = df2 - -#================= -# stability + affinity effect -#================= -effect_cols = c(gene_aff_cols, gene_stab_cols) - -give_col=function(x,y,df=df4){ - df[df$position==x,y] -} - -for (i in unique(df4$position) ){ - #print(i) - biggest = max(abs(give_col(i,effect_cols))) - - df4[df4$position==i,'abs_max_effect'] = biggest - df4[df4$position==i,'effect_type']= names( - give_col(i,effect_cols)[which( - abs( - give_col(i,effect_cols) - ) == biggest, arr.ind=T - )[, "col"]]) - - # effect_name = unique(df4[df4$position==i,'effect_type']) - effect_name = df4[df4$position==i,'effect_type'][1] # pick first one in case we have multiple exact values - - ind = rownames(which(abs(df4[df4$position==i,c('position',effect_name)][effect_name])== biggest, arr.ind=T)) - df4[df4$position==i,'effect_sign'] = sign(df4[effect_name][ind,]) -} - -df4$effect_type = sub("\\.[0-9]+", "", df4$effect_type) # cull duplicate effect types that happen when there are exact duplicate values -df4U = df4[!duplicated(df4[c('position')]), ] -table(df4U$effect_type) - -#%%============================================================ -# 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 -#=============================================================== diff --git a/scripts/plotting/mcsm_mean_affinity_ensemble.R b/scripts/plotting/mcsm_mean_affinity_ensemble.R deleted file mode 100644 index 068325b..0000000 --- a/scripts/plotting/mcsm_mean_affinity_ensemble.R +++ /dev/null @@ -1,316 +0,0 @@ -#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") -######################################################### -# TASK: Generate averaged affinity values -# across all affinity tools for a given structure -# as applicable... -######################################################### -#============= -# Input -#============= -df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv") -df3 = read.csv(df3_filename) -length(df3$mutationinformation) -all_colnames= colnames(df3) -#%%=============================================================== -# FIXME: ADD distance to NA when SP replies -dist_columns = c("ligand_distance", "interface_dist") -DistCutOff = 10 -common_cols = c("mutationinformation" - , "X5uhc_position" - , "X5uhc_offset" - , "position" - , "dst_mode" - , "mutation_info_labels" - , "sensitivity", dist_columns ) - -all_colnames[grep("scaled" , all_colnames)] -all_colnames[grep("outcome" , all_colnames)] - -#=================== -# stability cols -#=================== -raw_cols_stability = c("duet_stability_change" - , "deepddg" - , "ddg_dynamut2" - , "ddg_foldx") - -scaled_cols_stability = c("duet_scaled" - , "deepddg_scaled" - , "ddg_dynamut2_scaled" - , "foldx_scaled") - -outcome_cols_stability = c("duet_outcome" - , "deepddg_outcome" - , "ddg_dynamut2_outcome" - , "foldx_outcome") - -#=================== -# affinity cols -#=================== -raw_cols_affinity = c("ligand_affinity_change" - , "mmcsm_lig" - , "mcsm_ppi2_affinity" - , "mcsm_na_affinity") - -scaled_cols_affinity = c("affinity_scaled" - , "mmcsm_lig_scaled" - , "mcsm_ppi2_scaled" - , "mcsm_na_scaled" ) - -outcome_cols_affinity = c( "ligand_outcome" - , "mmcsm_lig_outcome" - , "mcsm_ppi2_outcome" - , "mcsm_na_outcome") -#=================== -# conservation cols -#=================== -# raw_cols_conservation = c("consurf_score" -# , "snap2_score" -# , "provean_score") -# -# scaled_cols_conservation = c("consurf_scaled" -# , "snap2_scaled" -# , "provean_scaled") -# -# # CANNOT strictly be used, as categories are not identical with conssurf missing altogether -# outcome_cols_conservation = c("provean_outcome" -# , "snap2_outcome" -# #consurf outcome doesn't exist -# ) -all_cols= c(common_cols - ,raw_cols_stability, scaled_cols_stability, outcome_cols_stability - , raw_cols_affinity, scaled_cols_affinity, outcome_cols_affinity) -#======= -# output -#======= -outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) - -#OutFile1 -outfile_mean_aff = paste0(outdir_images, "/", tolower(gene) - , "_mean_ligand.csv") -print(paste0("Output file:", outfile_mean_aff)) - -#OutFile2 -outfile_ppi2 = paste0(outdir_images, "/", tolower(gene) - , "_mean_ppi2.csv") -print(paste0("Output file:", outfile_ppi2)) - - - -#OutFile4 -#outfile_mean_aff_priorty = paste0(outdir_images, "/", tolower(gene) -# , "_mean_affinity_priority.csv") -#print(paste0("Output file:", outfile_mean_aff_priorty)) - -################################################################# -################################################################# -# mut positions -length(unique(df3$position)) - -# 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)) - -#=============== -# select columns specific to gene -#=============== -gene_aff_cols = colnames(df3)[colnames(df3)%in%c(outcome_cols_affinity - , scaled_cols_affinity)] - -gene_common_cols = colnames(df3)[colnames(df3)%in%common_cols] - -cols_to_extract = c(gene_common_cols - , gene_aff_cols) - -cat("\nExtracting", length(cols_to_extract), "columns") - -df3_plot = df3[, cols_to_extract] -table(df3_plot$mmcsm_lig_outcome) -table(df3_plot$ligand_outcome) - -############################################################## -# mCSM-lig, mCSM-NA, mCSM-ppi2, mmCSM-lig -######################################### -cols_to_numeric = c("ligand_outcome" - , "mcsm_na_outcome" - , "mcsm_ppi2_outcome" - , "mmcsm_lig_outcome") - -#===================================== -# mCSM-lig: Filter ligand distance <10 -#DistCutOff = 10 -#LigDist_colname = "ligand_distance" -# extract outcome cols and map numeric values to the categories -# Destabilising == 0, and stabilising == 1 so rescaling can let -1 be destabilising -#===================================== -df3_lig = df3[, c("mutationinformation" - , "position" - , "ligand_distance" - , "ligand_affinity_change" - , "affinity_scaled" - , "ligand_outcome")] - -df3_lig = df3_lig[df3_lig["ligand_distance"]% - dplyr::group_by(position) %>% - #dplyr::summarize(avg_lig = max(df3_lig_num)) - #dplyr::summarize(avg_lig = mean(ligand_outcome)) - #dplyr::summarize(avg_lig = mean(affinity_scaled, na.rm = T)) - dplyr::summarize(avg_lig = mean(ligand_affinity_change, na.rm = T)) - -class(mean_lig_by_position) - -# convert to a df -mean_lig_by_position = as.data.frame(mean_lig_by_position) -table(mean_lig_by_position$avg_lig) - -# REscale b/w -1 and 1 -lig_min = min(mean_lig_by_position['avg_lig']) -lig_max = max(mean_lig_by_position['avg_lig']) - -mean_lig_by_position['avg_lig_scaled'] = lapply(mean_lig_by_position['avg_lig'] - , function(x) { - scales::rescale_mid(x - , to = c(-1,1) - , from = c(lig_min,lig_max) - , mid = 0) - #, from = c(0,1)) - }) - -cat(paste0('Average (mcsm-lig+mmcsm-lig) scores:\n' - , head(mean_lig_by_position['avg_lig']) - , '\n---------------------------------------------------------------' - , '\nAverage (mcsm-lig+mmcsm-lig) scaled scores:\n' - , head(mean_lig_by_position['avg_lig_scaled']))) - -if ( nrow(mean_lig_by_position) == length(unique(df3_lig$position)) ){ - cat("\nPASS: Generated average values for ligand affinity" ) -}else{ - stop(paste0("\nAbort: length mismatch for ligand affinity data")) -} - -max(mean_lig_by_position$avg_lig); min(mean_lig_by_position$avg_lig) -max(mean_lig_by_position$avg_lig_scaled); min(mean_lig_by_position$avg_lig_scaled) - -################################################################# -# output -write.csv(mean_lig_by_position, outfile_mean_aff - , row.names = F) -cat("Finished writing file:\n" - , outfile_mean_aff - , "\nNo. of rows:", nrow(mean_lig_by_position) - , "\nNo. of cols:", ncol(mean_lig_by_position)) -################################################################## -################################################################## -#===================================== -# mCSM-ppi2: Filter interface_dist <10 -#DistCutOff = 10 - -#===================================== -df3_ppi2 = df3[, c("mutationinformation" - , "position" - , "interface_dist" - , "mcsm_ppi2_affinity" - , "mcsm_ppi2_scaled" - , "mcsm_ppi2_outcome")] - -df3_ppi2 = df3_ppi2[df3_ppi2["interface_dist"]% - dplyr::group_by(position) %>% - #dplyr::summarize(avg_ppi2 = max(df3_ppi2_num)) - #dplyr::summarize(avg_ppi2 = mean(mcsm_ppi2_outcome)) - #dplyr::summarize(avg_ppi2 = mean(mcsm_ppi2_scaled, na.rm = T)) - dplyr::summarize(avg_ppi2 = mean(mcsm_ppi2_affinity, na.rm = T)) - -class(mean_ppi2_by_position) - -# convert to a df -mean_ppi2_by_position = as.data.frame(mean_ppi2_by_position) -table(mean_ppi2_by_position$avg_ppi2) - -# REscale b/w -1 and 1 -lig_min = min(mean_ppi2_by_position['avg_ppi2']) -lig_max = max(mean_ppi2_by_position['avg_ppi2']) - -mean_ppi2_by_position['avg_ppi2_scaled'] = lapply(mean_ppi2_by_position['avg_ppi2'] - , function(x) { - scales::rescale_mid(x - , to = c(-1,1) - , from = c(lig_min,lig_max) - , mid = 0) - #, from = c(0,1)) - }) - -cat(paste0('Average ppi2 scores:\n' - , head(mean_ppi2_by_position['avg_ppi2']) - , '\n---------------------------------------------------------------' - , '\nAverage ppi2 scaled scores:\n' - , head(mean_ppi2_by_position['avg_ppi2_scaled']))) - -if ( nrow(mean_ppi2_by_position) == length(unique(df3_ppi2$position)) ){ - cat("\nPASS: Generated average values for ppi2" ) -}else{ - stop(paste0("\nAbort: length mismatch for ppi2 data")) -} - -max(mean_ppi2_by_position$avg_ppi2); min(mean_ppi2_by_position$avg_ppi2) -max(mean_ppi2_by_position$avg_ppi2_scaled); min(mean_ppi2_by_position$avg_ppi2_scaled) - - -write.csv(mean_ppi2_by_position, outfile_ppi2 - , row.names = F) -cat("Finished writing file:\n" - , outfile_ppi2 - , "\nNo. of rows:", nrow(mean_ppi2_by_position) - , "\nNo. of cols:", ncol(mean_ppi2_by_position)) - - -# end of script -#=============================================================== diff --git a/scripts/plotting/mcsm_mean_stability.R b/scripts/plotting/mcsm_mean_stability.R deleted file mode 100755 index 7e35e79..0000000 --- a/scripts/plotting/mcsm_mean_stability.R +++ /dev/null @@ -1,163 +0,0 @@ -getwd() -setwd("~/git/LSHTM_analysis/scripts/plotting") -getwd() - -######################################################### -# TASK: - -######################################################### -#source("~/git/LSHTM_analysis/scripts/Header_TT.R") -#require(data.table) -#require(dplyr) - -source("plotting_data.R") -# should return -#my_df -#my_df_u -#dup_muts - -# cmd parse arguments -#require('getopt', quietly = TRUE) -#======================================================== - - -#======================================================== -# Read file: call script for combining df for PS - -#source("../combining_two_df.R") - -#======================================================== - -# plotting_data.R imports all the dir names, etc - -#======= -# output -#======= -out_filename_mean_stability = paste0(tolower(gene), "_mean_stability.csv") -outfile_mean_stability = paste0(outdir, "/", out_filename_mean_stability) -print(paste0("Output file:", outfile_mean_stability)) - -#%%=============================================================== - -#================ -# Data for plots -#================ -# REASSIGNMENT as necessary -df = my_df_u -rm(my_df) - -########################### -# Data for bfactor figure -# PS (duet) average -# Ligand affinity average -########################### -head(df$position); head(df$mutationinformation) -head(df$duet_stability_change) - -# order data frame -#df = df[order(df$position),] #already done -#head(df$position); head(df$mutationinformation) -#head(df$duet_stability_change) - -#*********** -# PS(duet): average by position and then scale b/w -1 and 1 -# column to average: duet_stability_change (NOT scaled!) -#*********** -mean_duet_by_position <- df %>% - group_by(position) %>% - summarize(averaged_duet = mean(duet_stability_change)) - -# scale b/w -1 and 1 -duet_min = min(mean_duet_by_position['averaged_duet']) -duet_max = max(mean_duet_by_position['averaged_duet']) - -# scale the averaged_duet values -mean_duet_by_position['averaged_duet_scaled'] = lapply(mean_duet_by_position['averaged_duet'] - , function(x) ifelse(x < 0, x/abs(duet_min), x/duet_max)) - -cat(paste0('Average duet scores:\n', head(mean_duet_by_position['averaged_duet']) - , '\n---------------------------------------------------------------' - , '\nScaled duet scores:\n', head(mean_duet_by_position['averaged_duet_scaled']))) - -# sanity checks -l_bound_duet = min(mean_duet_by_position['averaged_duet_scaled']) -u_bound_duet = max(mean_duet_by_position['averaged_duet_scaled']) - -if ( (l_bound_duet == -1) && (u_bound_duet == 1) ){ - cat(paste0("PASS: duet scores averaged by position and then scaled" - , "\nmin averaged duet: ", l_bound_duet - , "\nmax averaged duet: ", u_bound_duet)) -}else{ - cat(paste0("FAIL: avergaed duet scores could not be scaled b/w -1 and 1" - , "\nmin averaged duet: ", l_bound_duet - , "\nmax averaged duet: ", u_bound_duet)) - quit() -} - -#*********** -# Lig: average by position and then scale b/w -1 and 1 -# column: ligand_affinity_change (NOT scaled!) -#*********** -mean_affinity_by_position <- df %>% - group_by(position) %>% - summarize(averaged_affinity = mean(ligand_affinity_change)) - -# scale b/w -1 and 1 -affinity_min = min(mean_affinity_by_position['averaged_affinity']) -affinity_max = max(mean_affinity_by_position['averaged_affinity']) - -# scale the averaged_affinity values -mean_affinity_by_position['averaged_affinity_scaled'] = lapply(mean_affinity_by_position['averaged_affinity'] - , function(x) ifelse(x < 0, x/abs(affinity_min), x/affinity_max)) - -cat(paste0('Average affinity scores:\n', head(mean_affinity_by_position['averaged_affinity']) - , '\n---------------------------------------------------------------' - , '\nScaled affinity scores:\n', head(mean_affinity_by_position['averaged_affinity_scaled']))) - -# sanity checks -l_bound_affinity = min(mean_affinity_by_position['averaged_affinity_scaled']) -u_bound_affinity = max(mean_affinity_by_position['averaged_affinity_scaled']) - -if ( (l_bound_affinity == -1) && (u_bound_affinity == 1) ){ - cat(paste0("PASS: affinity scores averaged by position and then scaled" - , "\nmin averaged affintiy: ", l_bound_affinity - , "\nmax averaged affintiy: ", u_bound_affinity)) -}else{ - cat(paste0("FAIL: avergaed affinity scores could not be scaled b/w -1 and 1" - , "\nmin averaged affintiy: ", l_bound_affinity - , "\nmax averaged affintiy: ", u_bound_affinity)) - quit() -} - -#*********** -# merge: mean_duet_by_position and mean_affinity_by_position -#*********** -common_cols = intersect(colnames(mean_duet_by_position), colnames(mean_affinity_by_position)) - -if (dim(mean_duet_by_position) && dim(mean_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_duet_by_position - , mean_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_stability - , row.names = F) -cat("Finished writing file:\n" - , outfile_mean_stability - , "\nNo. of rows:", nrow(combined_df) - , "\nNo. of cols:", ncol(combined_df)) - -# end of script -#=============================================================== diff --git a/scripts/plotting/mcsm_mean_stability_ensemble.R b/scripts/plotting/mcsm_mean_stability_ensemble.R deleted file mode 100644 index aeb5d0a..0000000 --- a/scripts/plotting/mcsm_mean_stability_ensemble.R +++ /dev/null @@ -1,212 +0,0 @@ -#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") -######################################################### -# 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_stability.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) -length(df3$mutationinformation) - -# 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" - , "interface_dist") - -all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)] -all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)] - -#=================== -# stability cols -#=================== -raw_cols_stability = c("duet_stability_change" - , "deepddg" - , "ddg_dynamut2" - , "ddg_foldx") - -scaled_cols_stability = c("duet_scaled" - , "deepddg_scaled" - , "ddg_dynamut2_scaled" - , "foldx_scaled") - -outcome_cols_stability = c("duet_outcome" - , "deepddg_outcome" - , "ddg_dynamut2_outcome" - , "foldx_outcome") - -#=================== -# affinity cols -#=================== -raw_cols_affinity = c("ligand_affinity_change" - , "mmcsm_lig" - , "mcsm_ppi2_affinity" - , "mcsm_na_affinity") - -scaled_cols_affinity = c("affinity_scaled" - , "mmcsm_lig_scaled" - , "mcsm_ppi2_scaled" - , "mcsm_na_scaled" ) - -outcome_cols_affinity = c( "ligand_outcome" - , "mmcsm_lig_outcome" - , "mcsm_ppi2_outcome" - , "mcsm_na_outcome") - -#=================== -# conservation cols -#=================== -# raw_cols_conservation = c("consurf_score" -# , "snap2_score" -# , "provean_score") -# -# scaled_cols_conservation = c("consurf_scaled" -# , "snap2_scaled" -# , "provean_scaled") -# -# # CANNOT strictly be used, as categories are not identical with conssurf missing altogether -# outcome_cols_conservation = c("provean_outcome" -# , "snap2_outcome" -# #consurf outcome doesn't exist -# ) - -########################################################### -cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols - , raw_cols_stability - , scaled_cols_stability - , outcome_cols_stability - , raw_cols_affinity - , scaled_cols_affinity - , outcome_cols_affinity)] - -cols_to_extract = cols_to_consider[cols_to_consider%in%c(common_cols - , outcome_cols_stability)] -############################################################## -##################### -# Ensemble stability: outcome_cols_stability -##################### -# extract outcome cols and map numeric values to the categories -# 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_stability] <- sapply(df3_plot[, outcome_cols_stability] - , 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 -# position ==> scale b/w -1 and 1 - -# column to average: ens_stability -#===================================== -cols_to_average = which(colnames(df3_plot)%in%outcome_cols_stability) - -# 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() -# } -################################################################## -# output -#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(mean_ens_stability_by_position) - , "\nNo. of cols:", ncol(mean_ens_stability_by_position)) - -# end of script -#=============================================================== diff --git a/scripts/plotting/mcsm_mean_stability_ensemble_5uhc_rpob.R b/scripts/plotting/mcsm_mean_stability_ensemble_5uhc_rpob.R deleted file mode 100644 index c75a6ba..0000000 --- a/scripts/plotting/mcsm_mean_stability_ensemble_5uhc_rpob.R +++ /dev/null @@ -1,176 +0,0 @@ -#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") -######################################################### -# 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, "/5uhc_", tolower(gene) - , "_mean_ens_stability.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) -length(df3$mutationinformation) - -# 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" - , "X5uhc_position" - , "dst_mode" - , "mutation_info_labels" - , "sensitivity" - , "X5uhc_position" - , "X5uhc_offset" - , "ligand_distance" - , "interface_dist") - -all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)] -all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)] - -#=================== -# stability cols -#=================== -raw_cols_stability = c("duet_stability_change" - , "deepddg" - , "ddg_dynamut2" - , "ddg_foldx") - -scaled_cols_stability = c("duet_scaled" - , "deepddg_scaled" - , "ddg_dynamut2_scaled" - , "foldx_scaled") - -outcome_cols_stability = c("duet_outcome" - , "deepddg_outcome" - , "ddg_dynamut2_outcome" - , "foldx_outcome") - -########################################################### -cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols - , raw_cols_stability - , scaled_cols_stability - , outcome_cols_stability)] - -cols_to_extract = cols_to_consider[cols_to_consider%in%c(common_cols - , outcome_cols_stability)] -############################################################## -##################### -# Ensemble stability: outcome_cols_stability -##################### -# extract outcome cols and map numeric values to the categories -# 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_stability] <- sapply(df3_plot[, outcome_cols_stability] - , 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 -# X5uhc_position ==> scale b/w -1 and 1 - -# column to average: ens_stability -#===================================== -cols_to_average = which(colnames(df3_plot)%in%outcome_cols_stability) - -# ensemble average across predictors -df3_plot$ens_stability = rowMeans(df3_plot[,cols_to_average]) - -head(df3_plot$X5uhc_position); head(df3_plot$mutationinformation) -head(df3_plot$ens_stability) -table(df3_plot$ens_stability) - -# ensemble average of predictors by X5uhc_position -mean_ens_stability_by_position <- df3_plot %>% - dplyr::group_by(X5uhc_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 X5uhc_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() -# } -################################################################## -# output -#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(mean_ens_stability_by_position) - , "\nNo. of cols:", ncol(mean_ens_stability_by_position)) - -# end of script -#=============================================================== diff --git a/scripts/plotting/mut_landscape.R b/scripts/plotting/mut_landscape.R deleted file mode 100644 index aeeb5c2..0000000 --- a/scripts/plotting/mut_landscape.R +++ /dev/null @@ -1,155 +0,0 @@ -source("~/git/LSHTM_analysis/config/alr.R") -source("~/git/LSHTM_analysis/config/embb.R") -source("~/git/LSHTM_analysis/config/gid.R") -source("~/git/LSHTM_analysis/config/katg.R") -source("~/git/LSHTM_analysis/config/pnca.R") -source("~/git/LSHTM_analysis/config/rpob.R") - -#================================ -# output files -# In total: 6 files are written -#================================ -outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) - -# mutational positions: all -outfile_mutpos = paste0(outdir_images, "/", tolower(gene), "_mutpos_all.txt") -outfile_meta1 = paste0(outdir_images, "/", tolower(gene), "_mutpos_cu.txt") - -# mutational positions with sensitivity: S, R and common -outfile_mutpos_S = paste0(outdir_images, "/", tolower(gene), "_mutpos_S.txt") -outfile_mutpos_R = paste0(outdir_images, "/", tolower(gene), "_mutpos_R.txt") -outfile_mutpos_common = paste0(outdir_images, "/", tolower(gene), "_mutpos_common.txt") -outfile_meta2 = paste0(outdir_images, "/", tolower(gene), "_mutpos_annot_cu.txt") - -#============= -# Input -#============= -df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv") -df3 = read.csv(df3_filename) - -# Determine for each gene -if (tolower(gene) == "embb"){ - chain_suffix = ".B" -} else{ - chain_suffix = ".A" -} - -# 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)) - -############################################################ -cols_to_extract = c("mutationinformation" - , "wild_type" - , "chain" - , "mutant_type" - , "position" - , "dst_mode" - , "mutation_info_labels_orig" - , "mutation_info_labels" - , "sensitivity") - -df3_plot = df3[, cols_to_extract] -# create pos_chain column: allows easier colouring in chimera -df3_plot$pos_chain = paste(df3_plot$position, df3_plot$chain, sep = ".") -pos_cu = length(unique(df3_plot$position)) - -#=========================== -# positions with mutations -#=========================== -#df3_all_mut_pos = df3_plot[, c("mutationinformation", "position", "chain")] -#df3_all_mut_pos$pos_chain = paste(df3_all_mut_pos$position, df3_all_mut_pos$chain, sep = ".") - -df3_all_mut_pos = df3_plot[, c("position", "pos_chain")] -gene_mut_pos_u = unique(df3_all_mut_pos$pos_chain) -class(gene_mut_pos_u) -paste(gene_mut_pos_u, collapse=',') - -if (length(gene_mut_pos_u) == pos_cu){ - cat("\nPASS: all mutation positions extracted" - , "\nWriting file:", outfile_mutpos) -} else{ - stop("\nAbort: mutation position count mismatch") -} - -write.table(paste(gene_mut_pos_u, collapse=',') - , outfile_mutpos - , row.names = F - , col.names = F) - -write.table(paste("Count of positions with mutations in gene" - , tolower(gene), ":", pos_cu) - , outfile_meta1 - , row.names = F - , col.names = F) -#======================================== -# positions with sensitivity annotations -#======================================== -df3_muts_annot = df3_plot[, c("mutationinformation", "position", "pos_chain", "sensitivity")] - -# aggregrate position count by sensitivity -result <- aggregate(sensitivity ~ position, data = df3_muts_annot, paste, collapse = "") - -sensitive_pos = result$position[grep("(^S+$)", result$sensitivity)] -sensitive_pos = paste0(sensitive_pos, chain_suffix) - -resistant_pos = result$position[grep("(^R+$)", result$sensitivity)] -resistant_pos = paste0(resistant_pos, chain_suffix) - -common_pos = result$position[grep("SR|RS" , result$sensitivity)] -common_pos = paste0(common_pos, chain_suffix) - -if (tolower(gene)!= "alr"){ - length_check = length(sensitive_pos) + length(resistant_pos) + length(common_pos) - cpl = length(common_pos) -}else{ - length_check = length(sensitive_pos) + length(resistant_pos) - cpl = 0 -} - -if (length_check == pos_cu){ - cat("\nPASS: position with mutational sensitivity extracted" - , "\nWriting files: sensitive, resistant and common position counts" ) -} else{ - stop("\nAbort: position with mutational sensitivity count mismatch") -} - -write.table(paste(sensitive_pos, collapse = ',') - , outfile_mutpos_S - , row.names = F, col.names = F) - -write.table(paste(resistant_pos, collapse = ',') - , outfile_mutpos_R - , row.names = F, col.names = F) - -write.table(paste(common_pos, collapse = ',') - , outfile_mutpos_common - , row.names = F, col.names = F) - -write.table(paste("Count of positions with mutations in gene:" - , tolower(gene) - , "\nTotal mutational positions:", pos_cu - , "\nsensitive:", length(sensitive_pos) - , "\nresistant:", length(resistant_pos) - , "\ncommon:" , cpl) - , outfile_meta2 - , row.names = F - , col.names = F) diff --git a/scripts/plotting/mut_landscape_5uhc_rpob.R b/scripts/plotting/mut_landscape_5uhc_rpob.R deleted file mode 100644 index 53f95d0..0000000 --- a/scripts/plotting/mut_landscape_5uhc_rpob.R +++ /dev/null @@ -1,191 +0,0 @@ -source("~/git/LSHTM_analysis/config/rpob.R") - -#================================ -# output files -# In total: 6 files are written -#================================ -outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) - -# mutational positions: all -outfile_mutpos = paste0(outdir_images, "/5uhc_", tolower(gene), "_mutpos_all.txt") -outfile_meta1 = paste0(outdir_images, "/5uhc_", tolower(gene), "_mutpos_cu.txt") - -# mutational positions with sensitivity: S, R and common -outfile_mutpos_S = paste0(outdir_images, "/5uhc_", tolower(gene), "_mutpos_S.txt") -outfile_mutpos_R = paste0(outdir_images, "/5uhc_", tolower(gene), "_mutpos_R.txt") -outfile_mutpos_common = paste0(outdir_images, "/5uhc_", tolower(gene), "_mutpos_common.txt") -outfile_meta2 = paste0(outdir_images, "/5uhc_", tolower(gene), "_mutpos_annot_cu.txt") - -#============= -# Input -#============= -df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv") -df3 = read.csv(df3_filename) - -chain_suffix = ".C" - -# 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)) - -############################################################ -cols_to_extract = c("mutationinformation" - , "wild_type" - , "chain" - , "mutant_type" - , "position" - , "X5uhc_position" - , "X5uhc_offset" - , "dst_mode" - , "mutation_info_labels_orig" - , "mutation_info_labels" - , "sensitivity") - -df3_plot = df3[, cols_to_extract] - -# use the x5uhc_position column - -# create pos_chain column: allows easier colouring in chimera -df3_plot$pos_chain = paste(df3_plot$X5uhc_position, chain_suffix, sep = ".") -pos_cu = length(unique(df3_plot$X5uhc_position)) -X5uhc_pos = unique(df3_plot$X5uhc_position) -X5uhc_pos = paste0(X5uhc_pos, chain_suffix) - -#=========================== -# positions with mutations -#=========================== -df3_all_mut_pos = df3_plot[, c("X5uhc_position", "pos_chain")] -gene_mut_pos_u = unique(df3_all_mut_pos$pos_chain) -class(gene_mut_pos_u) -paste(gene_mut_pos_u, collapse=',') - -if (length(gene_mut_pos_u) == pos_cu){ - cat("\nPASS: all mutation positions extracted" - , "\nWriting file:", outfile_mutpos) -} else{ - stop("\nAbort: mutation position count mismatch") -} - -write.table(paste(gene_mut_pos_u, collapse=',') - , outfile_mutpos - , row.names = F - , col.names = F) - -write.table(paste("Count of positions with mutations in gene" - , tolower(gene), ":", pos_cu) - , outfile_meta1 - , row.names = F - , col.names = F) -#======================================== -# positions with sensitivity annotations -#======================================== -df3_muts_annot = df3_plot[, c("mutationinformation", "X5uhc_position", "pos_chain", "sensitivity")] - -# aggregrate position count by sensitivity -result <- aggregate(sensitivity ~ X5uhc_position, data = df3_muts_annot, paste, collapse = "") - -sensitive_pos = result$X5uhc_position[grep("(^S+$)", result$sensitivity)] -sensitive_pos = paste0(sensitive_pos, chain_suffix) - -resistant_pos = result$X5uhc_position[grep("(^R+$)", result$sensitivity)] -resistant_pos = paste0(resistant_pos, chain_suffix) - -common_pos = result$X5uhc_position[grep("SR|RS" , result$sensitivity)] -common_pos = paste0(common_pos, chain_suffix) - -if (tolower(gene)!= "alr"){ - length_check = length(sensitive_pos) + length(resistant_pos) + length(common_pos) - cpl = length(common_pos) -}else{ - length_check = length(sensitive_pos) + length(resistant_pos) - cpl = 0 -} - -if (length_check == pos_cu){ - cat("\nPASS: position with mutational sensitivity extracted" - , "\nWriting files: sensitive, resistant and common position counts" ) -} else{ - stop("\nAbort: position with mutational sensitivity count mismatch") -} -# spl handling for rpob 5uhc -revised_gene_mut_pos_u = c(sensitive_pos, resistant_pos, common_pos) -revised_pos_cu = length(unique(revised_gene_mut_pos_u)) -if (length(revised_gene_mut_pos_u) == revised_pos_cu){ - cat("\nPASS: all mutation positions extracted" - , "\nWriting file:", outfile_mutpos) -} else{ - stop("\nAbort: mutation position count mismatch") -} - -write.table(paste(revised_gene_mut_pos_u, collapse=',') - , outfile_mutpos - , row.names = F - , col.names = F) - -write.table(paste("Count of positions with mutations in gene" - , tolower(gene), ":", revised_pos_cu) - , outfile_meta1 - , row.names = F - , col.names = F) - -# mut_annot -write.table(paste(sensitive_pos, collapse = ',') - , outfile_mutpos_S - , row.names = F, col.names = F) - -write.table(paste(resistant_pos, collapse = ',') - , outfile_mutpos_R - , row.names = F, col.names = F) - -write.table(paste(common_pos, collapse = ',') - , outfile_mutpos_common - , row.names = F, col.names = F) - -write.table(paste("Count of positions with mutations in gene:" - , tolower(gene) - , "\nTotal mutational positions:", revised_pos_cu - , "\nsensitive:", length(sensitive_pos) - , "\nresistant:", length(resistant_pos) - , "\ncommon:" , cpl) - , outfile_meta2 - , row.names = F - , col.names = F) - -#Quick check to find out the discrepancy -revised_gene_mut_pos_u -gene_mut_pos_u - -library("qpcR") -foo <- data.frame(qpcR:::cbind.na(gene_mut_pos_u, revised_gene_mut_pos_u)) - - -table(!gene_mut_pos_u%in%revised_gene_mut_pos_u) -table(!revised_gene_mut_pos_u%in%gene_mut_pos_u) - -X5uhc_pos -#table(!gene_mut_pos_u%in%X5uhc_pos) -table(X5uhc_pos%in%gene_mut_pos_u) - -X5uhc_pos[!X5uhc_pos%in%gene_mut_pos_u] -X5uhc_pos[!gene_mut_pos_u%in%X5uhc_pos] - -#TODO: NOTE -#D1148G (i.e D1154) is NOT Present in 5UHC diff --git a/scripts/plotting/replaceBfactor_pdb.R b/scripts/plotting/replaceBfactor_pdb.R deleted file mode 100755 index a840c19..0000000 --- a/scripts/plotting/replaceBfactor_pdb.R +++ /dev/null @@ -1,332 +0,0 @@ -#!/usr/bin/env Rscript - -######################################################### -# TASK: Replace B-factors in the pdb file with the mean -# normalised stability values. - -# read pdb file -# make two copies so you can replace B factors for 1)duet -# 2)affinity values and output 2 separate pdbs for -# rendering on chimera - -# read mcsm mean stability value files -# extract the respective mean values and assign to the -# b-factor column within their respective pdbs - -# generate some distribution plots for inspection - -######################################################### -# working dir and loading libraries -getwd() -setwd("~/git/LSHTM_analysis/scripts/plotting") -cat(c(getwd(),"\n")) - -#source("~/git/LSHTM_analysis/scripts/Header_TT.R") -library(bio3d) -require("getopt", quietly = TRUE) # cmd parse arguments -#======================================================== -#drug = "pyrazinamide" -#gene = "pncA" - -# command line args -spec = matrix(c( - "drug" , "d", 1, "character", - "gene" , "g", 1, "character" -), byrow = TRUE, ncol = 4) - -opt = getopt(spec) - -drug = opt$drug -gene = opt$gene - -if(is.null(drug)|is.null(gene)) { - stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)") -} -#======================================================== -gene_match = paste0(gene,"_p.") -cat(gene_match) - -#============= -# directories -#============= -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/Writing/thesis/images/results/", tolower(gene)) - -#====== -# input -#====== -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_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_b_stab_ms.pdb") -outfile_duet_mspdb = paste0(outdir_plots, "/", out_filename_duet_mspdb) -print(paste0("Output file:", outfile_duet_mspdb)) - -out_filename_lig_mspdb = paste0(tolower(gene), "_complex_blig_ms.pdb") -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 -########################### -my_df <- read.csv(infile_mean_stability, header = T) -str(my_df) - -############# -# Read pdb -############# -# list of 8 -my_pdb = read.pdb(infile_pdb - , maxlines = -1 - , multi = FALSE - , rm.insert = FALSE - , rm.alt = TRUE - , ATOM.only = FALSE - , hex = FALSE - , verbose = TRUE) - -rm(in_filename_mean_stability, in_filename_pdb) - -# assign separately for duet and ligand -my_pdb_duet = my_pdb -my_pdb_lig = my_pdb - -#========================================================= -# Replacing B factor with mean stability scores -# within the respective dfs -#========================================================== -# extract atom list into a variable -# since in the list this corresponds to data frame, variable will be a df -#df_duet = my_pdb_duet[[1]] -df_duet= my_pdb_duet[['atom']] -df_lig = my_pdb_lig[['atom']] - -# make a copy: required for downstream sanity checks -d2_duet = df_duet -d2_lig = df_lig - -# sanity checks: B factor -max(df_duet$b); min(df_duet$b) -max(df_lig$b); min(df_lig$b) - -#******************************************* -# histograms and density plots for inspection -# 1: original B-factors -# 2: original mean stability values -# 3: replaced B-factors with mean stability values -#********************************************* -# Set the margin on all sides -par(oma = c(3,2,3,0) - , mar = c(1,3,5,2) - #, mfrow = c(3,2) - , mfrow = c(3,4)) - -#============= -# Row 1 plots: original B-factors -# duet and affinity -#============= -hist(df_duet$b - , xlab = "" - , main = "Bfactor stability") - -plot(density(df_duet$b) - , xlab = "" - , main = "Bfactor stability") - - -hist(df_lig$b - , xlab = "" - , main = "Bfactor affinity") - -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$avg_ens_stability_scaled - , xlab = "" - , main = "mean stability values") - -#plot(density(my_df$averaged_duet) -plot(density(my_df$avg_ens_stability_scaled) - , xlab = "" - , main = "mean stability values") - -#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$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 -#========= -# sanity check: match and assign to a separate column to double check -# colnames(my_df) -# df_duet$duet_scaled = my_df$averge_duet_scaled[match(df_duet$resno, my_df$position)] - -#========= -# step 1_P1 -#========= -# 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$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 - -# count number of 0"s in Bactor -sum(df_duet$b == 0) -sum(df_lig$b == 0) - -# replace all NA in b factor with 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 == 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() -# } - -# 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 -#========= -# sanity check: dim should be same before reassignment -if ( (dim(df_duet)[1] == dim(d2_duet)[1]) & (dim(df_lig)[1] == dim(d2_lig)[1]) & - (dim(df_duet)[2] == dim(d2_duet)[2]) & (dim(df_lig)[2] == dim(d2_lig)[2]) - ){ - print("PASS: Dims of both dfs as expected") -} else { - print ("FAIL: Dims mismatch") - quit()} - -#========= -# step 4_P1: -# VERY important -#========= -# 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) - -#========= -# step 5_P1 -#========= -cat(paste0("output file duet mean stability pdb:", outfile_duet_mspdb)) -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 - , xlab = "" - , main = "repalcedB duet") - -plot(density(df_duet$b) - , xlab = "" - , main = "replacedB duet") - - -hist(df_lig$b - , xlab = "" - , main = "repalcedB affinity") - -plot(density(df_lig$b) - , xlab = "" - , main = "replacedB affinity") - -# graph titles -mtext(text = "Frequency" - , side = 2 - , line = 0 - , outer = TRUE) - -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/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??? -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - diff --git a/scripts/plotting/replaceBfactor_pdb_affinity.R b/scripts/plotting/replaceBfactor_pdb_affinity.R deleted file mode 100644 index 1ba5921..0000000 --- a/scripts/plotting/replaceBfactor_pdb_affinity.R +++ /dev/null @@ -1,281 +0,0 @@ -#!/usr/bin/env Rscript - -######################################################### -# TASK: Replace B-factors in the pdb file with the mean -# normalised stability values. - -# read pdb file - -# read mcsm mean stability value files -# extract the respective mean values and assign to the -# b-factor column within their respective pdbs - -# generate some distribution plots for inspection - -######################################################### -# working dir and loading libraries -getwd() -setwd("~/git/LSHTM_analysis/scripts/plotting") -cat(c(getwd(),"\n")) - -#source("~/git/LSHTM_analysis/scripts/Header_TT.R") -library(bio3d) -require("getopt", quietly = TRUE) # cmd parse arguments -#======================================================== -#drug = "pyrazinamide" -#gene = "pncA" - -# command line args -spec = matrix(c( - "drug" , "d", 1, "character", - "gene" , "g", 1, "character" -), byrow = TRUE, ncol = 4) - -opt = getopt(spec) - -drug = opt$drug -gene = opt$gene - -if(is.null(drug)|is.null(gene)) { - stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)") -} -#======================================================== -gene_match = paste0(gene,"_p.") -cat(gene_match) - -#============= -# directories -#============= -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/Writing/thesis/images/results/", tolower(gene)) - -#====== -# input -#====== -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_affinity = paste0(tolower(gene), "_mean_ligand.csv") -infile_mean_affinity = paste0(outdir_plots, "/", in_filename_mean_affinity) - -cat(paste0("Input file:", infile_mean_affinity) ) - -#======= -# output -#======= -#out_filename_duet_mspdb = paste0(tolower(gene), "_complex_bduet_ms.pdb") -out_filename_lig_mspdb = paste0(tolower(gene), "_complex_b_lig_ms.pdb") -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 -########################### -my_df <- read.csv(infile_mean_stability, header = T) -str(my_df) - -############# -# Read pdb -############# -# list of 8 -my_pdb = read.pdb(infile_pdb - , maxlines = -1 - , multi = FALSE - , rm.insert = FALSE - , rm.alt = TRUE - , ATOM.only = FALSE - , hex = FALSE - , verbose = TRUE) - -rm(in_filename_mean_affinity, in_filename_pdb) - -# assign separately for duet and ligand -my_pdb_duet = my_pdb - -#========================================================= -# Replacing B factor with mean stability scores -# within the respective dfs -#========================================================== -# extract atom list into a variable -# since in the list this corresponds to data frame, variable will be a df -#df_duet = my_pdb_duet[[1]] -df_duet= my_pdb_duet[['atom']] - -# make a copy: required for downstream sanity checks -d2_duet = df_duet - -# sanity checks: B factor -max(df_duet$b); min(df_duet$b) - -#================================================== -# histograms and density plots for inspection -# 1: original B-factors -# 2: original mean stability values -# 3: replaced B-factors with mean stability values -#================================================== -# Set the margin on all sides -par(oma = c(3,2,3,0) - , mar = c(1,3,5,2) - #, mfrow = c(3,2) - #, mfrow = c(3,4)) - , mfrow = c(3,2)) - -#============= -# Row 1 plots: original B-factors -# duet and affinity -#============= -hist(df_duet$b - , xlab = "" - , main = "Bfactor affinity") - -plot(density(df_duet$b) - , xlab = "" - , main = "Bfactor affinity") - -#============= -# Row 2 plots: original mean stability values -# duet and affinity -#============= - -#hist(my_df$averaged_duet -hist(my_df$avg_lig_scaled - , xlab = "" - , main = "mean affinity values") - -#plot(density(my_df$averaged_duet) -plot(density(my_df$avg_lig_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 -#========= -# sanity check: match and assign to a separate column to double check -# colnames(my_df) -# df_duet$duet_scaled = my_df$averge_duet_scaled[match(df_duet$resno, my_df$position)] - -#========= -# step 1_P1 -#========= -# 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_duet$b = my_df$avg_lig_scaled[match(df_duet$resno, my_df$position)] - -#========= -# step 2_P1 -#========= -# count NA in Bfactor -b_na_duet = sum(is.na(df_duet$b)) ; b_na_duet - -# count number of 0"s in Bactor -sum(df_duet$b == 0) - -# replace all NA in b factor with 0 -na_rep = 2 -df_duet$b[is.na(df_duet$b)] = na_rep - -# # 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() -# } - -# 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 -#========= -# sanity check: dim should be same before reassignment -if ( (dim(df_duet)[1] == dim(d2_duet)[1]) & - (dim(df_duet)[2] == dim(d2_duet)[2]) - ){ - print("PASS: Dims of both dfs as expected") -} else { - print ("FAIL: Dims mismatch") - quit()} - -#========= -# step 4_P1: -# VERY important -#========= -# 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)) - -#========= -# step 5_P1 -#========= -cat(paste0("output file duet mean stability pdb:" - , outfile_lig_mspdb)) -write.pdb(my_pdb_duet, 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 - , xlab = "" - , main = "repalcedB duet") - -plot(density(df_duet$b) - , xlab = "" - , main = "replacedB duet") - -# graph titles -mtext(text = "Frequency" - , side = 2 - , line = 0 - , outer = TRUE) - -mtext(text = paste0(tolower(gene), ": afinity 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/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??? -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! \ No newline at end of file diff --git a/scripts/plotting/replaceBfactor_pdb_ppi2.R b/scripts/plotting/replaceBfactor_pdb_ppi2.R deleted file mode 100644 index 03efa38..0000000 --- a/scripts/plotting/replaceBfactor_pdb_ppi2.R +++ /dev/null @@ -1,277 +0,0 @@ -#!/usr/bin/env Rscript - -######################################################### -# TASK: Replace B-factors in the pdb file with the mean -# normalised stability values. - -# read pdb file - -# read mcsm mean stability value files -# extract the respective mean values and assign to the -# b-factor column within their respective pdbs - -# generate some distribution plots for inspection - -######################################################### -# working dir and loading libraries -getwd() -setwd("~/git/LSHTM_analysis/scripts/plotting") -cat(c(getwd(),"\n")) - -#source("~/git/LSHTM_analysis/scripts/Header_TT.R") -library(bio3d) -require("getopt", quietly = TRUE) # cmd parse arguments -#======================================================== -#drug = "pyrazinamide" -#gene = "pncA" - -# command line args -spec = matrix(c( - "drug" , "d", 1, "character", - "gene" , "g", 1, "character" -), byrow = TRUE, ncol = 4) - -opt = getopt(spec) - -drug = opt$drug -gene = opt$gene - -if(is.null(drug)|is.null(gene)) { - stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)") -} -#======================================================== -gene_match = paste0(gene,"_p.") -cat(gene_match) - -#============= -# directories -#============= -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/Writing/thesis/images/results/", tolower(gene)) - -#====== -# input -#====== -in_filename_pdb = paste0(tolower(gene), "_complex.pdb") -infile_pdb = paste0(indir, "/", in_filename_pdb) -cat(paste0("Input file:", infile_pdb) ) - -# mean ppi2 -in_filename_mean_ppi2 = paste0(tolower(gene), "_mean_ppi2.csv") -infile_mean_ppi2 = paste0(outdir_plots, "/", in_filename_mean_ppi2) - -cat(paste0("Input file:", infile_mean_ppi2) ) - -#======= -# output -#======= -#out_filename_duet_mspdb = paste0(tolower(gene), "_complex_bduet_ms.pdb") -out_filename_ppi2_mspdb = paste0(tolower(gene), "_complex_b_ppi2_ms.pdb") -outfile_ppi2_mspdb = paste0(outdir_plots, "/", out_filename_ppi2_mspdb) -print(paste0("Output file:", outfile_ppi2_mspdb)) - -#%%=============================================================== -#NOTE: duet here refers to the ensemble stability values - -########################### -# Read file: average stability values -# or mcsm_normalised file -########################### -my_df <- read.csv(infile_mean_ppi2, header = T) -str(my_df) - -############# -# Read pdb -############# -# list of 8 -my_pdb = read.pdb(infile_pdb - , maxlines = -1 - , multi = FALSE - , rm.insert = FALSE - , rm.alt = TRUE - , ATOM.only = FALSE - , hex = FALSE - , verbose = TRUE) - -# assign separately for duet and ligand -my_pdb_duet = my_pdb - -#========================================================= -# Replacing B factor with mean stability scores -# within the respective dfs -#========================================================== -# extract atom list into a variable -# since in the list this corresponds to data frame, variable will be a df -#df_duet = my_pdb_duet[[1]] -df_duet= my_pdb_duet[['atom']] - -# make a copy: required for downstream sanity checks -d2_duet = df_duet - -# sanity checks: B factor -max(df_duet$b); min(df_duet$b) - -#================================================== -# histograms and density plots for inspection -# 1: original B-factors -# 2: original mean stability values -# 3: replaced B-factors with mean stability values -#================================================== -# Set the margin on all sides -par(oma = c(3,2,3,0) - , mar = c(1,3,5,2) - #, mfrow = c(3,2) - #, mfrow = c(3,4)) - , mfrow = c(3,2)) - -#============= -# Row 1 plots: original B-factors -# duet and affinity -#============= -hist(df_duet$b - , xlab = "" - , main = "Bfactor ppi2") - -plot(density(df_duet$b) - , xlab = "" - , main = "Bfactor ppi2") - -#============= -# Row 2 plots: original mean stability values -# duet and affinity -#============= - -#hist(my_df$averaged_duet -hist(my_df$avg_ppi2_scaled - , xlab = "" - , main = "mean ppi2 values") - -#plot(density(my_df$averaged_duet) -plot(density(my_df$avg_ppi2_scaled) - , xlab = "" - , main = "mean ppi2 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 -#========= -# sanity check: match and assign to a separate column to double check -# colnames(my_df) -# df_duet$duet_scaled = my_df$averge_duet_scaled[match(df_duet$resno, my_df$position)] - -#========= -# step 1_P1 -#========= -# 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_duet$b = my_df$avg_ppi2_scaled[match(df_duet$resno, my_df$position)] - -#========= -# step 2_P1 -#========= -# count NA in Bfactor -b_na_duet = sum(is.na(df_duet$b)) ; b_na_duet - -# count number of 0"s in Bactor -sum(df_duet$b == 0) - -# replace all NA in b factor with 0 -na_rep = 2 -df_duet$b[is.na(df_duet$b)] = na_rep - -# # 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() -# } - -# 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 -#========= -# sanity check: dim should be same before reassignment -if ( (dim(df_duet)[1] == dim(d2_duet)[1]) & - (dim(df_duet)[2] == dim(d2_duet)[2]) - ){ - print("PASS: Dims of both dfs as expected") -} else { - print ("FAIL: Dims mismatch") - quit()} - -#========= -# step 4_P1: -# VERY important -#========= -# 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)) - -#========= -# step 5_P1 -#========= -cat(paste0("output file mean ppi2 pdb:" - , outfile_ppi2_mspdb)) -write.pdb(my_pdb_duet, outfile_ppi2_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 - , xlab = "" - , main = "repalcedB duet") - -plot(density(df_duet$b) - , xlab = "" - , main = "replacedB duet") - -# graph titles -mtext(text = "Frequency" - , side = 2 - , line = 0 - , outer = TRUE) - -mtext(text = paste0(tolower(gene), ": ppi2 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/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??? -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! \ No newline at end of file diff --git a/scripts/plotting/replaceBfactor_pdb_stability.R b/scripts/plotting/replaceBfactor_pdb_stability.R deleted file mode 100644 index e3be30f..0000000 --- a/scripts/plotting/replaceBfactor_pdb_stability.R +++ /dev/null @@ -1,281 +0,0 @@ -#!/usr/bin/env Rscript - -######################################################### -# TASK: Replace B-factors in the pdb file with the mean -# normalised stability values. - -# read pdb file - -# read mcsm mean stability value files -# extract the respective mean values and assign to the -# b-factor column within their respective pdbs - -# generate some distribution plots for inspection - -######################################################### -# working dir and loading libraries -getwd() -setwd("~/git/LSHTM_analysis/scripts/plotting") -cat(c(getwd(),"\n")) - -#source("~/git/LSHTM_analysis/scripts/Header_TT.R") -library(bio3d) -require("getopt", quietly = TRUE) # cmd parse arguments -#======================================================== -#drug = "pyrazinamide" -#gene = "pncA" - -# command line args -spec = matrix(c( - "drug" , "d", 1, "character", - "gene" , "g", 1, "character" -), byrow = TRUE, ncol = 4) - -opt = getopt(spec) - -drug = opt$drug -gene = opt$gene - -if(is.null(drug)|is.null(gene)) { - stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)") -} -#======================================================== -gene_match = paste0(gene,"_p.") -cat(gene_match) - -#============= -# directories -#============= -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/Writing/thesis/images/results/", tolower(gene)) - -#====== -# input -#====== -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_ens_stability.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_b_stab_ms.pdb") -outfile_duet_mspdb = paste0(outdir_plots, "/", out_filename_duet_mspdb) -print(paste0("Output file:", outfile_duet_mspdb)) - -#%%=============================================================== -#NOTE: duet here refers to the ensemble stability values - -########################### -# Read file: average stability values -# or mcsm_normalised file -########################### -my_df <- read.csv(infile_mean_stability, header = T) -str(my_df) - -############# -# Read pdb -############# -# list of 8 -my_pdb = read.pdb(infile_pdb - , maxlines = -1 - , multi = FALSE - , rm.insert = FALSE - , rm.alt = TRUE - , ATOM.only = FALSE - , hex = FALSE - , verbose = TRUE) - -rm(in_filename_mean_stability, in_filename_pdb) - -# assign separately for duet and ligand -my_pdb_duet = my_pdb - -#========================================================= -# Replacing B factor with mean stability scores -# within the respective dfs -#========================================================== -# extract atom list into a variable -# since in the list this corresponds to data frame, variable will be a df -#df_duet = my_pdb_duet[[1]] -df_duet= my_pdb_duet[['atom']] - -# make a copy: required for downstream sanity checks -d2_duet = df_duet - -# sanity checks: B factor -max(df_duet$b); min(df_duet$b) - -#================================================== -# histograms and density plots for inspection -# 1: original B-factors -# 2: original mean stability values -# 3: replaced B-factors with mean stability values -#================================================== -# Set the margin on all sides -par(oma = c(3,2,3,0) - , mar = c(1,3,5,2) - #, mfrow = c(3,2) - #, mfrow = c(3,4)) - , mfrow = c(3,2)) - - -#============= -# Row 1 plots: original B-factors -# duet and affinity -#============= -hist(df_duet$b - , xlab = "" - , main = "Bfactor stability") - -plot(density(df_duet$b) - , xlab = "" - , main = "Bfactor stability") - -#============= -# Row 2 plots: original mean stability values -# duet and affinity -#============= - -#hist(my_df$averaged_duet -hist(my_df$avg_ens_stability_scaled - , xlab = "" - , main = "mean stability values") - -#plot(density(my_df$averaged_duet) -plot(density(my_df$avg_ens_stability_scaled) - , xlab = "" - , main = "mean stability 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 -#========= -# sanity check: match and assign to a separate column to double check -# colnames(my_df) -# df_duet$duet_scaled = my_df$averge_duet_scaled[match(df_duet$resno, my_df$position)] - -#========= -# step 1_P1 -#========= -# 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_duet$b = my_df$avg_ens_stability_scaled[match(df_duet$resno, my_df$position)] - -#========= -# step 2_P1 -#========= -# count NA in Bfactor -b_na_duet = sum(is.na(df_duet$b)) ; b_na_duet - -# count number of 0"s in Bactor -sum(df_duet$b == 0) - -# replace all NA in b factor with 0 -na_rep = 2 -df_duet$b[is.na(df_duet$b)] = na_rep - -# # 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() -# } - -# 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 -#========= -# sanity check: dim should be same before reassignment -if ( (dim(df_duet)[1] == dim(d2_duet)[1]) & - (dim(df_duet)[2] == dim(d2_duet)[2]) - ){ - print("PASS: Dims of both dfs as expected") -} else { - print ("FAIL: Dims mismatch") - quit()} - -#========= -# step 4_P1: -# VERY important -#========= -# 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)) - -#========= -# step 5_P1 -#========= -cat(paste0("output file duet mean stability pdb:", outfile_duet_mspdb)) -write.pdb(my_pdb_duet, outfile_duet_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 - , xlab = "" - , main = "repalcedB duet") - -plot(density(df_duet$b) - , xlab = "" - , main = "replacedB duet") - -# graph titles -mtext(text = "Frequency" - , side = 2 - , line = 0 - , outer = TRUE) - -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/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??? -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! \ No newline at end of file diff --git a/scripts/plotting/replaceBfactor_pdb_stability_5uhc_rpob.R b/scripts/plotting/replaceBfactor_pdb_stability_5uhc_rpob.R deleted file mode 100644 index d263820..0000000 --- a/scripts/plotting/replaceBfactor_pdb_stability_5uhc_rpob.R +++ /dev/null @@ -1,280 +0,0 @@ -#!/usr/bin/env Rscript - -######################################################### -# TASK: Replace B-factors in the pdb file with the mean -# normalised stability values. - -# read pdb file - -# read mcsm mean stability value files -# extract the respective mean values and assign to the -# b-factor column within their respective pdbs - -# generate some distribution plots for inspection - -######################################################### -# working dir and loading libraries -getwd() -setwd("~/git/LSHTM_analysis/scripts/plotting") -cat(c(getwd(),"\n")) - -#source("~/git/LSHTM_analysis/scripts/Header_TT.R") -library(bio3d) -require("getopt", quietly = TRUE) # cmd parse arguments -#======================================================== -#drug = "pyrazinamide" -#gene = "pncA" - -# command line args -spec = matrix(c( - "drug" , "d", 1, "character", - "gene" , "g", 1, "character" -), byrow = TRUE, ncol = 4) - -opt = getopt(spec) - -drug = opt$drug -gene = opt$gene - -if(is.null(drug)|is.null(gene)) { - stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)") -} -#======================================================== -gene_match = paste0(gene,"_p.") -cat(gene_match) - -#============= -# directories -#============= -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/Writing/thesis/images/results/", tolower(gene)) - -#====== -# input -#====== -#in_filename_pdb = paste0(tolower(gene), "_complex.pdb") -in_filename_pdb = "5uhc.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("/5uhc_", tolower(gene), "_mean_ens_stability.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("/5uhc_", tolower(gene), "_complex_b_stab_ms.pdb") -outfile_duet_mspdb = paste0(outdir_plots, out_filename_duet_mspdb) -print(paste0("Output file:", outfile_duet_mspdb)) - -#%%=============================================================== -#NOTE: duet here refers to the ensemble stability values - -########################### -# Read file: average stability values -# or mcsm_normalised file -########################### -my_df <- read.csv(infile_mean_stability, header = T) -str(my_df) -my_df = na.omit(my_df) - -############# -# Read pdb -############# -# list of 8 -my_pdb = read.pdb(infile_pdb - , maxlines = -1 - , multi = FALSE - , rm.insert = FALSE - , rm.alt = TRUE - , ATOM.only = FALSE - , hex = FALSE - , verbose = TRUE) - -rm(in_filename_mean_stability, in_filename_pdb) - -# assign separately for duet and ligand -my_pdb_duet = my_pdb - -#========================================================= -# Replacing B factor with mean stability scores -# within the respective dfs -#========================================================== -# extract atom list into a variable -# since in the list this corresponds to data frame, variable will be a df -#df_duet = my_pdb_duet[[1]] -df_duet= my_pdb_duet[['atom']] - -# make a copy: required for downstream sanity checks -d2_duet = df_duet - -# sanity checks: B factor -max(df_duet$b); min(df_duet$b) - -#================================================== -# histograms and density plots for inspection -# 1: original B-factors -# 2: original mean stability values -# 3: replaced B-factors with mean stability values -#================================================== -# Set the margin on all sides -par(oma = c(3,2,3,0) - , mar = c(1,3,5,2) - #, mfrow = c(3,2) - #, mfrow = c(3,4)) - , mfrow = c(3,2)) -#============= -# Row 1 plots: original B-factors -# duet and affinity -#============= -hist(df_duet$b - , xlab = "" - , main = "Bfactor stability") - -plot(density(df_duet$b) - , xlab = "" - , main = "Bfactor stability") - -#============= -# Row 2 plots: original mean stability values -# duet and affinity -#============= - -#hist(my_df$averaged_duet -hist(my_df$avg_ens_stability_scaled - , xlab = "" - , main = "mean stability values") - -#plot(density(my_df$averaged_duet) -plot(density(my_df$avg_ens_stability_scaled) - , xlab = "" - , main = "mean stability 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 -#========= -# sanity check: match and assign to a separate column to double check -# colnames(my_df) -# df_duet$duet_scaled = my_df$averge_duet_scaled[match(df_duet$resno, my_df$position)] - -#========= -# step 1_P1 -#========= -# 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_duet$b = my_df$avg_ens_stability_scaled[match(df_duet$resno, my_df$X5uhc_position)] - -#========= -# step 2_P1 -#========= -# count NA in Bfactor -b_na_duet = sum(is.na(df_duet$b)) ; b_na_duet - -# count number of 0"s in Bactor -sum(df_duet$b == 0) - -# replace all NA in b factor with 0 -na_rep = 2 -df_duet$b[is.na(df_duet$b)] = na_rep - -# # 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() -# } - -# 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 -#========= -# sanity check: dim should be same before reassignment -if ( (dim(df_duet)[1] == dim(d2_duet)[1]) & - (dim(df_duet)[2] == dim(d2_duet)[2]) - ){ - print("PASS: Dims of both dfs as expected") -} else { - print ("FAIL: Dims mismatch") - quit()} - -#========= -# step 4_P1: -# VERY important -#========= -# 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)) - -#========= -# step 5_P1 -#========= -cat(paste0("output file duet mean stability pdb:", outfile_duet_mspdb)) -write.pdb(my_pdb_duet, outfile_duet_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 - , xlab = "" - , main = "repalcedB duet") - -plot(density(df_duet$b) - , xlab = "" - , main = "replacedB duet") - -# graph titles -mtext(text = "Frequency" - , side = 2 - , line = 0 - , outer = TRUE) - -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/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??? -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! \ No newline at end of file