added scripts
This commit is contained in:
parent
dd69da01f6
commit
23b4f06017
10 changed files with 147 additions and 1014 deletions
|
|
@ -1,46 +1,23 @@
|
|||
#!/usr/bin/env Rscript
|
||||
source("~/git/LSHTM_analysis/config/gid.R")
|
||||
|
||||
source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
|
||||
|
||||
#########################################################
|
||||
# TASK: Replace B-factors in the pdb file with the mean
|
||||
# normalised stability values.
|
||||
|
||||
# read pdb file
|
||||
|
||||
# read mcsm mean stability value files
|
||||
# extract the respective mean values and assign to the
|
||||
# b-factor column within their respective pdbs
|
||||
|
||||
# generate some distribution plots for inspection
|
||||
# normalised affinity values
|
||||
|
||||
#########################################################
|
||||
# working dir and loading libraries
|
||||
getwd()
|
||||
setwd("~/git/LSHTM_analysis/scripts/plotting")
|
||||
cat(c(getwd(),"\n"))
|
||||
|
||||
#source("~/git/LSHTM_analysis/scripts/Header_TT.R")
|
||||
library(bio3d)
|
||||
require("getopt", quietly = TRUE) # cmd parse arguments
|
||||
#========================================================
|
||||
#drug = "pyrazinamide"
|
||||
#gene = "pncA"
|
||||
|
||||
# command line args
|
||||
spec = matrix(c(
|
||||
"drug" , "d", 1, "character",
|
||||
"gene" , "g", 1, "character"
|
||||
), byrow = TRUE, ncol = 4)
|
||||
|
||||
opt = getopt(spec)
|
||||
|
||||
drug = opt$drug
|
||||
gene = opt$gene
|
||||
|
||||
if(is.null(drug)|is.null(gene)) {
|
||||
stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)")
|
||||
}
|
||||
#========================================================
|
||||
gene_match = paste0(gene,"_p.")
|
||||
cat(gene)
|
||||
gene_match = paste0(gene,"_p."); cat(gene_match)
|
||||
cat(gene_match)
|
||||
|
||||
#=============
|
||||
|
|
@ -49,9 +26,13 @@ cat(gene_match)
|
|||
datadir = paste0("~/git/Data")
|
||||
indir = paste0(datadir, "/", drug, "/input")
|
||||
outdir = paste0("~/git/Data", "/", drug, "/output")
|
||||
#outdir_plots = paste0("~/git/Data", "/", drug, "/output/plots")
|
||||
outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
|
||||
#outdir_plots = paste0("~/git/Writing/thesis/images/results/", tolower(gene))
|
||||
|
||||
#=======
|
||||
# output
|
||||
#=======
|
||||
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene), "/")
|
||||
cat("plots will output to:", outdir_images)
|
||||
#======
|
||||
# input
|
||||
#======
|
||||
|
|
@ -59,31 +40,31 @@ in_filename_pdb = paste0(tolower(gene), "_complex.pdb")
|
|||
infile_pdb = paste0(indir, "/", in_filename_pdb)
|
||||
cat(paste0("Input file:", infile_pdb) )
|
||||
|
||||
#in_filename_mean_stability = paste0(tolower(gene), "_mean_stability.csv")
|
||||
#infile_mean_stability = paste0(outdir, "/", in_filename_mean_stability)
|
||||
|
||||
in_filename_mean_affinity = paste0(tolower(gene), "_mean_ligand.csv")
|
||||
infile_mean_affinity = paste0(outdir_plots, "/", in_filename_mean_affinity)
|
||||
|
||||
cat(paste0("Input file:", infile_mean_affinity) )
|
||||
|
||||
#=======
|
||||
# output
|
||||
#=======
|
||||
#out_filename_duet_mspdb = paste0(tolower(gene), "_complex_bduet_ms.pdb")
|
||||
out_filename_lig_mspdb = paste0(tolower(gene), "_complex_b_lig_ms.pdb")
|
||||
outfile_lig_mspdb = paste0(outdir_plots, "/", out_filename_lig_mspdb)
|
||||
outfile_lig_mspdb = paste0(outdir_images,out_filename_lig_mspdb)
|
||||
print(paste0("Output file:", outfile_lig_mspdb))
|
||||
|
||||
#%%===============================================================
|
||||
#NOTE: duet here refers to the ensemble stability values
|
||||
#NOTE: duet here refers to the ensemble affinity values
|
||||
|
||||
###########################
|
||||
# Read file: average stability values
|
||||
# Read file: average affinity values
|
||||
# or mcsm_normalised file
|
||||
###########################
|
||||
my_df <- read.csv(infile_mean_stability, header = T)
|
||||
str(my_df)
|
||||
my_df_raw = merged_df3[, c("position", "ligand_distance", "avg_lig_affinity_scaled", "avg_lig_affinity")]
|
||||
my_df_raw = my_df_raw[my_df_raw$ligand_distance<10,]
|
||||
|
||||
# avg by position on the SCALED values
|
||||
my_df <- my_df_raw %>%
|
||||
group_by(position) %>%
|
||||
summarize(avg_ligaff_sc_pos = mean(avg_lig_affinity_scaled))
|
||||
|
||||
max(my_df$avg_ligaff_sc_pos)
|
||||
min(my_df$avg_ligaff_sc_pos)
|
||||
|
||||
#############
|
||||
# Read pdb
|
||||
|
|
@ -98,13 +79,11 @@ my_pdb = read.pdb(infile_pdb
|
|||
, hex = FALSE
|
||||
, verbose = TRUE)
|
||||
|
||||
rm(in_filename_mean_affinity, in_filename_pdb)
|
||||
|
||||
# assign separately for duet and ligand
|
||||
my_pdb_duet = my_pdb
|
||||
|
||||
#=========================================================
|
||||
# Replacing B factor with mean stability scores
|
||||
# Replacing B factor with mean affinity scores
|
||||
# within the respective dfs
|
||||
#==========================================================
|
||||
# extract atom list into a variable
|
||||
|
|
@ -121,8 +100,8 @@ max(df_duet$b); min(df_duet$b)
|
|||
#==================================================
|
||||
# histograms and density plots for inspection
|
||||
# 1: original B-factors
|
||||
# 2: original mean stability values
|
||||
# 3: replaced B-factors with mean stability values
|
||||
# 2: original mean affinity values
|
||||
# 3: replaced B-factors with mean affinity values
|
||||
#==================================================
|
||||
# Set the margin on all sides
|
||||
par(oma = c(3,2,3,0)
|
||||
|
|
@ -131,6 +110,7 @@ par(oma = c(3,2,3,0)
|
|||
#, mfrow = c(3,4))
|
||||
, mfrow = c(3,2))
|
||||
|
||||
|
||||
#=============
|
||||
# Row 1 plots: original B-factors
|
||||
# duet and affinity
|
||||
|
|
@ -144,40 +124,28 @@ plot(density(df_duet$b)
|
|||
, main = "Bfactor affinity")
|
||||
|
||||
#=============
|
||||
# Row 2 plots: original mean stability values
|
||||
# duet and affinity
|
||||
# Row 2 plots: original mean affinity values
|
||||
# affinity
|
||||
#=============
|
||||
|
||||
#hist(my_df$averaged_duet
|
||||
hist(my_df$avg_lig_scaled
|
||||
hist(my_df$avg_ligaff_sc_pos
|
||||
, xlab = ""
|
||||
, main = "mean affinity values")
|
||||
|
||||
#plot(density(my_df$averaged_duet)
|
||||
plot(density(my_df$avg_lig_scaled)
|
||||
plot(density(my_df$avg_ligaff_sc_pos)
|
||||
, xlab = ""
|
||||
, main = "mean affinity values")
|
||||
|
||||
#==============
|
||||
# Row 3 plots: replaced B-factors with mean stability values
|
||||
# After actual replacement in the b factor column
|
||||
#===============
|
||||
################################################################
|
||||
#=========
|
||||
# step 0_P1: DONT RUN once you have double checked the matched output
|
||||
#=========
|
||||
# sanity check: match and assign to a separate column to double check
|
||||
# colnames(my_df)
|
||||
# df_duet$duet_scaled = my_df$averge_duet_scaled[match(df_duet$resno, my_df$position)]
|
||||
|
||||
|
||||
#=========
|
||||
# step 1_P1
|
||||
#=========
|
||||
# Be brave and replace in place now (don"t run sanity check)
|
||||
# this makes all the B-factor values in the non-matched positions as NA
|
||||
|
||||
#df_duet$b = my_df$averaged_duet_scaled[match(df_duet$resno, my_df$position)]
|
||||
df_duet$b = my_df$avg_lig_scaled[match(df_duet$resno, my_df$position)]
|
||||
df_duet$b = my_df$avg_ligaff_sc_pos[match(df_duet$resno, my_df$position)]
|
||||
|
||||
#=========
|
||||
# step 2_P1
|
||||
|
|
@ -192,32 +160,6 @@ sum(df_duet$b == 0)
|
|||
na_rep = 2
|
||||
df_duet$b[is.na(df_duet$b)] = na_rep
|
||||
|
||||
# # sanity check: should be 0 and True
|
||||
# # duet and lig
|
||||
# if ( (sum(df_duet$b == na_rep) == b_na_duet) && (sum(df_lig$b == na_rep) == b_na_lig) ) {
|
||||
# print ("PASS: NA's replaced with 0s successfully in df_duet and df_lig")
|
||||
# } else {
|
||||
# print("FAIL: NA replacement in df_duet NOT successful")
|
||||
# quit()
|
||||
# }
|
||||
#
|
||||
# max(df_duet$b); min(df_duet$b)
|
||||
#
|
||||
# # sanity checks: should be True
|
||||
# if( (max(df_duet$b) == max(my_df$avg_ens_stability_scaled)) & (min(df_duet$b) == min(my_df$avg_ens_stability_scaled)) ){
|
||||
# print("PASS: B-factors replaced correctly in df_duet")
|
||||
# } else {
|
||||
# print ("FAIL: To replace B-factors in df_duet")
|
||||
# quit()
|
||||
# }
|
||||
|
||||
# if( (max(df_lig$b) == max(my_df$avg_ens_affinity_scaled)) & (min(df_lig$b) == min(my_df$avg_ens_affinity_scaled)) ){
|
||||
# print("PASS: B-factors replaced correctly in df_lig")
|
||||
# } else {
|
||||
# print ("FAIL: To replace B-factors in df_lig")
|
||||
# quit()
|
||||
# }
|
||||
|
||||
#=========
|
||||
# step 3_P1
|
||||
#=========
|
||||
|
|
@ -241,17 +183,23 @@ table(df_duet$b)
|
|||
sum(is.na(df_duet$b))
|
||||
|
||||
#=========
|
||||
# step 5_P1
|
||||
# step 5_P1: OUTPUT
|
||||
#=========
|
||||
cat(paste0("output file duet mean stability pdb:"
|
||||
, outfile_lig_mspdb))
|
||||
cat(paste0("output file duet mean affinity pdb:", outfile_lig_mspdb))
|
||||
write.pdb(my_pdb_duet, outfile_lig_mspdb)
|
||||
|
||||
# OUTPUT: position file
|
||||
poscsvF = paste0(outdir_images, tolower(gene), "_ligaff_positions.csv")
|
||||
cat(paste0("output file duet mean NA affinity POSITIONS:", poscsvF))
|
||||
|
||||
filtered_pos = toString(my_df$position)
|
||||
write.table(filtered_pos, poscsvF, row.names = F, col.names = F )
|
||||
|
||||
#============================
|
||||
# Add the 3rd histogram and density plots for comparisons
|
||||
#============================
|
||||
# Plots continued...
|
||||
# Row 3 plots: hist and density of replaced B-factors with stability values
|
||||
# Row 3 plots: hist and density of replaced B-factors with affinity values
|
||||
hist(df_duet$b
|
||||
, xlab = ""
|
||||
, main = "repalcedB duet")
|
||||
|
|
@ -266,16 +214,8 @@ mtext(text = "Frequency"
|
|||
, line = 0
|
||||
, outer = TRUE)
|
||||
|
||||
mtext(text = paste0(tolower(gene), ": afinity distribution")
|
||||
mtext(text = paste0(tolower(gene), ": affinity distribution")
|
||||
, side = 3
|
||||
, line = 0
|
||||
, outer = TRUE)
|
||||
#============================================
|
||||
|
||||
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
# NOTE: This replaced B-factor distribution has the same
|
||||
# x-axis as the PredAff normalised values, but the distribution
|
||||
# is affected since 0 is overinflated/or hs an additional blip because
|
||||
# of the positions not associated with resistance. This is because all the positions
|
||||
# where there are no SNPs have been assigned 0???
|
||||
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||
Loading…
Add table
Add a link
Reference in a new issue