From 23b4f06017da58ea63fdaa4c704d30fee084127b Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Tue, 23 Aug 2022 16:30:42 +0100 Subject: [PATCH] added scripts --- config/gid.R | 6 +- .../gid/basic_barplots_gid_layout.R | 8 +- .../gid/sensitivity_count_gid.R | 7 +- .../AFFINITY_TEST_prominent_effects.R | 138 -------- .../mcsm_mean_affinity_ensemble.R | 316 ------------------ .../mcsm_mean_stability_ensemble.R | 74 ---- .../mcsm_mean_stability_ensemble_5uhc_rpob.R | 176 ---------- .../replaceBfactor_pdb_affinity.R | 152 +++------ .../replaceBfactor_pdb_ppi2.R | 158 +++------ .../replaceBfactor_pdb_stability.R | 126 +++---- 10 files changed, 147 insertions(+), 1014 deletions(-) delete mode 100644 scripts/plotting/structure_figures/AFFINITY_TEST_prominent_effects.R delete mode 100644 scripts/plotting/structure_figures/mcsm_mean_affinity_ensemble.R delete mode 100644 scripts/plotting/structure_figures/mcsm_mean_stability_ensemble.R delete mode 100644 scripts/plotting/structure_figures/mcsm_mean_stability_ensemble_5uhc_rpob.R diff --git a/config/gid.R b/config/gid.R index 353607c..763a429 100644 --- a/config/gid.R +++ b/config/gid.R @@ -135,5 +135,9 @@ aa_pos_lig2 = aa_pos_rna aa_pos_lig3 = aa_pos_amp tile_map=data.frame(tile=c("SRY","SAM","RNA","AMP"), - tile_colour=c("green","darkslategrey","navyblue","purple")) + tile_colour=c("green","darkslategrey","darkred","navyblue")) +# green: #00ff00 +# darkslategrey : #2f4f4f +# darkred : #8b0000 +# navyblue :#000080 \ No newline at end of file diff --git a/scripts/plotting/plotting_thesis/gid/basic_barplots_gid_layout.R b/scripts/plotting/plotting_thesis/gid/basic_barplots_gid_layout.R index 11b1c5b..c8429ba 100644 --- a/scripts/plotting/plotting_thesis/gid/basic_barplots_gid_layout.R +++ b/scripts/plotting/plotting_thesis/gid/basic_barplots_gid_layout.R @@ -2,12 +2,18 @@ # Data: Input #============== source("~/git/LSHTM_analysis/config/gid.R") -source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") + source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") #cat("\nSourced plotting cols as well:", length(plotting_cols)) source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/gid/basic_barplots_gid.R") source("/home/tanu/git/LSHTM_analysis/scripts/plotting/plotting_thesis/gid/pe_sens_site_count_gid.R") +#======= +# output +#======= +outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/") +cat("plots will output to:", outdir_images) +######################################################### if ( tolower(gene)%in%c("gid") ){ cat("\nPlots available for layout are:") diff --git a/scripts/plotting/plotting_thesis/gid/sensitivity_count_gid.R b/scripts/plotting/plotting_thesis/gid/sensitivity_count_gid.R index 66caef2..b11eedf 100644 --- a/scripts/plotting/plotting_thesis/gid/sensitivity_count_gid.R +++ b/scripts/plotting/plotting_thesis/gid/sensitivity_count_gid.R @@ -12,11 +12,12 @@ sensP_df = merged_df3[,c("mutationinformation", head(sensP_df) table(sensP_df$sensitivity) -#--------------- +#-------------------------- # Total unique positions -#---------------- +#-------------------------- tot_mut_pos = length(unique(sensP_df[[pos_colname_c]])) cat("\nNo of Tot muts sites:", tot_mut_pos) +cat("\nThese are:", unique(sensP_df[[pos_colname_c]])) # resistant mut pos sens_site_allR = sensP_df[[pos_colname_c]][sensP_df$sensitivity=="R"] @@ -34,7 +35,7 @@ length(sens_site_UR) common_pos = intersect(sens_site_UR,sens_site_US) site_Cc = length(common_pos) cat("\nNo of Common sites:", site_Cc - , "\nThese are:", common_pos) + , "\nThese are:", sort(unique(common_pos))) #--------------- # Resistant muts diff --git a/scripts/plotting/structure_figures/AFFINITY_TEST_prominent_effects.R b/scripts/plotting/structure_figures/AFFINITY_TEST_prominent_effects.R deleted file mode 100644 index 2ac94b8..0000000 --- a/scripts/plotting/structure_figures/AFFINITY_TEST_prominent_effects.R +++ /dev/null @@ -1,138 +0,0 @@ - -foo = df3_affinity_filtered[df3_affinity_filtered$ligand_distance<10,] -bar = df3_affinity_filtered[df3_affinity_filtered$interface_dist<10,] - -wilcox.test(foo$mmcsm_lig_scaled~foo$sensitivity) -wilcox.test(foo$mmcsm_lig~foo$sensitivity) - -wilcox.test(foo$affinity_scaled~foo$sensitivity) -wilcox.test(foo$ligand_affinity_change~foo$sensitivity) - -wilcox.test(bar$mcsm_na_scaled~bar$sensitivity) -wilcox.test(bar$mcsm_na_affinity~bar$sensitivity) - -wilcox.test(bar$mcsm_ppi2_scaled~bar$sensitivity) -wilcox.test(bar$mcsm_ppi2_affinity~bar$sensitivity) - - -# find the most "impactful" effect value -biggest=max(abs((a[c('affinity_scaled','mmcsm_lig_scaled','mcsm_ppi2_scaled','mcsm_na_scaled')]))) - -abs((a[c('affinity_scaled','mmcsm_lig_scaled','mcsm_ppi2_scaled','mcsm_na_scaled')]))==biggest - -abs((a[c('affinity_scaled','mmcsm_lig_scaled','mcsm_ppi2_scaled','mcsm_na_scaled')]))==c(,biggest) - -max(abs((a[c('affinity_scaled','mmcsm_lig_scaled','mcsm_ppi2_scaled','mcsm_na_scaled')]))) - - -a2 = (a[c('affinity_scaled','mmcsm_lig_scaled','mcsm_ppi2_scaled','mcsm_na_scaled')]) -a2 -# -# biggest = max(abs(a2[1,])) -# -# #hmm -# #which(abs(a2) == biggest) -# #names(a2)[apply(a2, 1:4, function(i) which(i == max()))] -# -# # get row max -# a2$row_maximum = apply(abs(a2[,-1]), 1, max) -# -# # get colname for abs(max_value) -# #https://stackoverflow.com/questions/36960010/get-column-name-that-matches-specific-row-value-in-dataframe -# #names(df)[which(df == 1, arr.ind=T)[, "col"]] -# # yayy -# names(a2)[which(abs(a2) == biggest, arr.ind=T)[, "col"]] -# -# #another:https://statisticsglobe.com/return-column-name-of-largest-value-for-each-row-in-r -# colnames(a2)[max.col(abs(a2), ties.method = "first")] # Apply colnames & max.col functions -# ################################################# -# # use whole df -# #gene_aff_cols = c('affinity_scaled','mmcsm_lig_scaled','mcsm_ppi2_scaled','mcsm_na_scaled') -# -# biggest = max(abs(a[gene_aff_cols])) -# a$max_es = biggest -# a$effect = names(a[gene_aff_cols])[which(abs(a[gene_aff_cols]) == biggest, arr.ind=T)[, "col"]] -# -# effect_name = unique(a$effect) -# #get index of value of max effect -# ind = (which(abs(a[effect_name]) == biggest, arr.ind=T)) -# a[effect_name][ind] -# # extract sign -# a$effect_sign = sign(a[effect_name][ind]) -######################################################## -# maxn <- function(n) function(x) order(x, decreasing = TRUE)[n] -# second_big = abs(a[gene_aff_cols])[maxn(2)(abs(a[gene_aff_cols])] -# apply(df, 1, function(x)x[maxn(1)(x)]) -# apply(a[gene_aff_cols], 1, function(x) abs(a[gene_aff_cols])[maxn(2)(abs(a[gene_aff_cols]))]) -######################################################### -# loop -a2 = df2[df2$position%in%c(167, 423, 427),] -test <- a2 %>% - dplyr::group_by(position) %>% - biggest = max(abs(a2[gene_aff_cols])) - a2$max_es = max(abs(a2[gene_aff_cols])) - a2$effect = names(a2[gene_aff_cols])[which(abs(a2[gene_aff_cols]) == biggest, arr.ind=T)[, "col"]] - effect_name = unique(a2$effect) - - #get index of value of max effect - ind = (which(abs(a2[effect_name]) == biggest, arr.ind=T)) - a2[effect_name][ind] - # extract sign - a2$effect_dir = sign(a2[effect_name][ind]) -################################# -df2_short = df2[df2$position%in%c(167, 423, 427),] - -for (i in unique(df2_short$position) ){ - #print(i) - #print(paste0("\nNo. of unique positions:", length(unique(df2$position))) ) - #cat(length(unique(df2$position))) - a2 = df2_short[df2_short$position==i,] - biggest = max(abs(a2[gene_aff_cols])) - a2$max_es = max(abs(a2[gene_aff_cols])) - a2$effect = names(a2[gene_aff_cols])[which(abs(a2[gene_aff_cols]) == biggest, arr.ind=T)[, "col"]] - effect_name = unique(a2$effect) - - #get index of value of max effect - ind = (which(abs(a2[effect_name]) == biggest, arr.ind=T)) - a2[effect_name][ind] - # extract sign - a2$effect_sign = sign(a2[effect_name][ind]) -} - -#======================== -df2_short = df3[df3$position%in%c(167, 423, 427),] -df2_short = df3[df3$position%in%c(170, 167, 493, 453, 435, 433, 480, 456, 445),] -df2_short = df3[df3$position%in%c(435, 480),] -df2_short = df3[df3$position%in%c(435, 480),] - -give_col=function(x,y,df=df2_short){ - df[df$position==x,y] -} - -for (i in unique(df2_short$position) ){ - #print(i) - #print(paste0("\nNo. of unique positions:", length(unique(df2$position))) ) - #cat(length(unique(df2$position))) - #df2_short[df2_short$position==i,gene_aff_cols] - - biggest = max(abs(give_col(i,gene_aff_cols))) - - df2_short[df2_short$position==i,'abs_max_effect'] = biggest - df2_short[df2_short$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 = df2_short[df2_short$position==i,'effect_type'][1] # pick first one in case we have multiple exact values - - # get index/rowname for value of max effect, and then use it to get the original sign - # here - #df2_short[df2_short$position==i,c(effect_name)] - #which(abs(df2_short[df2_short$position==i,c('position',effect_name)][effect_name])==biggest, arr.ind=T) - ind = rownames(which(abs(df2_short[df2_short$position==i,c('position',effect_name)][effect_name])== biggest, arr.ind=T)) - df2_short[df2_short$position==i,'effect_sign'] = sign(df2_short[effect_name][ind,]) -} - -df2_short$effect_type = sub("\\.[0-9]+", "", df2_short$effect_type) # cull duplicate effect types that happen when there are exact duplicate values \ No newline at end of file diff --git a/scripts/plotting/structure_figures/mcsm_mean_affinity_ensemble.R b/scripts/plotting/structure_figures/mcsm_mean_affinity_ensemble.R deleted file mode 100644 index 068325b..0000000 --- a/scripts/plotting/structure_figures/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/structure_figures/mcsm_mean_stability_ensemble.R b/scripts/plotting/structure_figures/mcsm_mean_stability_ensemble.R deleted file mode 100644 index 26beac4..0000000 --- a/scripts/plotting/structure_figures/mcsm_mean_stability_ensemble.R +++ /dev/null @@ -1,74 +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") - -######################################################### -# TASK: Generate averaged stability values by position -# calculated 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: merged_df3 -#============= -source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") -#merged_df3= paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv") - -cols_to_extract_ms = c("mutationinformation", "position", "avg_stability_scaled") - -df3 = merged_df3[, cols_to_extract_ms] -length(df3$mutationinformation) - -# ensemble average of predictors by position -avg_stability_by_position <- df3 %>% - dplyr::group_by(position) %>% - dplyr::summarize(avg_stability_scaled_pos = mean(avg_stability_scaled)) - -min(avg_stability_by_position$avg_stability_scaled_pos) -max(avg_stability_by_position$avg_stability_scaled_pos) - -avg_stability_by_position['avg_stability_scaled_pos_scaled'] = lapply(avg_stability_by_position['avg_stability_scaled_pos'] - , function(x) { - scales::rescale_mid(x, to = c(-1,1) - #, from = c(en_stab_min,en_stab_max)) - , mid = 0 - , from = c(0,1)) - }) -cat(paste0('Average stability scores:\n' - , head(avg_stability_by_position['avg_stability_scaled_pos']) - , '\n---------------------------------------------------------------' - , '\nAverage stability scaled scores:\n' - , head(avg_stability_by_position['avg_stability_scaled_pos_scaled']) - )) - -all(avg_stability_by_position['avg_stability_scaled_pos'] == avg_stability_by_position['avg_stability_scaled_pos_scaled']) - -# convert to a data frame -avg_stability_by_position = as.data.frame(avg_stability_by_position) - -################################################################## -# output -#write.csv(combined_df, outfile_mean_ens_st_aff -write.csv(avg_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(avg_stability_by_position) - , "\nNo. of cols:", ncol(avg_stability_by_position)) - -# end of script -#=============================================================== diff --git a/scripts/plotting/structure_figures/mcsm_mean_stability_ensemble_5uhc_rpob.R b/scripts/plotting/structure_figures/mcsm_mean_stability_ensemble_5uhc_rpob.R deleted file mode 100644 index c75a6ba..0000000 --- a/scripts/plotting/structure_figures/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/structure_figures/replaceBfactor_pdb_affinity.R b/scripts/plotting/structure_figures/replaceBfactor_pdb_affinity.R index 1ba5921..4775e1f 100644 --- a/scripts/plotting/structure_figures/replaceBfactor_pdb_affinity.R +++ b/scripts/plotting/structure_figures/replaceBfactor_pdb_affinity.R @@ -1,46 +1,23 @@ #!/usr/bin/env Rscript +source("~/git/LSHTM_analysis/config/gid.R") + +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") ######################################################### # 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 +# normalised affinity values ######################################################### # 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) +gene_match = paste0(gene,"_p."); cat(gene_match) cat(gene_match) #============= @@ -49,9 +26,13 @@ cat(gene_match) datadir = paste0("~/git/Data") indir = paste0(datadir, "/", drug, "/input") outdir = paste0("~/git/Data", "/", drug, "/output") -#outdir_plots = paste0("~/git/Data", "/", drug, "/output/plots") -outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) +#outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) +#======= +# output +#======= +outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/") +cat("plots will output to:", outdir_images) #====== # input #====== @@ -59,31 +40,31 @@ 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) +outfile_lig_mspdb = paste0(outdir_images,out_filename_lig_mspdb) print(paste0("Output file:", outfile_lig_mspdb)) #%%=============================================================== -#NOTE: duet here refers to the ensemble stability values +#NOTE: duet here refers to the ensemble affinity values ########################### -# Read file: average stability values +# Read file: average affinity values # or mcsm_normalised file ########################### -my_df <- read.csv(infile_mean_stability, header = T) -str(my_df) +my_df_raw = merged_df3[, c("position", "ligand_distance", "avg_lig_affinity_scaled", "avg_lig_affinity")] +my_df_raw = my_df_raw[my_df_raw$ligand_distance<10,] + +# avg by position on the SCALED values +my_df <- my_df_raw %>% + group_by(position) %>% + summarize(avg_ligaff_sc_pos = mean(avg_lig_affinity_scaled)) + +max(my_df$avg_ligaff_sc_pos) +min(my_df$avg_ligaff_sc_pos) ############# # Read pdb @@ -98,13 +79,11 @@ my_pdb = read.pdb(infile_pdb , 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 +# Replacing B factor with mean affinity scores # within the respective dfs #========================================================== # extract atom list into a variable @@ -121,8 +100,8 @@ 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 +# 2: original mean affinity values +# 3: replaced B-factors with mean affinity values #================================================== # Set the margin on all sides par(oma = c(3,2,3,0) @@ -131,6 +110,7 @@ par(oma = c(3,2,3,0) #, mfrow = c(3,4)) , mfrow = c(3,2)) + #============= # Row 1 plots: original B-factors # duet and affinity @@ -144,40 +124,28 @@ plot(density(df_duet$b) , main = "Bfactor affinity") #============= -# Row 2 plots: original mean stability values -# duet and affinity +# Row 2 plots: original mean affinity values +# affinity #============= #hist(my_df$averaged_duet -hist(my_df$avg_lig_scaled +hist(my_df$avg_ligaff_sc_pos , xlab = "" , main = "mean affinity values") #plot(density(my_df$averaged_duet) -plot(density(my_df$avg_lig_scaled) +plot(density(my_df$avg_ligaff_sc_pos) , 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)] +df_duet$b = my_df$avg_ligaff_sc_pos[match(df_duet$resno, my_df$position)] #========= # step 2_P1 @@ -192,32 +160,6 @@ sum(df_duet$b == 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 #========= @@ -241,17 +183,23 @@ table(df_duet$b) sum(is.na(df_duet$b)) #========= -# step 5_P1 +# step 5_P1: OUTPUT #========= -cat(paste0("output file duet mean stability pdb:" - , outfile_lig_mspdb)) +cat(paste0("output file duet mean affinity pdb:", outfile_lig_mspdb)) write.pdb(my_pdb_duet, outfile_lig_mspdb) +# OUTPUT: position file +poscsvF = paste0(outdir_images, tolower(gene), "_ligaff_positions.csv") +cat(paste0("output file duet mean NA affinity POSITIONS:", poscsvF)) + +filtered_pos = toString(my_df$position) +write.table(filtered_pos, poscsvF, row.names = F, col.names = F ) + #============================ # Add the 3rd histogram and density plots for comparisons #============================ # Plots continued... -# Row 3 plots: hist and density of replaced B-factors with stability values +# Row 3 plots: hist and density of replaced B-factors with affinity values hist(df_duet$b , xlab = "" , main = "repalcedB duet") @@ -266,16 +214,8 @@ mtext(text = "Frequency" , line = 0 , outer = TRUE) -mtext(text = paste0(tolower(gene), ": afinity distribution") +mtext(text = paste0(tolower(gene), ": affinity 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/structure_figures/replaceBfactor_pdb_ppi2.R b/scripts/plotting/structure_figures/replaceBfactor_pdb_ppi2.R index 03efa38..55ed85a 100644 --- a/scripts/plotting/structure_figures/replaceBfactor_pdb_ppi2.R +++ b/scripts/plotting/structure_figures/replaceBfactor_pdb_ppi2.R @@ -1,46 +1,19 @@ #!/usr/bin/env Rscript +source("~/git/LSHTM_analysis/config/gid.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") ######################################################### # 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 - +# normalised ppi2 values. ######################################################### -# 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) +gene_match = paste0(gene,"_p."); cat(gene_match) cat(gene_match) #============= @@ -49,9 +22,13 @@ cat(gene_match) datadir = paste0("~/git/Data") indir = paste0(datadir, "/", drug, "/input") outdir = paste0("~/git/Data", "/", drug, "/output") -#outdir_plots = paste0("~/git/Data", "/", drug, "/output/plots") -outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) +#outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) +#======= +# output +#======= +outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/") +cat("plots will output to:", outdir_images) #====== # input #====== @@ -59,30 +36,33 @@ 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) +outfile_ppi2_mspdb = paste0(outdir_images,out_filename_ppi2_mspdb) print(paste0("Output file:", outfile_ppi2_mspdb)) #%%=============================================================== -#NOTE: duet here refers to the ensemble stability values +#NOTE: duet here refers to the ensemble ppi2 values ########################### -# Read file: average stability values +# Read file: average ppi2 values # or mcsm_normalised file ########################### -my_df <- read.csv(infile_mean_ppi2, header = T) -str(my_df) +my_df_raw = merged_df3[, c("position", "mcsm_ppi2_scaled", "interface_dist")] +head(my_df_raw) +my_df_raw = my_df_raw[my_df_raw$interface_dist<10,] +my_df_raw$position +# avg by position on the SCALED values +my_df <- my_df_raw %>% + group_by(position) %>% + summarize(avg_ppi2_sc_pos = mean(mcsm_ppi2_scaled)) + +max(my_df$avg_ppi2_sc_pos) +min(my_df$avg_ppi2_sc_pos) +#============================================================ ############# # Read pdb ############# @@ -100,7 +80,7 @@ my_pdb = read.pdb(infile_pdb my_pdb_duet = my_pdb #========================================================= -# Replacing B factor with mean stability scores +# Replacing B factor with mean ppi2 scores # within the respective dfs #========================================================== # extract atom list into a variable @@ -117,8 +97,8 @@ 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 +# 2: original mean ppi2 values +# 3: replaced B-factors with mean ppi2 values #================================================== # Set the margin on all sides par(oma = c(3,2,3,0) @@ -129,7 +109,7 @@ par(oma = c(3,2,3,0) #============= # Row 1 plots: original B-factors -# duet and affinity +# duet and ppi2 #============= hist(df_duet$b , xlab = "" @@ -140,40 +120,24 @@ plot(density(df_duet$b) , main = "Bfactor ppi2") #============= -# Row 2 plots: original mean stability values -# duet and affinity +# Row 2 plots: original mean ppi2 values +# ppi2 #============= #hist(my_df$averaged_duet -hist(my_df$avg_ppi2_scaled +hist(my_df$avg_ppi2_sc_pos , xlab = "" , main = "mean ppi2 values") #plot(density(my_df$averaged_duet) -plot(density(my_df$avg_ppi2_scaled) +plot(density(my_df$avg_ppi2_sc_pos) , 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)] +df_duet$b = my_df$avg_ppi2_sc_pos[match(df_duet$resno, my_df$position)] #========= # step 2_P1 @@ -188,32 +152,6 @@ sum(df_duet$b == 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 #========= @@ -237,17 +175,23 @@ table(df_duet$b) sum(is.na(df_duet$b)) #========= -# step 5_P1 +# step 5_P1: OUTPUT #========= -cat(paste0("output file mean ppi2 pdb:" - , outfile_ppi2_mspdb)) +cat(paste0("output file duet mean ppi2 pdb:", outfile_ppi2_mspdb)) write.pdb(my_pdb_duet, outfile_ppi2_mspdb) +# OUTPUT: position file +poscsvF = paste0(outdir_images, tolower(gene), "_ppi2_positions.csv") +cat(paste0("output file duet mean ppi2 POSITIONS:", poscsvF)) + +filtered_pos = toString(my_df$position) +write.table(filtered_pos, poscsvF, row.names = F, col.names = F ) + #============================ # Add the 3rd histogram and density plots for comparisons #============================ # Plots continued... -# Row 3 plots: hist and density of replaced B-factors with stability values +# Row 3 plots: hist and density of replaced B-factors with ppi2 values hist(df_duet$b , xlab = "" , main = "repalcedB duet") @@ -266,12 +210,4 @@ 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 +#============================================ \ No newline at end of file diff --git a/scripts/plotting/structure_figures/replaceBfactor_pdb_stability.R b/scripts/plotting/structure_figures/replaceBfactor_pdb_stability.R index fa85162..ba6999b 100644 --- a/scripts/plotting/structure_figures/replaceBfactor_pdb_stability.R +++ b/scripts/plotting/structure_figures/replaceBfactor_pdb_stability.R @@ -1,11 +1,7 @@ #!/usr/bin/env Rscript - -#source("~/git/LSHTM_analysis/config/alr.R") -source("~/git/LSHTM_analysis/config/embb.R") -#source("~/git/LSHTM_analysis/config/katg.R") -#source("~/git/LSHTM_analysis/config/gid.R") -#source("~/git/LSHTM_analysis/config/pnca.R") -#source("~/git/LSHTM_analysis/config/rpob.R") +source("~/git/LSHTM_analysis/config/gid.R") +source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R") + ######################################################### # TASK: Replace B-factors in the pdb file with the mean # normalised stability values. @@ -20,31 +16,12 @@ source("~/git/LSHTM_analysis/config/embb.R") ######################################################### # 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 = "" -#gene = "" -# # 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)") -# } #======================================================== cat(gene) gene_match = paste0(gene,"_p."); cat(gene_match) @@ -56,29 +33,25 @@ cat(gene_match) datadir = paste0("~/git/Data") indir = paste0(datadir, "/", drug, "/input") outdir = paste0("~/git/Data", "/", drug, "/output") -#outdir_plots = paste0("~/git/Data", "/", drug, "/output/plots") -outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) +#outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene)) #====== # input #====== in_filename_pdb = paste0(tolower(gene), "_complex.pdb") +#in_filename_pdb = "/home/tanu/git/Writing/thesis/images/results/gid/str_figures/gid_complex_copy_arpeg.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 #======= +outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/") +cat("plots will output to:", outdir_images) + #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) +outfile_duet_mspdb = paste0(outdir_images, out_filename_duet_mspdb) print(paste0("Output file:", outfile_duet_mspdb)) #%%=============================================================== @@ -88,8 +61,31 @@ print(paste0("Output file:", outfile_duet_mspdb)) # Read file: average stability values # or mcsm_normalised file ########################### -my_df <- read.csv(infile_mean_stability, header = T) -str(my_df) +my_df_raw = merged_df3[, c("position", "avg_stability", "avg_stability_scaled")] + +# avg by position on the SCALED values +my_df <- my_df_raw %>% + group_by(position) %>% + summarize(avg_stab_sc_pos = mean(avg_stability_scaled)) + +max(my_df$avg_stab_sc_pos) +min(my_df$avg_stab_sc_pos) +#============================================================ +# # scale b/w -1 and 1 +# duet_min = min(my_df_by_position['avg_stab_sc_pos']) +# duet_max = max(my_df_by_position['avg_stab_sc_pos']) +# +# # scale the averaged_duet values +# my_df_by_position['avg_stab_sc_pos_scaled'] = lapply(my_df_by_position['avg_stab_sc_pos'] +# , function(x) ifelse(x < 0, x/abs(duet_min), x/duet_max)) +# +# cat(paste0('Average duet scores:\n', head(my_df_by_position['avg_stab_sc_pos_scaled']) +# , '\n---------------------------------------------------------------' +# , '\nScaled duet scores:\n', head(my_df_by_position['avg_stab_sc_pos_scaled']))) +# +# min(my_df_by_position['avg_stab_sc_pos_scaled']) +# max(my_df_by_position['avg_stab_sc_pos_scaled']) +#============================================================ ############# # Read pdb @@ -104,8 +100,6 @@ my_pdb = read.pdb(infile_pdb , hex = FALSE , verbose = TRUE) -rm(in_filename_mean_stability, in_filename_pdb) - # assign separately for duet and ligand my_pdb_duet = my_pdb @@ -113,9 +107,6 @@ 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 @@ -156,35 +147,22 @@ plot(density(df_duet$b) #============= #hist(my_df$averaged_duet -hist(my_df$avg_stability_scaled_pos_scaled +hist(my_df$avg_stab_sc_pos , xlab = "" , main = "mean stability values") #plot(density(my_df$averaged_duet) -plot(density(my_df$avg_stability_scaled_pos_scaled) +plot(density(my_df$avg_stab_sc_pos) , 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_stability_scaled_pos_scaled[match(df_duet$resno, my_df$position)] +df_duet$b = my_df$avg_stab_sc_pos[match(df_duet$resno, my_df$position)] #========= # step 2_P1 @@ -198,26 +176,6 @@ 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 -# if ( (sum(df_duet$b == na_rep) == b_na_duet) { -# 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_stability_scaled_pos_scaled)) & (min(df_duet$b) == min(my_df$avg_stability_scaled_pos_scaled)) ){ -# print("PASS: B-factors replaced correctly in df_duet") -# } else { -# print ("FAIL: To replace B-factors in df_duet") -# quit() -# } - #========= # step 3_P1 #========= @@ -269,12 +227,4 @@ 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 +#============================================ \ No newline at end of file