From a448d9276b2131e9e7031a33328ef10d9b929ab4 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Fri, 21 Aug 2020 13:25:01 +0100 Subject: [PATCH] added plotting scripts from old run --- scripts/plotting/combining_two_df_FIXME.R | 442 ++++++++++++++++++++++ scripts/plotting/replaceBfactor_pdb.R | 318 ++++++++++++++++ 2 files changed, 760 insertions(+) create mode 100644 scripts/plotting/combining_two_df_FIXME.R create mode 100644 scripts/plotting/replaceBfactor_pdb.R diff --git a/scripts/plotting/combining_two_df_FIXME.R b/scripts/plotting/combining_two_df_FIXME.R new file mode 100644 index 0000000..0d34e06 --- /dev/null +++ b/scripts/plotting/combining_two_df_FIXME.R @@ -0,0 +1,442 @@ +getwd() +setwd("~/git/LSHTM_analysis/scripts/plotting/") +getwd() + +######################################################### +# TASK: To combine struct params and meta data for plotting +# Input csv files: +# 1) _all_params.csv +# 2) _meta_data.csv + +# Output: +# 1) muts with opposite effects on stability +# 2) large combined df including NAs for AF, OR,etc +# Dim: same no. of rows as gene associated meta_data_with_AFandOR +# 3) small combined df including NAs for AF, OR, etc. +# Dim: same as mcsm data +# 4) large combined df excluding NAs +# Dim: dim(#1) - no. of NAs(AF|OR) + 1 +# 5) small combined df excluding NAs +# Dim: dim(#2) - no. of unique NAs - 1 +# This script is sourced from other .R scripts for plotting +######################################################### + +########################################################## +# Installing and loading required packages +########################################################## +source("Header_TT.R") +#require(data.table) +#require(arsenal) +#require(compare) +#library(tidyverse) + + +#%% variable assignment: input and output paths & filenames +drug = "pyrazinamide" +gene = "pncA" +gene_match = paste0(gene,"_p.") +cat(gene_match) + +#============= +# directories +#============= +datadir = paste0("~/git/Data") +indir = paste0(datadir, "/", drug, "/input") +outdir = paste0("~/git/Data", "/", drug, "/output") + +#=========== +# input +#=========== +#in_filename = "mcsm_complex1_normalised.csv" +in_filename_params = paste0(tolower(gene), "_all_params.csv") +infile_params = paste0(outdir, "/", in_filename_params) +cat(paste0("Input file 1:", infile_params) ) + +# infile 2: gene associated meta data +#in_filename_gene_metadata = paste0(tolower(gene), "_meta_data_with_AFandOR.csv") +in_filename_gene_metadata = paste0(tolower(gene), "_metadata.csv") +infile_gene_metadata = paste0(outdir, "/", in_filename_gene_metadata) +cat(paste0("Input infile 2:", infile_gene_metadata)) + +#=========== +# output +#=========== +# mutations with opposite effects +out_filename_opp_muts = paste0(tolower(gene), "_muts_opp_effects.csv") +outfile_opp_muts = paste0(outdir, "/", out_filename_opp_muts) + + +#%%=============================================================== +########################### +# Read file: struct params +########################### +cat("Reading struct params including mcsm:" + , in_filename_params) + +mcsm_data = read.csv(infile_params + #, row.names = 1 + , stringsAsFactors = F + , header = T) + +cat("Input dimensions:", dim(mcsm_data)) #416, 86 + +# clear variables +rm(in_filename_params, infile_params) + +str(mcsm_data) + +table(mcsm_data$duet_outcome); sum(table(mcsm_data$duet_outcome) ) + +# spelling Correction 1: DUET incase American spelling needed! +#mcsm_data$duet_outcome[mcsm_data$duet_outcome=="Stabilising"] <- "Stabilizing" +#mcsm_data$duet_outcome[mcsm_data$duet_outcome=="Destabilising"] <- "Destabilizing" + +# checks: should be the same as above +table(mcsm_data$duet_outcome); sum(table(mcsm_data$duet_outcome) ) +head(mcsm_data$duet_outcome); tail(mcsm_data$duet_outcome) + +# spelling Correction 2: Ligand incase American spelling needed! +table(mcsm_data$ligand_outcome); sum(table(mcsm_data$ligand_outcome) ) +#mcsm_data$ligand_outcome[mcsm_data$ligand_outcome=="Stabilising"] <- "Stabilizing" +#mcsm_data$ligand_outcome[mcsm_data$ligand_outcome=="Destabilising"] <- "Destabilizing" + +# checks: should be the same as above +table(mcsm_data$ligand_outcome); sum(table(mcsm_data$ligand_outcome) ) +head(mcsm_data$ligand_outcome); tail(mcsm_data$ligand_outcome) + +# muts with opposing effects on protomer and ligand stability +table(mcsm_data$duet_outcome != mcsm_data$ligand_outcome) +changes = mcsm_data[which(mcsm_data$duet_outcome != mcsm_data$ligand_outcome),] + +# sanity check: redundant, but uber cautious! +dl_i = which(mcsm_data$duet_outcome != mcsm_data$ligand_outcome) +ld_i = which(mcsm_data$ligand_outcome != mcsm_data$duet_outcome) + +cat("Identifying muts with opposite stability effects") +if(nrow(changes) == (table(mcsm_data$duet_outcome != mcsm_data$ligand_outcome)[[2]]) & identical(dl_i,ld_i)) { + cat("PASS: muts with opposite effects on stability and affinity correctly identified" + , "\nNo. of such muts: ", nrow(changes)) +}else { + cat("FAIL: unsuccessful in extracting muts with changed stability effects") +} + +#*************************** +# write file: changed muts +write.csv(changes, outfile_opp_muts) + +cat("Finished writing file for muts with opp effects:" + , "\nFilename: ", outfile_opp_muts + , "\nDim:", dim(changes)) + +# clear variables +rm(out_filename_opp_muts, outfile_opp_muts) +rm(changes, dl_i, ld_i) + +#*************************** +# count na in each column +na_count = sapply(mcsm_data, function(y) sum(length(which(is.na(y))))); na_count + +# sort by mutationinformation +##mcsm_data = mcsm_data[order(mcsm_data$mutationinformation),] +##head(mcsm_data$mutationinformation) + +df_ncols = ncol(mcsm_data) + +# REMOVE as this is dangerous due to dup muts +# get freq count of positions and add to the df +#setDT(mcsm_data)[, occurrence := .N, by = .(position)] + +#cat("Added 1 col: position frequency to see which posn has how many muts" +# , "\nNo. of cols now", ncol(mcsm_data) +# , "\nNo. of cols before: ", df_ncols) + +#pos_count_check = data.frame(mcsm_data$position, mcsm_data$occurrence) + +# check duplicate muts +if (length(unique(mcsm_data$mutationinformation)) == length(mcsm_data$mutationinformation)){ + cat("No duplicate mutations in mcsm data") +}else{ + dup_muts = mcsm_data[duplicated(mcsm_data$mutationinformation),] + dup_muts_nu = length(unique(dup_muts$mutationinformation)) + cat(paste0("CAUTION:", nrow(dup_muts), " Duplicate mutations identified" + , "\nOf these, no. of unique mutations are:", dup_muts_nu + , "\nExtracting df with unique mutations only")) + mcsm_data_u = mcsm_data[!duplicated(mcsm_data$mutationinformation),] +} + +if (nrow(mcsm_data_u) == length(unique(mcsm_data$mutationinformation))){ + cat("Df without duplicate mutations successfully extracted") +} else{ + cat("FAIL: could not extract clean df!") + quit() +} + +########################### +# 2: Read file: _meta data.csv +########################### +cat("Reading meta data file:", infile_gene_metadata) + +gene_metadata <- read.csv(infile_gene_metadata + , stringsAsFactors = F + , header = T) +cat("Dim:", dim(gene_metadata)) + +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +# FIXME: remove +# counting NAs in AF, OR cols: +if (identical(sum(is.na(gene_metadata$OR)) + , sum(is.na(gene_metadata$pvalue)) + , sum(is.na(gene_metadata$AF)))){ + cat("PASS: NA count match for OR, pvalue and AF\n") + na_count = sum(is.na(gene_metadata$AF)) + cat("No. of NAs: ", sum(is.na(gene_metadata$OR))) +} else{ + cat("FAIL: NA count mismatch" + , "\nNA in OR: ", sum(is.na(gene_metadata$OR)) + , "\nNA in pvalue: ", sum(is.na(gene_metadata$pvalue)) + , "\nNA in AF:", sum(is.na(gene_metadata$AF))) +} +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +# clear variables +rm(in_filename_gene_metadata, infile_gene_metadata) + +str(gene_metadata) + +# sort by position (same as mcsm_data) +# earlier it was mutationinformation +#head(gene_metadata$mutationinformation) +#gene_metadata = gene_metadata[order(gene_metadata$mutationinformation),] +##head(gene_metadata$mutationinformation) + +head(gene_metadata$position) +gene_metadata = gene_metadata[order(gene_metadata$position),] +head(gene_metadata$position) + +########################### +# Merge 1: two dfs with NA +# merged_df2 +########################### +head(mcsm_data$mutationinformation) +head(gene_metadata$mutationinformation) + +# Find common columns b/w two df +merging_cols = intersect(colnames(mcsm_data), colnames(gene_metadata)) + +cat(paste0("Merging dfs with NAs: big df (1-many relationship b/w id & mut)" + , "\nNo. of merging cols:", length(merging_cols) + , "\nMerging columns identified:")) +print(merging_cols) + +#============= +# merged_df2): gene_metadata + mcsm_data +#============== +merged_df2 = merge(x = gene_metadata + , y = mcsm_data + , by = merging_cols + , all.y = T) + +cat("Dim of merged_df2: ", dim(merged_df2) #4520, 11 + ) +head(merged_df2$position) + +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +# FIXME: count how many unique muts have entries in meta data +# sanity check +cat("Checking nrows in merged_df2") +if(nrow(gene_metadata) == nrow(merged_df2)){ + cat("nrow(merged_df2) = nrow (gene associated gene_metadata)" + ,"\nExpected no. of rows: ",nrow(gene_metadata) + ,"\nGot no. of rows: ", nrow(merged_df2)) +} else{ + cat("nrow(merged_df2)!= nrow(gene associated gene_metadata)" + , "\nExpected no. of rows after merge: ", nrow(gene_metadata) + , "\nGot no. of rows: ", nrow(merged_df2) + , "\nFinding discrepancy") + merged_muts_u = unique(merged_df2$mutationinformation) + meta_muts_u = unique(gene_metadata$mutationinformation) + # find the index where it differs + unique(meta_muts_u[! meta_muts_u %in% merged_muts_u]) +} + +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +# sort by position +head(merged_df2$position) +merged_df2 = merged_df2[order(merged_df2$position),] +head(merged_df2$position) + +merged_df2v3 = merge(x = gene_metadata + , y = mcsm_data + , by = merging_cols + , all = T) + +merged_df2v2 = merge(x = gene_metadata + , y = mcsm_data + , by = merging_cols + , all.x = T) +#!=!=!=!=!=!=!=! +# COMMENT: used all.y since position 186 is not part of the struc, +# hence doesn"t have a mcsm value +# but 186 is associated with mutation +#!=!=!=!=!=!=!=! + +# should be False +identical(merged_df2, merged_df2v2) +table(merged_df2$position%in%merged_df2v2$position) + +rm(merged_df2v2) + +#!!!!!!!!!!! check why these differ + +######### +# merge 3b (merged_df3):remove duplicate mutation information +######### +cat("Merging dfs without NAs: small df (removing muts with no AF|OR associated)" + ,"\nCannot trust lineage info from this" + ,"\nlinking col: Mutationinforamtion" + ,"\nfilename: merged_df3") + +#==#=#=#=#=#=# +# Cannot trust lineage, country from this df as the same mutation +# can have many different lineages +# but this should be good for the numerical corr plots +#=#=#=#=#=#=#= +merged_df3 = merged_df2[!duplicated(merged_df2$mutationinformation),] +head(merged_df3$position); tail(merged_df3$position) # should be sorted + +# sanity check +cat("Checking nrows in merged_df3") +if(nrow(mcsm_data) == nrow(merged_df3)){ + cat("PASS: No. of rows match with mcsm_data" + ,"\nExpected no. of rows: ", nrow(mcsm_data) + ,"\nGot no. of rows: ", nrow(merged_df3)) +} else { + cat("FAIL: No. of rows mismatch" + , "\nNo. of rows mcsm_data: ", nrow(mcsm_data) + , "\nNo. of rows merged_df3: ", nrow(merged_df3)) +} + +# counting NAs in AF, OR cols in merged_df3 +# this is becuase mcsm has no AF, OR cols, +# so you cannot count NAs +if (identical(sum(is.na(merged_df3$OR)) + , sum(is.na(merged_df3$pvalue)) + , sum(is.na(merged_df3$AF)))){ + cat("PASS: NA count match for OR, pvalue and AF\n") + na_count_df3 = sum(is.na(merged_df3$AF)) + cat("No. of NAs: ", sum(is.na(merged_df3$OR))) +} else{ + cat("FAIL: NA count mismatch" + , "\nNA in OR: ", sum(is.na(merged_df3$OR)) + , "\nNA in pvalue: ", sum(is.na(merged_df3$pvalue)) + , "\nNA in AF:", sum(is.na(merged_df3$AF))) +} + +########################### +# 4: merging two dfs: without NA +########################### +######### +# merge 4a (merged_df2_comp): same as merge 1 but excluding NA +######### +cat("Merging dfs without any NAs: big df (1-many relationship b/w id & mut)" + ,"\nlinking col: Mutationinforamtion" + ,"\nfilename: merged_df2_comp") + +merged_df2_comp = merged_df2[!is.na(merged_df2$AF),] +#merged_df2_comp = merged_df2[!duplicated(merged_df2$mutationinformation),] + +# sanity check +cat("Checking nrows in merged_df2_comp") +if(nrow(merged_df2_comp) == (nrow(merged_df2) - na_count + 1)){ + cat("PASS: No. of rows match" + ,"\nDim of merged_df2_comp: " + ,"\nExpected no. of rows: ", nrow(merged_df2) - na_count + 1 + , "\nNo. of rows: ", nrow(merged_df2_comp) + , "\nNo. of cols: ", ncol(merged_df2_comp)) +}else{ + cat("FAIL: No. of rows mismatch" + ,"\nExpected no. of rows: ", nrow(merged_df2) - na_count + 1 + ,"\nGot no. of rows: ", nrow(merged_df2_comp)) +} + +######### +# merge 4b (merged_df3_comp): remove duplicate mutation information +######### +merged_df3_comp = merged_df2_comp[!duplicated(merged_df2_comp$mutationinformation),] + +cat("Dim of merged_df3_comp: " + , "\nNo. of rows: ", nrow(merged_df3_comp) + , "\nNo. of cols: ", ncol(merged_df3_comp)) + +# alternate way of deriving merged_df3_comp +foo = merged_df3[!is.na(merged_df3$AF),] +# compare dfs: foo and merged_df3_com +all.equal(foo, merged_df3) + +summary(comparedf(foo, merged_df3)) + +# sanity check +cat("Checking nrows in merged_df3_comp") +if(nrow(merged_df3_comp) == nrow(merged_df3)){ + cat("NO NAs detected in merged_df3 in AF|OR cols" + ,"\nNo. of rows are identical: ", nrow(merged_df3)) +} else{ + if(nrow(merged_df3_comp) == nrow(merged_df3) - na_count_df3) { + cat("PASS: NAs detected in merged_df3 in AF|OR cols" + , "\nNo. of NAs: ", na_count_df3 + , "\nExpected no. of rows in merged_df3_comp: ", nrow(merged_df3) - na_count_df3 + , "\nGot no. of rows: ", nrow(merged_df3_comp)) + } +} + +#=============== end of combining df +#********************* +# writing 1 file in the style of a loop: merged_df3 +# print(output dir) +#i = "merged_df3" +#out_filename = paste0(i, ".csv") +#outfile = paste0(outdir, "/", out_filename) + +#cat("Writing output file: " +# ,"\nFilename: ", out_filename +# ,"\nPath: ", outdir) + +#template: write.csv(merged_df3, "merged_df3.csv") +#write.csv(get(i), outfile, row.names = FALSE) +#cat("Finished writing: ", outfile +# , "\nNo. of rows: ", nrow(get(i)) +# , "\nNo. of cols: ", ncol(get(i))) + +#%% write_output files; all 4 files: +outvars = c("merged_df2" + , "merged_df3" + , "merged_df2_comp" + , "merged_df3_comp") + +cat("Writing output files: " + , "\nPath:", outdir) + +for (i in outvars){ +# cat(i, "\n") + out_filename = paste0(i, ".csv") +# cat(out_filename, "\n") +# cat("getting value of variable: ", get(i)) + outfile = paste0(outdir, "/", out_filename) +# cat("Full output path: ", outfile, "\n") + cat("Writing output file:" + ,"\nFilename: ", out_filename,"\n") + write.csv(get(i), outfile, row.names = FALSE) + cat("Finished writing: ", outfile + , "\nNo. of rows: ", nrow(get(i)) + , "\nNo. of cols: ", ncol(get(i)), "\n") +} + +# alternate way to replace with implicit loop +# FIXME +#sapply(outvars, function(x, y) write.csv(get(outvars), paste0(outdir, "/", outvars, ".csv"))) +#************************* +# clear variables +rm(mcsm_data, gene_metadata, foo, drug, gene, gene_match, indir, merged_muts_u, meta_muts_u, na_count, df_ncols, outdir) +rm(pos_count_check) +#============================= end of script + diff --git a/scripts/plotting/replaceBfactor_pdb.R b/scripts/plotting/replaceBfactor_pdb.R new file mode 100644 index 0000000..0104059 --- /dev/null +++ b/scripts/plotting/replaceBfactor_pdb.R @@ -0,0 +1,318 @@ +#!/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 respecitve 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("Header_TT.R") +library(bio3d) +require("getopt", quietly = TRUE) # cmd parse arguments +#======================================================== +# 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)") +} +#======================================================== +#%% variable assignment: input and output paths & filenames +drug = "pyrazinamide" +gene = "pncA" +gene_match = paste0(gene,"_p.") +cat(gene_match) + +#============= +# directories +#============= +datadir = paste0("~/git/Data") +indir = paste0(datadir, "/", drug, "/input") +outdir = paste0("~/git/Data", "/", drug, "/output") + +#====== +# 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) +cat(paste0("Input file:", infile_mean_stability) ) + +#======= +# output +#======= +out_filename_duet_mspdb = paste0(tolower(gene), "_complex_b_duetms.pdb") +outfile_duet_mspdb = paste0(outdir, "/", out_filename_duet_mspdb) +print(paste0("Output file:", outfile_duet_mspdb)) + +out_filename_lig_mspdb = paste0(tolower(gene), "_complex_b_ligms.pdb") +outfile_lig_mspdb = paste0(outdir, "/", out_filename_lig_mspdb) +print(paste0("Output file:", outfile_lig_mspdb)) + +#%%=============================================================== +########################### +# 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_lig = my_pdb_lig[[1]] + +# 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 duet") + +plot(density(df_duet$b) + , xlab = "" + , main = "Bfactor duet") + + +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$average_duet_scaled + , xlab = "" + , main = "mean duet scaled") + +plot(density(my_df$average_duet_scaled) + , xlab = "" + , main = "mean duet scaled") + + +hist(my_df$average_affinity_scaled + , xlab = "" + , main = "mean affinity scaled") + +plot(density(my_df$average_affinity_scaled) + , xlab = "" + , main = "mean affinity scaled") + +#************ +# 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$average_duet_scaled[match(df_duet$resno, my_df$position)] +df_lig$b = my_df$average_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 +df_duet$b[is.na(df_duet$b)] = 0 +df_lig$b[is.na(df_lig$b)] = 0 + +# sanity check: should be 0 and True +# duet +if (sum(df_duet$b == 0) == b_na_duet){ + print ("PASS: NA"s replaced with 0"s successfully in df_duet") +} else { + print("FAIL: NA replacement in df_duet NOT successful") + quit() +} +max(df_duet$b); min(df_duet$b) + +# lig +if (sum(df_lig$b == 0) == b_na_lig){ + print ("PASS: NA"s replaced with 0"s successfully df_lig") +} else { + print("FAIL: NA replacement in df_lig NOT successful") + quit() +} +max(df_lig$b); min(df_lig$b) + +# sanity checks: should be True +if( (max(df_duet$b) == max(my_df$average_duet_scaled)) & (min(df_duet$b) == min(my_df$average_duet_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$average_affinity_scaled)) & (min(df_lig$b) == min(my_df$average_affinity_scaled)) ){ + print("PASS: B-factors replaced correctly in lig_duet") +} else { + print ("FAIL: To replace B-factors in lig_duet") + 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[[1]] = df_duet +max(df_duet$b); min(df_duet$b) + +my_pdb_lig[[1]] = 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 = "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. This is because all the positions +# where there are no SNPs have been assigned 0??? +#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +