422 lines
14 KiB
R
422 lines
14 KiB
R
#!/usr/bin/env Rscript
|
|
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)
|
|
|
|
source("plotting_data.R")
|
|
|
|
# should return the following dfs, directories and variables
|
|
# my_df
|
|
# my_df_u
|
|
# my_df_u_lig
|
|
# dup_muts
|
|
|
|
cat(paste0("Directories imported:"
|
|
, "\ndatadir:", datadir
|
|
, "\nindir:", indir
|
|
, "\noutdir:", outdir
|
|
, "\nplotdir:", plotdir))
|
|
|
|
cat(paste0("Variables imported:"
|
|
, "\ndrug:", drug
|
|
, "\ngene:", gene
|
|
, "\ngene_match:", gene_match
|
|
, "\nLength of upos:", length(upos)
|
|
, "\nAngstrom symbol:", angstroms_symbol))
|
|
|
|
# clear excess variable
|
|
rm(my_df, upos, dup_muts, my_df_u_lig)
|
|
#========================================================
|
|
|
|
#========================================================
|
|
#%% 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")
|
|
plotdir = paste0("~/git/Data", "/", drug, "/output/plots")
|
|
#===========
|
|
# input
|
|
#===========
|
|
#in_file1: output of plotting_data.R
|
|
|
|
|
|
# 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)
|
|
|
|
|
|
#%%===============================================================
|
|
table(my_df_u$duet_outcome); sum(table(my_df_u$duet_outcome) )
|
|
|
|
# spelling Correction 1: DUET incase American spelling needed!
|
|
#my_df_u$duet_outcome[my_df_u$duet_outcome=="Stabilising"] <- "Stabilizing"
|
|
#my_df_u$duet_outcome[my_df_u$duet_outcome=="Destabilising"] <- "Destabilizing"
|
|
|
|
|
|
# spelling Correction 2: Ligand incase American spelling needed!
|
|
table(my_df_u$ligand_outcome); sum(table(my_df_u$ligand_outcome) )
|
|
#my_df_u$ligand_outcome[my_df_u$ligand_outcome=="Stabilising"] <- "Stabilizing"
|
|
#my_df_u$ligand_outcome[my_df_u$ligand_outcome=="Destabilising"] <- "Destabilizing"
|
|
|
|
|
|
# muts with opposing effects on protomer and ligand stability
|
|
table(my_df_u$duet_outcome != my_df_u$ligand_outcome)
|
|
changes = my_df_u[which(my_df_u$duet_outcome != my_df_u$ligand_outcome),]
|
|
|
|
# sanity check: redundant, but uber cautious!
|
|
dl_i = which(my_df_u$duet_outcome != my_df_u$ligand_outcome)
|
|
ld_i = which(my_df_u$ligand_outcome != my_df_u$duet_outcome)
|
|
|
|
cat("Identifying muts with opposite stability effects")
|
|
if(nrow(changes) == (table(my_df_u$duet_outcome != my_df_u$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(my_df_u, function(y) sum(length(which(is.na(y))))); na_count
|
|
df_ncols = ncol(my_df_u)
|
|
|
|
###########################
|
|
# 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(my_df_u$or_mychisq))
|
|
, sum(is.na(my_df_u$pval_fisher))
|
|
, sum(is.na(my_df_u$af)))){
|
|
cat("\nPASS: NA count match for OR, pvalue and AF\n")
|
|
na_count = sum(is.na(my_df_u$af))
|
|
cat("\nNo. of NAs: ", sum(is.na(my_df_u$or_mychisq)))
|
|
} else{
|
|
cat("\nFAIL: NA count mismatch"
|
|
, "\nNA in OR: ", sum(is.na(my_df_u$or_mychisq))
|
|
, "\nNA in pvalue: ", sum(is.na(my_df_u$pval_fisher))
|
|
, "\nNA in AF:", sum(is.na(my_df_u$af)))
|
|
}
|
|
|
|
|
|
if (identical(sum(is.na(my_df_u$or_kin))
|
|
, sum(is.na(my_df_u$pwald_kin))
|
|
, sum(is.na(my_df_u$af_kin)))){
|
|
cat("\nPASS: NA count match for OR, pvalue and AF\n from Kinship matrix calculations")
|
|
na_count = sum(is.na(my_df_u$af_kin))
|
|
cat("\nNo. of NAs: ", sum(is.na(my_df_u$or_kin)))
|
|
} else{
|
|
cat("\nFAIL: NA count mismatch"
|
|
, "\nNA in OR: ", sum(is.na(my_df_u$or_kin))
|
|
, "\nNA in pvalue: ", sum(is.na(my_df_u$pwald_kin))
|
|
, "\nNA in AF:", sum(is.na(my_df_u$af_kin)))
|
|
}
|
|
|
|
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
# clear variables
|
|
rm(in_filename_gene_metadata, infile_gene_metadata)
|
|
|
|
str(gene_metadata)
|
|
|
|
# sort by position (same as my_df)
|
|
# 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(my_df_u$mutationinformation)
|
|
head(gene_metadata$mutationinformation)
|
|
|
|
# Find common columns b/w two df
|
|
# FIXME: mutation has empty cell for some muts
|
|
merging_cols = intersect(colnames(my_df_u), 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)
|
|
|
|
# important checks!
|
|
table(nchar(my_df_u$mutationinformation))
|
|
table(nchar(my_df_u$wild_type))
|
|
table(nchar(my_df_u$mutant_type))
|
|
table(nchar(my_df_u$position))
|
|
|
|
#=============
|
|
# merged_df2: gene_metadata + my_df
|
|
#==============
|
|
# all.y because x might contain non-structural positions!
|
|
merged_df2 = merge(x = gene_metadata
|
|
, y = my_df_u
|
|
, by = merging_cols
|
|
, all.y = T)
|
|
|
|
cat("Dim of merged_df2: ", dim(merged_df2))
|
|
head(merged_df2$position)
|
|
|
|
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
# FIXME: count how many unique muts have entries in meta data
|
|
# should PASS
|
|
cat("Checking nrows in merged_df2")
|
|
if(nrow(gene_metadata) == nrow(merged_df2)){
|
|
cat("PASS: nrow(merged_df2) = nrow (gene associated gene_metadata)"
|
|
,"\nExpected no. of rows: ",nrow(gene_metadata)
|
|
,"\nGot no. of rows: ", nrow(merged_df2))
|
|
} else{
|
|
cat("FAIL: 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 = my_df_u
|
|
, by = merging_cols
|
|
, all = T)
|
|
|
|
merged_df2v2 = merge(x = gene_metadata
|
|
, y = my_df_u
|
|
, by = merging_cols
|
|
, all.x = T)
|
|
#!=!=!=!=!=!=!=!
|
|
#identical(merged_df2, merged_df2v2)
|
|
|
|
nrow(merged_df2[merged_df2$position==186,])
|
|
#!=!=!=!=!=!=!=!
|
|
|
|
# should be False
|
|
identical(merged_df2, merged_df2v2)
|
|
table(merged_df2$position%in%merged_df2v2$position)
|
|
|
|
#!!!!!!!!!!! check why these differ
|
|
|
|
#########
|
|
# merge 3b (merged_df3):remove duplicated mutations
|
|
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(my_df_u) == nrow(merged_df3)){
|
|
cat("PASS: No. of rows match with my_df"
|
|
,"\nExpected no. of rows: ", nrow(my_df_u)
|
|
,"\nGot no. of rows: ", nrow(merged_df3))
|
|
} else {
|
|
cat("FAIL: No. of rows mismatch"
|
|
, "\nNo. of rows my_df: ", nrow(my_df_u)
|
|
, "\nNo. of rows merged_df3: ", nrow(merged_df3))
|
|
quit()
|
|
}
|
|
|
|
# counting NAs in AF, OR cols in merged_df3
|
|
# this is because mcsm has no AF, OR cols,
|
|
# so you cannot count NAs
|
|
if (identical(sum(is.na(merged_df3$or_kin))
|
|
, sum(is.na(merged_df3$pwald_kin))
|
|
, sum(is.na(merged_df3$af_kin)))){
|
|
cat("PASS: NA count match for OR, pvalue and AF\n")
|
|
na_count_df3 = sum(is.na(merged_df3$af_kin))
|
|
cat("No. of NAs: ", sum(is.na(merged_df3$or_kin)))
|
|
} else{
|
|
cat("FAIL: NA count mismatch"
|
|
, "\nNA in OR: ", sum(is.na(merged_df3$or_kin))
|
|
, "\nNA in pvalue: ", sum(is.na(merged_df3$pwald_kin))
|
|
, "\nNA in AF:", sum(is.na(merged_df3$af_kin)))
|
|
}
|
|
|
|
# check if the same or and afs are missing for
|
|
if ( identical( which(is.na(merged_df2$or_mychisq)), which(is.na(merged_df2$or_kin)))
|
|
&& identical( which(is.na(merged_df2$af)), which(is.na(merged_df2$af_kin)))
|
|
&& identical( which(is.na(merged_df2$pval_fisher)), which(is.na(merged_df2$pwald_kin))) ){
|
|
cat('PASS: Indices match for mychisq and kin ors missing values')
|
|
} else{
|
|
cat('Index mismatch: mychisq and kin ors missing indices match')
|
|
quit()
|
|
}
|
|
|
|
###########################
|
|
# 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")
|
|
|
|
if ( identical( which(is.na(merged_df2$af)), which(is.na(merged_df2$af_kin))) ){
|
|
print('mychisq and kin ors missing indices match. Procedding with omitting NAs')
|
|
na_count_df2 = sum(is.na(merged_df2$af))
|
|
merged_df2_comp = merged_df2[!is.na(merged_df2$af),]
|
|
# sanity check: no +-1 gymnastics
|
|
cat("Checking nrows in merged_df2_comp")
|
|
if(nrow(merged_df2_comp) == (nrow(merged_df2) - na_count_df2)){
|
|
cat("\nPASS: No. of rows match"
|
|
,"\nDim of merged_df2_comp: "
|
|
,"\nExpected no. of rows: ", nrow(merged_df2) - na_count_df2
|
|
, "\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_df2
|
|
,"\nGot no. of rows: ", nrow(merged_df2_comp))
|
|
}
|
|
}else{
|
|
print('Index mismatch for mychisq and kin ors. Aborting NA ommission')
|
|
}
|
|
|
|
#########
|
|
# merge 4b (merged_df3_comp): remove duplicate mutation information
|
|
#########
|
|
if ( identical( which(is.na(merged_df3$af)), which(is.na(merged_df3$af_kin))) ){
|
|
print('mychisq and kin ors missing indices match. Procedding with omitting NAs')
|
|
na_count_df3 = sum(is.na(merged_df3$af))
|
|
#merged_df3_comp = merged_df3_comp[!duplicated(merged_df3_comp$mutationinformation),] # a way
|
|
merged_df3_comp = merged_df3[!is.na(merged_df3$af),] # another way
|
|
cat("Checking nrows in merged_df3_comp")
|
|
if(nrow(merged_df3_comp) == (nrow(merged_df3) - na_count_df3)){
|
|
cat("\nPASS: No. of rows match"
|
|
,"\nDim of merged_df3_comp: "
|
|
,"\nExpected no. of rows: ", nrow(merged_df3) - na_count_df3
|
|
, "\nNo. of rows: ", nrow(merged_df3_comp)
|
|
, "\nNo. of cols: ", ncol(merged_df3_comp))
|
|
}else{
|
|
cat("FAIL: No. of rows mismatch"
|
|
,"\nExpected no. of rows: ", nrow(merged_df3) - na_count_df3
|
|
,"\nGot no. of rows: ", nrow(merged_df3_comp))
|
|
}
|
|
} else{
|
|
print('Index mismatch for mychisq and kin ors. Aborting NA ommission')
|
|
}
|
|
|
|
# alternate way of deriving merged_df3_comp
|
|
foo = merged_df3[!is.na(merged_df3$af),]
|
|
bar = merged_df3_comp[!duplicated(merged_df3_comp$mutationinformation),]
|
|
# compare dfs: foo and merged_df3_com
|
|
all.equal(foo, bar)
|
|
#summary(comparedf(foo, bar))
|
|
|
|
#=============== end of combining df
|
|
#==============================================================
|
|
#################
|
|
# OPTIONAL: write ALL 4 output files
|
|
#################
|
|
#outvars = c("merged_df2"
|
|
# , "merged_df3"
|
|
# , "merged_df2_comp"
|
|
# , "merged_df3_comp")
|
|
|
|
#cat("Writing output files: "
|
|
# , "\nPath:", outdir)
|
|
|
|
#for (i in outvars){
|
|
# out_filename = paste0(i, ".csv")
|
|
# outfile = paste0(outdir, "/", out_filename)
|
|
# 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")
|
|
#}
|
|
|
|
#*************************
|
|
# clear variables
|
|
rm(foo, bar, gene_metadata
|
|
, in_filename_params, infile_params, merging_cols
|
|
, merged_df2v2, merged_df2v3)
|
|
|
|
#============================= end of script
|
|
|