diff --git a/scripts/plotting/structure_figures/AFFINITY_TEST.R b/scripts/plotting/structure_figures/AFFINITY_TEST.R new file mode 100644 index 0000000..08a1331 --- /dev/null +++ b/scripts/plotting/structure_figures/AFFINITY_TEST.R @@ -0,0 +1,138 @@ + +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_affinity_data_only.R b/scripts/plotting/structure_figures/mcsm_affinity_data_only.R new file mode 100644 index 0000000..f962935 --- /dev/null +++ b/scripts/plotting/structure_figures/mcsm_affinity_data_only.R @@ -0,0 +1,241 @@ +#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/structure_figures/mcsm_mean_affinity_ensemble.R b/scripts/plotting/structure_figures/mcsm_mean_affinity_ensemble.R new file mode 100644 index 0000000..068325b --- /dev/null +++ b/scripts/plotting/structure_figures/mcsm_mean_affinity_ensemble.R @@ -0,0 +1,316 @@ +#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.R b/scripts/plotting/structure_figures/mcsm_mean_stability.R new file mode 100755 index 0000000..7e35e79 --- /dev/null +++ b/scripts/plotting/structure_figures/mcsm_mean_stability.R @@ -0,0 +1,163 @@ +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/structure_figures/mcsm_mean_stability_ensemble.R b/scripts/plotting/structure_figures/mcsm_mean_stability_ensemble.R new file mode 100644 index 0000000..aeb5d0a --- /dev/null +++ b/scripts/plotting/structure_figures/mcsm_mean_stability_ensemble.R @@ -0,0 +1,212 @@ +#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/structure_figures/mcsm_mean_stability_ensemble_5uhc_rpob.R b/scripts/plotting/structure_figures/mcsm_mean_stability_ensemble_5uhc_rpob.R new file mode 100644 index 0000000..c75a6ba --- /dev/null +++ b/scripts/plotting/structure_figures/mcsm_mean_stability_ensemble_5uhc_rpob.R @@ -0,0 +1,176 @@ +#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/mut_landscape.R b/scripts/plotting/structure_figures/mut_landscape.R new file mode 100644 index 0000000..aeeb5c2 --- /dev/null +++ b/scripts/plotting/structure_figures/mut_landscape.R @@ -0,0 +1,155 @@ +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/structure_figures/mut_landscape_5uhc_rpob.R b/scripts/plotting/structure_figures/mut_landscape_5uhc_rpob.R new file mode 100644 index 0000000..53f95d0 --- /dev/null +++ b/scripts/plotting/structure_figures/mut_landscape_5uhc_rpob.R @@ -0,0 +1,191 @@ +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/structure_figures/replaceBfactor_pdb.R b/scripts/plotting/structure_figures/replaceBfactor_pdb.R new file mode 100755 index 0000000..a840c19 --- /dev/null +++ b/scripts/plotting/structure_figures/replaceBfactor_pdb.R @@ -0,0 +1,332 @@ +#!/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/structure_figures/replaceBfactor_pdb_affinity.R b/scripts/plotting/structure_figures/replaceBfactor_pdb_affinity.R new file mode 100644 index 0000000..1ba5921 --- /dev/null +++ b/scripts/plotting/structure_figures/replaceBfactor_pdb_affinity.R @@ -0,0 +1,281 @@ +#!/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/structure_figures/replaceBfactor_pdb_ppi2.R b/scripts/plotting/structure_figures/replaceBfactor_pdb_ppi2.R new file mode 100644 index 0000000..03efa38 --- /dev/null +++ b/scripts/plotting/structure_figures/replaceBfactor_pdb_ppi2.R @@ -0,0 +1,277 @@ +#!/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/structure_figures/replaceBfactor_pdb_stability.R b/scripts/plotting/structure_figures/replaceBfactor_pdb_stability.R new file mode 100644 index 0000000..e3be30f --- /dev/null +++ b/scripts/plotting/structure_figures/replaceBfactor_pdb_stability.R @@ -0,0 +1,281 @@ +#!/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/structure_figures/replaceBfactor_pdb_stability_5uhc_rpob.R b/scripts/plotting/structure_figures/replaceBfactor_pdb_stability_5uhc_rpob.R new file mode 100644 index 0000000..d263820 --- /dev/null +++ b/scripts/plotting/structure_figures/replaceBfactor_pdb_stability_5uhc_rpob.R @@ -0,0 +1,280 @@ +#!/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