213 lines
No EOL
5.6 KiB
R
213 lines
No EOL
5.6 KiB
R
#!/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 ppi2 values.
|
|
#########################################################
|
|
|
|
#source("~/git/LSHTM_analysis/scripts/Header_TT.R")
|
|
library(bio3d)
|
|
require("getopt", quietly = TRUE) # cmd parse arguments
|
|
|
|
#========================================================
|
|
cat(gene)
|
|
gene_match = paste0(gene,"_p."); cat(gene_match)
|
|
cat(gene_match)
|
|
|
|
#=============
|
|
# directories
|
|
#=============
|
|
datadir = paste0("~/git/Data")
|
|
indir = paste0(datadir, "/", drug, "/input")
|
|
outdir = paste0("~/git/Data", "/", drug, "/output")
|
|
#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
|
|
#======
|
|
in_filename_pdb = paste0(tolower(gene), "_complex.pdb")
|
|
infile_pdb = paste0(indir, "/", in_filename_pdb)
|
|
cat(paste0("Input file:", infile_pdb) )
|
|
|
|
#=======
|
|
# output
|
|
#=======
|
|
out_filename_ppi2_mspdb = paste0(tolower(gene), "_complex_b_ppi2_ms.pdb")
|
|
outfile_ppi2_mspdb = paste0(outdir_images,out_filename_ppi2_mspdb)
|
|
print(paste0("Output file:", outfile_ppi2_mspdb))
|
|
|
|
#%%===============================================================
|
|
#NOTE: duet here refers to the ensemble ppi2 values
|
|
|
|
###########################
|
|
# Read file: average ppi2 values
|
|
# or mcsm_normalised file
|
|
###########################
|
|
my_df_raw = merged_df3[, c("position", "mcsm_ppi2_scaled", "interface_dist")]
|
|
head(my_df_raw)
|
|
my_df_raw = my_df_raw[my_df_raw$interface_dist<10,]
|
|
my_df_raw$position
|
|
|
|
# avg by position on the SCALED values
|
|
my_df <- my_df_raw %>%
|
|
group_by(position) %>%
|
|
summarize(avg_ppi2_sc_pos = mean(mcsm_ppi2_scaled))
|
|
|
|
max(my_df$avg_ppi2_sc_pos)
|
|
min(my_df$avg_ppi2_sc_pos)
|
|
#============================================================
|
|
#############
|
|
# 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 ppi2 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 ppi2 values
|
|
# 3: replaced B-factors with mean ppi2 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 ppi2
|
|
#=============
|
|
hist(df_duet$b
|
|
, xlab = ""
|
|
, main = "Bfactor ppi2")
|
|
|
|
plot(density(df_duet$b)
|
|
, xlab = ""
|
|
, main = "Bfactor ppi2")
|
|
|
|
#=============
|
|
# Row 2 plots: original mean ppi2 values
|
|
# ppi2
|
|
#=============
|
|
|
|
#hist(my_df$averaged_duet
|
|
hist(my_df$avg_ppi2_sc_pos
|
|
, xlab = ""
|
|
, main = "mean ppi2 values")
|
|
|
|
#plot(density(my_df$averaged_duet)
|
|
plot(density(my_df$avg_ppi2_sc_pos)
|
|
, xlab = ""
|
|
, main = "mean ppi2 values")
|
|
|
|
#=========
|
|
# step 1_P1
|
|
#=========
|
|
df_duet$b = my_df$avg_ppi2_sc_pos[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
|
|
|
|
#=========
|
|
# 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: OUTPUT
|
|
#=========
|
|
cat(paste0("output file duet mean ppi2 pdb:", outfile_ppi2_mspdb))
|
|
write.pdb(my_pdb_duet, outfile_ppi2_mspdb)
|
|
|
|
# OUTPUT: position file
|
|
poscsvF = paste0(outdir_images, tolower(gene), "_ppi2_positions.csv")
|
|
cat(paste0("output file duet mean ppi2 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 ppi2 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)
|
|
#============================================ |