output from comb script & electrostatic mut changes calculated

This commit is contained in:
Tanushree Tunstall 2020-03-25 13:42:18 +00:00
parent 954eb88c45
commit d44ab57f5a
4 changed files with 250 additions and 167 deletions

View file

@ -1,6 +1,19 @@
#########################################################
# TASK: To combine mcsm and meta data with af and or
# This script doesn't output anything, but can do if needed.
# TASK: To combine mcsm and meta data with af and or files
# Input csv files:
# 1) mcsm output formatted
# 2) gene associated meta_data_with_AFandOR
# 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
#########################################################
getwd()
@ -10,7 +23,6 @@ getwd()
##########################################################
# Installing and loading required packages
##########################################################
source('Header_TT.R')
#require(data.table)
#require(arsenal)
@ -21,19 +33,23 @@ source('Header_TT.R')
# Read file: normalised file
# output of step 4 mcsm_pipeline
#################################
#%% variable assignment: input and output paths & filenames
drug = 'pyrazinamide'
gene = 'pncA'
gene_match = paste0(gene,'_p.')
cat(gene_match)
#===========
# data dir
#===========
datadir = paste0('~/git/Data')
#===========
# input
#===========
# infile1: mCSM data
#indir = '~/git/Data/pyrazinamide/input/processed/'
indir = paste0('~/git/Data', '/', drug, '/', 'output') # revised {TODO: change in mcsm pipeline}
indir = paste0(datadir, '/', drug, '/', 'output') # revised {TODO: change in mcsm pipeline}
in_filename = 'mcsm_complex1_normalised.csv'
infile = paste0(indir, '/', in_filename)
cat(paste0('Reading infile1: mCSM output file', ' ', infile) )
@ -105,8 +121,9 @@ changes = mcsm_data[which(mcsm_data$DUET_outcome != mcsm_data$Lig_outcome),]
dl_i = which(mcsm_data$DUET_outcome != mcsm_data$Lig_outcome)
ld_i = which(mcsm_data$Lig_outcome != mcsm_data$DUET_outcome)
cat('Identifying muts with opposite stability effects')
if(nrow(changes) == (table(mcsm_data$DUET_outcome != mcsm_data$Lig_outcome)[[2]]) & identical(dl_i,ld_i)) {
cat('PASS: muts with opposite effects on stability and affinity identified correctly'
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')
@ -134,6 +151,7 @@ mcsm_data = mcsm_data[order(mcsm_data$Mutationinformation),]
head(mcsm_data$Mutationinformation)
orig_col = ncol(mcsm_data)
# get freq count of positions and add to the df
setDT(mcsm_data)[, occurrence := .N, by = .(Position)]
@ -158,6 +176,20 @@ cat('Read mcsm_data file:'
, '\nNo.of rows: ', nrow(meta_with_afor)
, '\nNo. of cols:', ncol(meta_with_afor))
# counting NAs in AF, OR cols
if (identical(sum(is.na(meta_with_afor$OR))
, sum(is.na(meta_with_afor$pvalue))
, sum(is.na(meta_with_afor$AF)))){
cat('PASS: NA count match for OR, pvalue and AF\n')
na_count = sum(is.na(meta_with_afor$AF))
cat('No. of NAs: ', sum(is.na(meta_with_afor$OR)))
} else{
cat('FAIL: NA count mismatch'
, '\nNA in OR: ', sum(is.na(meta_with_afor$OR))
, '\nNA in pvalue: ', sum(is.na(meta_with_afor$pvalue))
, '\nNA in AF:', sum(is.na(meta_with_afor$AF)))
}
# clear variables
rm(in_filename_comb, infile_comb)
@ -172,15 +204,15 @@ head(meta_with_afor$Mutationinformation)
# 3: merging two dfs: with NA
###########################
# link col name = 'Mutationinforamtion'
head(mcsm_data$Mutationinformation)
head(meta_with_afor$Mutationinformation)
cat('Merging dfs with NAs: big df (1-many relationship b/w id & mut)'
,'\nlinking col: Mutationinforamtion'
,'\nfilename: merged_df2')
head(mcsm_data$Mutationinformation)
head(meta_with_afor$Mutationinformation)
#########
# merge 3a: meta data with mcsm
# merge 3a (merged_df2): meta data with mcsm
#########
merged_df2 = merge(x = meta_with_afor
,y = mcsm_data
@ -192,6 +224,8 @@ cat('Dim of merged_df2: '
, '\nNo. of cols: ', ncol(merged_df2))
head(merged_df2$Position)
# sanity check
cat('Checking nrows in merged_df2')
if(nrow(meta_with_afor) == nrow(merged_df2)){
cat('nrow(merged_df2) = nrow (gene associated metadata)'
,'\nExpected no. of rows: ',nrow(meta_with_afor)
@ -229,9 +263,9 @@ table(merged_df2$Position%in%merged_df2v2$Position)
rm(merged_df2v2)
#########
# merge 3b:remove duplicate mutation information
# merge 3b (merged_df3):remove duplicate mutation information
#########
cat('Merging dfs with NAs: small df (removing duplicate muts)'
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')
@ -244,8 +278,8 @@ cat('Merging dfs with NAs: small df (removing duplicate muts)'
merged_df3 = merged_df2[!duplicated(merged_df2$Mutationinformation),]
head(merged_df3$Position); tail(merged_df3$Position) # should be sorted
# sanity checks
# nrows of merged_df3 should be the same as the nrows of mcsm_data
# 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)
@ -256,41 +290,51 @@ if(nrow(mcsm_data) == nrow(merged_df3)){
, '\nNo. of rows merged_df3: ', nrow(merged_df3))
}
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# uncomment as necessary
# only need to run this if merged_df2v2 i.e non structural pos included
#mcsm = mcsm_data$Mutationinformation
#my_merged = merged_df3$Mutationinformation
# find the index where it differs
#diff_n = which(!my_merged%in%mcsm)
#check if it is indeed pos 186
#merged_df3[diff_n,]
# remove this entry
#merged_df3 = merged_df3[-diff_n,]]
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# 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')
#########
# merge 4a: same as merge 1 but excluding NA
#########
merged_df2_comp = merged_df2[!is.na(merged_df2$AF),]
#merged_df2_comp = merged_df2[!duplicated(merged_df2$Mutationinformation),]
cat('Dim of merged_df2_comp: '
, '\nNo. of rows: ', nrow(merged_df2_comp)
, '\nNo. of cols: ', ncol(merged_df2_comp))
# 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: remove duplicate mutation information
# merge 4b (merged_df3_comp): remove duplicate mutation information
#########
merged_df3_comp = merged_df2_comp[!duplicated(merged_df2_comp$Mutationinformation),]
@ -305,38 +349,65 @@ 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
#*********************
# write_output files
# output dir
#outdir = '~/git/Data/pyrazinamide/output/'
#uncomment as necessary
#FIXME
#out_filenames = c('merged_df2'
# , 'merged_df3'
# , 'meregd_df2_comp'
# , 'merged_df3_comp'
#)
# 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 files: '
# , '\nPath:', outdir)
#cat('Writing output file: '
# ,'\nFilename: ', out_filename
# ,'\nPath: ', outdir)
#for (i in out_filenames){
# print(i)
# print(get(i))
# outvar = get(i)
# print(outvar)
# outfile = paste0(outdir, '/', outvar, '.csv')
# cat('Writing output file:'
# ,'\nFilename: ', outfile
# ,'\n')
# write.csv(outvar, outfile)
# cat('Finished writing file:'
# ,'\nNo. of rows:', nrow(outvar)
# , '\nNo. of cols:', ncol(outvar))
#}
#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)))
#sapply(out_filenames, function(x) write.csv(x, 'x.csv'))
#%% 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, meta_with_afor, foo, drug, gene, gene_match, indir, merged_muts_u, meta_muts_u, na_count, orig_col, outdir)

View file

@ -1,8 +1,8 @@
#########################################################
# TASK: To combine mcsm and meta data with af and or
# by filtering for distance to ligand (<10Ang).
# This script doesn't output anything, but can do if needed.
# This script is sourced from other .R scripts for plotting
# This script doesn't output anything.
# This script is sourced from other .R scripts for plotting ligand plots
#########################################################
getwd()
setwd('~/git/LSHTM_analysis/mcsm_analysis/pyrazinamide/scripts/')