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