added plotting scripts from old run
This commit is contained in:
parent
d78626048c
commit
a448d9276b
2 changed files with 760 additions and 0 deletions
442
scripts/plotting/combining_two_df_FIXME.R
Normal file
442
scripts/plotting/combining_two_df_FIXME.R
Normal file
|
|
@ -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) <gene>_all_params.csv
|
||||||
|
# 2) <gene>_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: <gene>_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
|
||||||
|
|
||||||
318
scripts/plotting/replaceBfactor_pdb.R
Normal file
318
scripts/plotting/replaceBfactor_pdb.R
Normal file
|
|
@ -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???
|
||||||
|
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
||||||
|
|
||||||
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue