tanu/LSHTM_analysis
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

added scripts to generate mean stability for rpob

Browse source
This commit is contained in:
Tanushree Tunstall 2022-08-01 21:41:55 +01:00
parent ccc877e811
commit 66337c289c
3 changed files with 702 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

243
scripts/plotting/mcsm_affinity_data_only.R Normal file
View file

@ -0,0 +1,243 @@
#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))
common_cols = c("mutationinformation"
, "X5uhc_position"
, "X5uhc_offset"
, "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_affinity
, scaled_cols_affinity
, outcome_cols_affinity
# , raw_cols_stability
# , scaled_cols_stability
# , outcome_cols_stability
)]
cols_to_extract = cols_to_consider[cols_to_consider%in%c(common_cols
, raw_cols_affinity
, scaled_cols_affinity)]
df3_plot = df3[, cols_to_extract]
DistCutOff_colnames = c("ligand_distance", "interface_dist")
DistCutOff = 10
df3_affinity_filtered = df3_plot[ (df3_plot$ligand_distance<10 | df3_plot$interface_dist <10),]
c0u = unique(df3_affinity_filtered$position)
length(c0u)
foo = df3_affinity_filtered[df3_affinity_filtered$ligand_distance<10,]
bar = df3_affinity_filtered[df3_affinity_filtered$interface_dist<10,]
wilcox.test(foo$mmcsm_lig_scaled~foo$sensitivity)
wilcox.test(foo$mmcsm_lig~foo$sensitivity)
wilcox.test(foo$affinity_scaled~foo$sensitivity)
wilcox.test(foo$ligand_affinity_change~foo$sensitivity)
wilcox.test(bar$mcsm_na_scaled~bar$sensitivity)
wilcox.test(bar$mcsm_na_affinity~bar$sensitivity)
wilcox.test(bar$mcsm_ppi2_scaled~bar$sensitivity)
wilcox.test(bar$mcsm_ppi2_affinity~bar$sensitivity)
##############################################################
df = df3_affinity_filtered
sum(is.na(df))
df2 = na.omit(df) # Apply na.omit function
a = df2[df2$position==37,]
sel_cols = c("mutationinformation", "position", scaled_cols_affinity)
a = a[, sel_cols]
##############################################################
# FIXME: ADD distance to NA when SP replies
#####################
# Ensemble affinity: affinity_cols
# mcsm_lig, mmcsm_lig and mcsm_na
#####################
# extract outcome cols and map numeric values to the categories
# Destabilising == 0, and stabilising == 1 so rescaling can let -1 be destabilising
#########################################
#=====================================
# Affintiy (2 cols): average the scores
# across predictors ==> average by
# position ==> scale b/w -1 and 1
# column to average: ens_affinity
#=====================================
cols_mcsm_lig = c("mutationinformation"
, "position"
, "sensitivity"
, "X5uhc_position"
, "X5uhc_offset"
, "ligand_distance"
, "ligand_outcome"
, "mmcsm_lig_outcome")
######################################################################
##################
# merge: mean ensemble stability and affinity by_position
####################
# if ( class(mean_ens_stability_by_position) && class(mean_ens_affinity_by_position) != "data.frame"){
# cat("Y")
# }
common_cols = intersect(colnames(mean_ens_stability_by_position), colnames(mean_ens_affinity_by_position))
if (dim(mean_ens_stability_by_position) && dim(mean_ens_affinity_by_position)){
print(paste0("PASS: dim's match, mering dfs by column :", common_cols))
#combined = as.data.frame(cbind(mean_duet_by_position, mean_affinity_by_position ))
combined_df = as.data.frame(merge(mean_ens_stability_by_position
, mean_ens_affinity_by_position
, by = common_cols
, all = T))
cat(paste0("\nnrows combined_df:", nrow(combined_df)
, "\nnrows combined_df:", ncol(combined_df)))
}else{
cat(paste0("FAIL: dim's mismatch, aborting cbind!"
, "\nnrows df1:", nrow(mean_duet_by_position)
, "\nnrows df2:", nrow(mean_affinity_by_position)))
quit()
}
#%%============================================================
# output
write.csv(combined_df, outfile_mean_ens_st_aff
, 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
#===============================================================

176
scripts/plotting/mcsm_mean_stability_ensemble_5uhc_rpob.R Normal file
View file

@ -0,0 +1,176 @@
#source("~/git/LSHTM_analysis/config/pnca.R")
#source("~/git/LSHTM_analysis/config/alr.R")
#source("~/git/LSHTM_analysis/config/gid.R")
#source("~/git/LSHTM_analysis/config/embb.R")
#source("~/git/LSHTM_analysis/config/katg.R")
#source("~/git/LSHTM_analysis/config/rpob.R")
source("/home/tanu/git/LSHTM_analysis/my_header.R")
#########################################################
# TASK: Generate averaged stability values
# across all stability tools
# for a given structure
#########################################################
#=======
# output
#=======
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
outfile_mean_ens_st_aff = paste0(outdir_images, "/5uhc_", tolower(gene)
, "_mean_ens_stability.csv")
print(paste0("Output file:", outfile_mean_ens_st_aff))
#%%===============================================================
#=============
# Input
#=============
df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv")
df3 = read.csv(df3_filename)
length(df3$mutationinformation)
# mut_info checks
table(df3$mutation_info)
table(df3$mutation_info_orig)
table(df3$mutation_info_labels_orig)
# used in plots and analyses
table(df3$mutation_info_labels) # different, and matches dst_mode
table(df3$dst_mode)
# create column based on dst mode with different colname
table(is.na(df3$dst))
table(is.na(df3$dst_mode))
#===============
# Create column: sensitivity mapped to dst_mode
#===============
df3$sensitivity = ifelse(df3$dst_mode == 1, "R", "S")
table(df3$sensitivity)
length(unique((df3$mutationinformation)))
all_colnames = as.data.frame(colnames(df3))
common_cols = c("mutationinformation"
, "X5uhc_position"
, "dst_mode"
, "mutation_info_labels"
, "sensitivity"
, "X5uhc_position"
, "X5uhc_offset"
, "ligand_distance"
, "interface_dist")
all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)]
all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)]
#===================
# stability cols
#===================
raw_cols_stability = c("duet_stability_change"
, "deepddg"
, "ddg_dynamut2"
, "ddg_foldx")
scaled_cols_stability = c("duet_scaled"
, "deepddg_scaled"
, "ddg_dynamut2_scaled"
, "foldx_scaled")
outcome_cols_stability = c("duet_outcome"
, "deepddg_outcome"
, "ddg_dynamut2_outcome"
, "foldx_outcome")
###########################################################
cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols
, raw_cols_stability
, scaled_cols_stability
, outcome_cols_stability)]
cols_to_extract = cols_to_consider[cols_to_consider%in%c(common_cols
, outcome_cols_stability)]
##############################################################
#####################
# Ensemble stability: outcome_cols_stability
#####################
# extract outcome cols and map numeric values to the categories
# Destabilising == 0, and stabilising == 1, so rescaling can let -1 be destabilising
df3_plot = df3[, cols_to_extract]
# assign numeric values to outcome
df3_plot[, outcome_cols_stability] <- sapply(df3_plot[, outcome_cols_stability]
, function(x){ifelse(x == "Destabilising", 0, 1)})
table(df3$duet_outcome)
table(df3_plot$duet_outcome)
#=====================================
# Stability (4 cols): average the scores
# across predictors ==> average by
# X5uhc_position ==> scale b/w -1 and 1
# column to average: ens_stability
#=====================================
cols_to_average = which(colnames(df3_plot)%in%outcome_cols_stability)
# ensemble average across predictors
df3_plot$ens_stability = rowMeans(df3_plot[,cols_to_average])
head(df3_plot$X5uhc_position); head(df3_plot$mutationinformation)
head(df3_plot$ens_stability)
table(df3_plot$ens_stability)
# ensemble average of predictors by X5uhc_position
mean_ens_stability_by_position <- df3_plot %>%
dplyr::group_by(X5uhc_position) %>%
dplyr::summarize(avg_ens_stability = mean(ens_stability))
# REscale b/w -1 and 1
#en_stab_min = min(mean_ens_stability_by_position['avg_ens_stability'])
#en_stab_max = max(mean_ens_stability_by_position['avg_ens_stability'])
# scale the average stability value between -1 and 1
# mean_ens_by_position['averaged_stability3_scaled'] = lapply(mean_ens_by_position['averaged_stability3']
# , function(x) ifelse(x < 0, x/abs(en3_min), x/en3_max))
mean_ens_stability_by_position['avg_ens_stability_scaled'] = lapply(mean_ens_stability_by_position['avg_ens_stability']
, function(x) {
scales::rescale(x, to = c(-1,1)
#, from = c(en_stab_min,en_stab_max))
, from = c(0,1))
})
cat(paste0('Average stability scores:\n'
, head(mean_ens_stability_by_position['avg_ens_stability'])
, '\n---------------------------------------------------------------'
, '\nAverage stability scaled scores:\n'
, head(mean_ens_stability_by_position['avg_ens_stability_scaled'])))
# convert to a data frame
mean_ens_stability_by_position = as.data.frame(mean_ens_stability_by_position)
#FIXME: sanity checks
# TODO: predetermine the bounds
# l_bound_ens = min(mean_ens_stability_by_position['avg_ens_stability_scaled'])
# u_bound_ens = max(mean_ens_stability_by_position['avg_ens_stability_scaled'])
#
# if ( (l_bound_ens == -1) && (u_bound_ens == 1) ){
# cat(paste0("PASS: ensemble stability scores averaged by X5uhc_position and then scaled"
# , "\nmin ensemble averaged stability: ", l_bound_ens
# , "\nmax ensemble averaged stability: ", u_bound_ens))
# }else{
# cat(paste0("FAIL: avergaed duet scores could not be scaled b/w -1 and 1"
# , "\nmin ensemble averaged stability: ", l_bound_ens
# , "\nmax ensemble averaged stability: ", u_bound_ens))
# quit()
# }
##################################################################
# output
#write.csv(combined_df, outfile_mean_ens_st_aff
write.csv(mean_ens_stability_by_position
, outfile_mean_ens_st_aff
, row.names = F)
cat("Finished writing file:\n"
, outfile_mean_ens_st_aff
, "\nNo. of rows:", nrow(mean_ens_stability_by_position)
, "\nNo. of cols:", ncol(mean_ens_stability_by_position))
# end of script
#===============================================================

283
scripts/plotting/replaceBfactor_pdb_stability_5uhc_rpob.R Normal file
View file

@ -0,0 +1,283 @@
#!/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/5uhc_rpob_", 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(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???
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!