From 954eb88c452122164eb4b3aa71d9bf7764ecc603 Mon Sep 17 00:00:00 2001 From: Tanushree Tunstall Date: Tue, 24 Mar 2020 18:28:52 +0000 Subject: [PATCH] updated combining df scripts for duet & lig --- .../pyrazinamide/scripts/combining_two_df.R | 326 ++++++++++-------- .../scripts/combining_two_df_lig.R | 321 ++++++++--------- meta_data_analysis/AF_and_OR_calcs.R | 153 ++++++-- 3 files changed, 458 insertions(+), 342 deletions(-) diff --git a/mcsm_analysis/pyrazinamide/scripts/combining_two_df.R b/mcsm_analysis/pyrazinamide/scripts/combining_two_df.R index d0e13cd..47302b9 100644 --- a/mcsm_analysis/pyrazinamide/scripts/combining_two_df.R +++ b/mcsm_analysis/pyrazinamide/scripts/combining_two_df.R @@ -1,17 +1,17 @@ +######################################################### +# TASK: To combine mcsm and meta data with af and or +# This script doesn't output anything, but can do if needed. +# This script is sourced from other .R scripts for plotting +######################################################### getwd() -setwd("~/git/LSHTM_analysis/mcsm_analysis/pyrazinamide/scripts/") +setwd('~/git/LSHTM_analysis/mcsm_analysis/pyrazinamide/scripts/') getwd() -######################################################### -# TASK: To combine mcsm and meta data with af and or -######################################################### +########################################################## +# Installing and loading required packages +########################################################## - -######################################################################## -# Installing and loading required packages # -######################################################################## - -#source("Header_TT.R") +source('Header_TT.R') #require(data.table) #require(arsenal) #require(compare) @@ -22,14 +22,58 @@ getwd() # output of step 4 mcsm_pipeline ################################# -inDir = "~/git/Data/pyrazinamide/input/processed/" -inFile = paste0(inDir, "mcsm_complex1_normalised.csv"); inFile +#%% variable assignment: input and output paths & filenames +drug = 'pyrazinamide' +gene = 'pncA' +gene_match = paste0(gene,'_p.') +cat(gene_match) -mcsm_data = read.csv(inFile +#=========== +# input +#=========== +# infile1: mCSM data +#indir = '~/git/Data/pyrazinamide/input/processed/' +indir = paste0('~/git/Data', '/', 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) ) + +# infile2: gene associated meta data combined with AF and OR +#indir: same as above +in_filename_comb = paste0(tolower(gene), '_meta_data_with_AFandOR.csv') +infile_comb = paste0(indir, '/', in_filename_comb) +cat(paste0('Reading infile2: gene associated combined metadata:', infile_comb)) + +#=========== +# output +#=========== +# Uncomment if and when required to output +outdir = paste0('~/git/Data', '/', drug, '/', 'output') #same as indir +cat('Output dir: ', outdir) +#out_filename = paste0(tolower(gene), 'XXX') +#outfile = paste0(outdir, '/', out_filename) +#cat(paste0('Output file with full path:', outfile)) +#%% end of variable assignment for input and output files + +################################# +# Read file: normalised file +# output of step 4 mcsm_pipeline +################################# +cat('Reading mcsm_data:' + , '\nindir: ', indir + , '\ninfile_comb: ', in_filename) + +mcsm_data = read.csv(infile , row.names = 1 , stringsAsFactors = F - , header = T) -rm(inDir, inFile) + , header = T) + +cat('Read mcsm_data file:' + , '\nNo.of rows: ', nrow(mcsm_data) + , '\nNo. of cols:', ncol(mcsm_data)) + +# clear variables +rm(in_filename, infile) str(mcsm_data) @@ -53,6 +97,35 @@ mcsm_data$Lig_outcome[mcsm_data$Lig_outcome=='Destabilizing'] <- 'Destabilising' table(mcsm_data$Lig_outcome); sum(table(mcsm_data$Lig_outcome) ) head(mcsm_data$Lig_outcome); tail(mcsm_data$Lig_outcome) +# muts with opposing effects on protomer and ligand stability +table(mcsm_data$DUET_outcome != mcsm_data$Lig_outcome) +changes = mcsm_data[which(mcsm_data$DUET_outcome != mcsm_data$Lig_outcome),] + +# sanity check: redundant, but uber cautious! +dl_i = which(mcsm_data$DUET_outcome != mcsm_data$Lig_outcome) +ld_i = which(mcsm_data$Lig_outcome != mcsm_data$DUET_outcome) + +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' + , '\nNo. of such muts: ', nrow(changes)) +}else { + cat('FAIL: unsuccessful in extracting muts with changed stability effects') +} + +#*************************** +# write file: changed muts +out_filename = 'muts_opp_effects.csv' +outfile = paste0(outdir, '/', out_filename) +cat('Writing file for muts with opp effects:' + , '\nFilename: ', outfile + , '\nPath: ', outdir) + +write.csv(changes, outfile) +#**************************** +# clear variables +rm(out_filename, outfile) +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 @@ -60,23 +133,33 @@ na_count = sapply(mcsm_data, function(y) sum(length(which(is.na(y))))); na_count 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)] +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: ', orig_col) + pos_count_check = data.frame(mcsm_data$Position, mcsm_data$occurrence) ########################### # 2: Read file: meta data with AFandOR ########################### +cat('Reading combined meta data and AFandOR file:' + , '\nindir: ', indir + , '\ninfile_comb: ', in_filename_comb) -inDir = "~/git/Data/pyrazinamide/input/processed/" -inFile2 = paste0(inDir, "meta_data_with_AFandOR.csv"); inFile2 - -meta_with_afor <- read.csv(inFile2 +meta_with_afor <- read.csv(infile_comb , stringsAsFactors = F , header = T) -rm(inDir, inFile2) +cat('Read mcsm_data file:' + , '\nNo.of rows: ', nrow(meta_with_afor) + , '\nNo. of cols:', ncol(meta_with_afor)) + +# clear variables +rm(in_filename_comb, infile_comb) str(meta_with_afor) @@ -85,113 +168,45 @@ head(meta_with_afor$Mutationinformation) meta_with_afor = meta_with_afor[order(meta_with_afor$Mutationinformation),] head(meta_with_afor$Mutationinformation) -# sanity check: should be True for all the mentioned columns -#is.numeric(meta_with_afor$OR) -na_var = c("AF", "OR", "pvalue", "logor", "neglog10pvalue") - -c1 = NULL -for (i in na_var){ - print(i) - c0 = is.numeric(meta_with_afor[,i]) - c1 = c(c0, c1) - if ( all(c1) ){ - print("Sanity check passed: These are all numeric cols") - } else{ - print("Error: Please check your respective data types") - } -} - -# If OR, and P value are not numeric, then convert to numeric and then count -# else they will say 0 -na_count = sapply(meta_with_afor, function(y) sum(length(which(is.na(y))))); na_count -str(na_count) - -# compare if the No of "NA" are the same for all these cols -na_len = NULL -for (i in na_var){ - temp = na_count[[i]] - na_len = c(na_len, temp) -} - -# extract how many NAs: -# should be all TRUE -# should be a single number since -# all the cols should have "equal" and "same" no. of NAs - -my_nrows = NULL -for ( i in 1: (length(na_len)-1) ){ - #print(compare(na_len[i]), na_len[i+1]) - c = compare(na_len[i], na_len[i+1]) - if ( c$result ) { - my_nrows = na_len[i] } - else { - print("Error: Please check your numbers") - } -} - -my_nrows - -#=#=#=#=#=#=#=#=# -# COMMENT: AF, OR, pvalue, logor and neglog10pvalue -# these are the same 7 ones -#=#=#=#=#=#=#=#=# - -# sanity check -#which(is.na(meta_with_afor$OR)) - -# initialise an empty df with nrows as extracted above -na_count_df = data.frame(matrix(vector(mode = 'numeric' -# , length = length(na_var) - ) - , nrow = my_nrows -# , ncol = length(na_var) - )) - -# populate the df with the indices of the cols that are NA -for (i in na_var){ - print(i) - na_i = which(is.na(meta_with_afor[i])) - na_count_df = cbind(na_count_df, na_i) - colnames(na_count_df)[which(na_var == i)] <- i -} - -# Now compare these indices to ensure these are the same -c2 = NULL -for ( i in 1: ( length(na_count_df)-1 ) ) { -# print(na_count_df[i] == na_count_df[i+1]) - c1 = identical(na_count_df[[i]], na_count_df[[i+1]]) - c2 = c(c1, c2) - if ( all(c2) ) { - print("Sanity check passed: The indices for AF, OR, etc are all the same") - } else { - print ("Error: Please check indices which are NA") - } -} - -rm( c, c0, c1, c2, i, my_nrows - , na_count, na_i, na_len - , na_var, temp - , na_count_df - , pos_count_check ) - ########################### -# 3:merging two dfs: with NA +# 3: merging two dfs: with NA ########################### +# link col name = 'Mutationinforamtion' +cat('Merging dfs with NAs: big df (1-many relationship b/w id & mut)' + ,'\nlinking col: Mutationinforamtion' + ,'\nfilename: merged_df2') -# link col name = Mutationinforamtion head(mcsm_data$Mutationinformation) head(meta_with_afor$Mutationinformation) ######### -# merge 1a: meta data with mcsm +# merge 3a: meta data with mcsm ######### merged_df2 = merge(x = meta_with_afor ,y = mcsm_data - , by = "Mutationinformation" + , by = 'Mutationinformation' , all.y = T) +cat('Dim of merged_df2: ' + , '\nNo. of rows: ', nrow(merged_df2) + , '\nNo. of cols: ', ncol(merged_df2)) head(merged_df2$Position) +if(nrow(meta_with_afor) == nrow(merged_df2)){ + cat('nrow(merged_df2) = nrow (gene associated metadata)' + ,'\nExpected no. of rows: ',nrow(meta_with_afor) + ,'\nGot no. of rows: ', nrow(merged_df2)) +} else{ + cat('nrow(merged_df2)!= nrow(gene associated metadata)' + , '\nExpected no. of rows after merge: ', nrow(meta_with_afor) + , '\nGot no. of rows: ', nrow(merged_df2) + , '\nFinding discrepancy') + merged_muts_u = unique(merged_df2$Mutationinformation) + meta_muts_u = unique(meta_with_afor$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),] @@ -199,12 +214,12 @@ head(merged_df2$Position) merged_df2v2 = merge(x = meta_with_afor ,y = mcsm_data - , by = "Mutationinformation" + , by = 'Mutationinformation' , 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 with mutation +# but 186 is associated with mutation #!=!=!=!=!=!=!=! # should be False @@ -214,8 +229,12 @@ table(merged_df2$Position%in%merged_df2v2$Position) rm(merged_df2v2) ######### -# merge 1b:remove duplicate mutation information +# merge 3b:remove duplicate mutation information ######### +cat('Merging dfs with NAs: small df (removing duplicate muts)' + ,'\nCannot trust lineage info from this' + ,'\nlinking col: Mutationinforamtion' + ,'\nfilename: merged_df3') #==#=#=#=#=#=# # Cannot trust lineage, country from this df as the same mutation @@ -228,9 +247,13 @@ 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 if(nrow(mcsm_data) == nrow(merged_df3)){ - print("sanity check: Passed") + cat('PASS: No. of rows match with mcsm_data' + ,'\nExpected no. of rows: ', nrow(mcsm_data) + ,'\nGot no. of rows: ', nrow(merged_df3)) } else { - print("Error!: check data, nrows is not as expected") + cat('FAIL: No. of rows mismatch' + , '\nNo. of rows mcsm_data: ', nrow(mcsm_data) + , '\nNo. of rows merged_df3: ', nrow(merged_df3)) } #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -250,19 +273,31 @@ if(nrow(mcsm_data) == nrow(merged_df3)){ #!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ########################### -# 3b :merging two dfs: without NA +# 4: merging two dfs: without NA ########################### +cat('Merging dfs without any NAs: big df (1-many relationship b/w id & mut)' + ,'\nlinking col: Mutationinforamtion' + ,'\nfilename: merged_df2_comp') ######### -# merge 2a:same as merge 1 but excluding NA +# 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)) ######### -# merge 2b: remove duplicate mutation information +# merge 4b: 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 @@ -271,29 +306,40 @@ all.equal(foo, merged_df3) summary(comparedf(foo, merged_df3)) #=============== end of combining df -#clear variables -rm(mcsm_data - , meta_with_afor - , foo) - -#rm(diff_n, my_merged, mcsm) - -#===================== +#********************* # write_output files -#===================== # output dir -outDir = "~/git/Data/pyrazinamide/output/" -getwd() +#outdir = '~/git/Data/pyrazinamide/output/' +#uncomment as necessary +#FIXME +#out_filenames = c('merged_df2' +# , 'merged_df3' +# , 'meregd_df2_comp' +# , 'merged_df3_comp' +#) -outFile1 = paste0(outDir, "merged_df3.csv"); outFile1 -#write.csv(merged_df3, outFile1) +#cat('Writing output files: ' +# , '\nPath:', outdir) -#outFile2 = paste0(outDir, "merged_df3_comp.csv"); outFile2 -#write.csv(merged_df3_comp, outFile2) +#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)) +#} -rm(outDir - , outFile1 -# , outFile2 -) +#sapply(out_filenames, function(x) write.csv(x, 'x.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) +rm(pos_count_check) #============================= end of script diff --git a/mcsm_analysis/pyrazinamide/scripts/combining_two_df_lig.R b/mcsm_analysis/pyrazinamide/scripts/combining_two_df_lig.R index 1dccd65..14867c9 100644 --- a/mcsm_analysis/pyrazinamide/scripts/combining_two_df_lig.R +++ b/mcsm_analysis/pyrazinamide/scripts/combining_two_df_lig.R @@ -1,17 +1,18 @@ -getwd() -setwd("~/git/LSHTM_analysis/mcsm_analysis/pyrazinamide/scripts/") -getwd() - ######################################################### # TASK: To combine mcsm and meta data with af and or -# by filtering for distance to ligand (<10Ang) +# 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 ######################################################### +getwd() +setwd('~/git/LSHTM_analysis/mcsm_analysis/pyrazinamide/scripts/') +getwd() -######################################################### -# Installing and loading required packages -######################################################### +########################################################## +# Installing and loading required packages +########################################################## -#source("Header_TT.R") +source('Header_TT.R') #require(data.table) #require(arsenal) #require(compare) @@ -22,14 +23,58 @@ getwd() # output of step 4 mcsm_pipeline ################################# -inDir = "~/git/Data/pyrazinamide/input/processed/" -inFile = paste0(inDir, "mcsm_complex1_normalised.csv"); inFile +#%% variable assignment: input and output paths & filenames +drug = 'pyrazinamide' +gene = 'pncA' +gene_match = paste0(gene,'_p.') +cat(gene_match) -mcsm_data = read.csv(inFile +#=========== +# input +#=========== +# infile1: mCSM data +#indir = '~/git/Data/pyrazinamide/input/processed/' +indir = paste0('~/git/Data', '/', 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) ) + +# infile2: gene associated meta data combined with AF and OR +#indir: same as above +in_filename_comb = paste0(tolower(gene), '_meta_data_with_AFandOR.csv') +infile_comb = paste0(indir, '/', in_filename_comb) +cat(paste0('Reading infile2: gene associated combined metadata:', infile_comb)) + +#=========== +# output +#=========== +# Uncomment if and when required to output +outdir = paste0('~/git/Data', '/', drug, '/', 'output') #same as indir +cat('Output dir: ', outdir) +#out_filename = paste0(tolower(gene), 'XXX') +#outfile = paste0(outdir, '/', out_filename) +#cat(paste0('Output file with full path:', outfile)) +#%% end of variable assignment for input and output files + +################################# +# Read file: normalised file +# output of step 4 mcsm_pipeline +################################# +cat('Reading mcsm_data:' + , '\nindir: ', indir + , '\ninfile_comb: ', in_filename) + +mcsm_data = read.csv(infile , row.names = 1 , stringsAsFactors = F - , header = T) -rm(inDir, inFile) + , header = T) + +cat('Read mcsm_data file:' + , '\nNo.of rows: ', nrow(mcsm_data) + , '\nNo. of cols:', ncol(mcsm_data)) + +# clear variables +rm(in_filename, infile) str(mcsm_data) @@ -39,7 +84,7 @@ table(mcsm_data$DUET_outcome); sum(table(mcsm_data$DUET_outcome) ) mcsm_data$DUET_outcome[mcsm_data$DUET_outcome=='Stabilizing'] <- 'Stabilising' mcsm_data$DUET_outcome[mcsm_data$DUET_outcome=='Destabilizing'] <- 'Destabilising' -# checks +# 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) @@ -53,6 +98,10 @@ mcsm_data$Lig_outcome[mcsm_data$Lig_outcome=='Destabilizing'] <- 'Destabilising' table(mcsm_data$Lig_outcome); sum(table(mcsm_data$Lig_outcome) ) head(mcsm_data$Lig_outcome); tail(mcsm_data$Lig_outcome) +# muts with opposing effects on protomer and ligand stability +# excluded from here as it is redundant. +# check 'combining_two_df.R' to refer if required. + ########################### !!! only for mcsm_lig # 4: Filter/subset data # Lig plots < 10Ang @@ -82,57 +131,60 @@ length(unique(mcsm_data2$Mutationinformation)) # count Destabilisinga and stabilising table(mcsm_data2$Lig_outcome) #{RESULT: no of mutations within 10Ang} -#<<<<<<<<<<<<<<<<<<<<<<<<<<< -# REASSIGNMENT: so as not to alter the script -mcsm_data = mcsm_data2 -#<<<<<<<<<<<<<<<<<<<<<<<<<<< - ############################# # Extra sanity check: # for mcsm_lig ONLY # Dis_lig_Ang should be <10 ############################# -if (max(mcsm_data$Dis_lig_Ang) < 10){ +if (max(mcsm_data2$Dis_lig_Ang) < 10){ print ("Sanity check passed: lig data is <10Ang") }else{ print ("Error: data should be filtered to be within 10Ang") } +#!!!!!!!!!!!!!!!!!!!!! +# REASSIGNMENT: so as not to alter the script +mcsm_data = mcsm_data2 +#!!!!!!!!!!!!!!!!!!!!! # clear variables rm(mcsm_data2) # count na in each column na_count = sapply(mcsm_data, function(y) sum(length(which(is.na(y))))); na_count -head(mcsm_data$Mutationinformation) -mcsm_data[mcsm_data$Mutationinformation=="Q10P",] -mcsm_data[mcsm_data$Mutationinformation=="L4S",] - # sort by Mutationinformation mcsm_data = mcsm_data[order(mcsm_data$Mutationinformation),] head(mcsm_data$Mutationinformation) -# check -mcsm_data[grep("Q10P", mcsm_data$Mutationinformation),] -mcsm_data[grep("A102T", 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)] +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: ', orig_col) pos_count_check = data.frame(mcsm_data$Position, mcsm_data$occurrence) ########################### # 2: Read file: meta data with AFandOR ########################### +cat('Reading combined meta data and AFandOR file:' + , '\nindir: ', indir + , '\ninfile_comb: ', in_filename_comb) -inDir = "~/git/Data/pyrazinamide/input/processed/" -inFile2 = paste0(inDir, "meta_data_with_AFandOR.csv"); inFile2 - -meta_with_afor <- read.csv(inFile2 +meta_with_afor <- read.csv(infile_comb , stringsAsFactors = F , header = T) +cat('Read mcsm_data file:' + , '\nNo.of rows: ', nrow(meta_with_afor) + , '\nNo. of cols:', ncol(meta_with_afor)) + +# clear variables +rm(in_filename_comb, infile_comb) + str(meta_with_afor) # sort by Mutationinformation @@ -140,114 +192,45 @@ head(meta_with_afor$Mutationinformation) meta_with_afor = meta_with_afor[order(meta_with_afor$Mutationinformation),] head(meta_with_afor$Mutationinformation) -# sanity check: should be True for all the mentioned columns -#is.numeric(meta_with_afor$OR) -na_var = c("AF", "OR", "pvalue", "logor", "neglog10pvalue") - -c1 = NULL -for (i in na_var){ - print(i) - c0 = is.numeric(meta_with_afor[,i]) - c1 = c(c0, c1) - if ( all(c1) ){ - print("Sanity check passed: These are all numeric cols") - } else{ - print("Error: Please check your respective data types") - } -} - -# If OR, and P value are not numeric, then convert to numeric and then count -# else they will say 0 - -# NOW count na in each column: if you did it before, then -# OR and Pvalue column would say 0 na since these were not numeric -na_count = sapply(meta_with_afor, function(y) sum(length(which(is.na(y))))); na_count -str(na_count) - -# compare if the No of "NA" are the same for all these cols -na_len = NULL -na_var = c("AF", "OR", "pvalue", "logor", "neglog10pvalue") -for (i in na_var){ - temp = na_count[[i]] - na_len = c(na_len, temp) -} - -my_nrows = NULL - -for ( i in 1: (length(na_len)-1) ){ - #print(compare(na_len[i]), na_len[i+1]) - c = compare(na_len[i], na_len[i+1]) - if ( c$result ) { - my_nrows = na_len[i] } - else { - print("Error: Please check your numbers") - } -} - -my_nrows - -#=#=#=#=#=#=#=#=# -# COMMENT: AF, OR, pvalue, logor and neglog10pvalue -# all have 81 NA, with pyrazinamide with 960 -# and these are the same 7 ones -#=#=#=#=#=#=#=#=# - -# sanity check -#which(is.na(meta_with_afor$OR)) - -# initialise an empty df with nrows as extracted above -na_count_df = data.frame(matrix(vector(mode = 'numeric' -# , length = length(na_var) - ) - , nrow = my_nrows -# , ncol = length(na_var) - )) - -# populate the df with the indices of the cols that are NA -for (i in na_var){ - print(i) - na_i = which(is.na(meta_with_afor[i])) - na_count_df = cbind(na_count_df, na_i) - colnames(na_count_df)[which(na_var == i)] <- i -} - -# Now compare these indices to ensure these are the same -c2 = NULL -for ( i in 1: ( length(na_count_df)-1 ) ) { - # print(na_count_df[i] == na_count_df[i+1]) - c1 = identical(na_count_df[[i]], na_count_df[[i+1]]) - c2 = c(c1, c2) - if ( all(c2) ) { - print("Sanity check passed: The indices for AF, OR, etc are all the same") - } else { - print ("Error: Please check indices which are NA") - } -} - -rm( c, c1, c2, i, my_nrows - , na_count, na_i, na_len - , na_var, temp - , na_count_df - , pos_count_check ) - ########################### -# 3:merging two dfs: with NA +# 3: merging two dfs: with NA ########################### +# link col name = 'Mutationinforamtion' +cat('Merging dfs with NAs: big df (1-many relationship b/w id & mut)' + ,'\nlinking col: Mutationinforamtion' + ,'\nfilename: merged_df2') -# link col name = Mutationinforamtion head(mcsm_data$Mutationinformation) head(meta_with_afor$Mutationinformation) ######### -# merge 1a: meta data with mcsm +# merge 3a: meta data with mcsm ######### merged_df2 = merge(x = meta_with_afor - , y = mcsm_data - , by = "Mutationinformation" + ,y = mcsm_data + , by = 'Mutationinformation' , all.y = T) +cat('Dim of merged_df2: ' + , '\nNo. of rows: ', nrow(merged_df2) + , '\nNo. of cols: ', ncol(merged_df2)) head(merged_df2$Position) +if(nrow(meta_with_afor) == nrow(merged_df2)){ + cat('nrow(merged_df2) = nrow (gene associated metadata)' + ,'\nExpected no. of rows: ',nrow(meta_with_afor) + ,'\nGot no. of rows: ', nrow(merged_df2)) +} else{ + cat('nrow(merged_df2)!= nrow(gene associated metadata)' + , '\nExpected no. of rows after merge: ', nrow(meta_with_afor) + , '\nGot no. of rows: ', nrow(merged_df2) + , '\nFinding discrepancy') + merged_muts_u = unique(merged_df2$Mutationinformation) + meta_muts_u = unique(meta_with_afor$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),] @@ -255,13 +238,12 @@ head(merged_df2$Position) merged_df2v2 = merge(x = meta_with_afor ,y = mcsm_data - , by = "Mutationinformation" + , by = 'Mutationinformation' , 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 with mutation +# but 186 is associated with mutation #!=!=!=!=!=!=!=! # should be False @@ -271,8 +253,12 @@ table(merged_df2$Position%in%merged_df2v2$Position) rm(merged_df2v2) ######### -# merge 1b:remove duplicate mutation information +# merge 3b:remove duplicate mutation information ######### +cat('Merging dfs with NAs: small df (removing duplicate muts)' + ,'\nCannot trust lineage info from this' + ,'\nlinking col: Mutationinforamtion' + ,'\nfilename: merged_df3') #==#=#=#=#=#=# # Cannot trust lineage, country from this df as the same mutation @@ -280,69 +266,60 @@ rm(merged_df2v2) # 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 +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 if(nrow(mcsm_data) == nrow(merged_df3)){ - print("sanity check: Passed") + cat('PASS: No. of rows match with mcsm_data' + ,'\nExpected no. of rows: ', nrow(mcsm_data) + ,'\nGot no. of rows: ', nrow(merged_df3)) } else { - print("Error!: check data, nrows is not as expected") + cat('FAIL: No. of rows mismatch' + , '\nNo. of rows mcsm_data: ', nrow(mcsm_data) + , '\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,] -#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ########################### -# 3b :merging two dfs: without NA +# 4: merging two dfs: without NA ########################### +cat('Merging dfs without any NAs: big df (1-many relationship b/w id & mut)' + ,'\nlinking col: Mutationinforamtion' + ,'\nfilename: merged_df2_comp') ######### -# merge 2a:same as merge 1 but excluding NA +# merge 4a: same as merge 1 but excluding NA ######### -merged_df2_comp = merged_df2[!is.na(merged_df2$AF),] +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)) ######### -# merge 2b: remove duplicate mutation information +# merge 4b: remove duplicate mutation information ######### -merged_df3_comp = merged_df2_comp[!duplicated(merged_df2_comp$Mutationinformation),] +merged_df3_comp = merged_df2_comp[!duplicated(merged_df2_comp$Mutationinformation),] -# FIXME: add this as a sanity check. I have manually checked! +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)) #=============== end of combining df -#clear variables -rm(mcsm_data - , meta_with_afor - , foo) - -#rm(diff_n, my_merged, mcsm) - -#===============end of script - -#===================== +#********************* # write_output files -#===================== - -# Not required as this is a subset of the "combining_two_df.R" script +# Not required as this is a subset of the combining_two_df.R +#************************* +# 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) +rm(pos_count_check) +#============================= end of script diff --git a/meta_data_analysis/AF_and_OR_calcs.R b/meta_data_analysis/AF_and_OR_calcs.R index 051e777..f5c9db6 100644 --- a/meta_data_analysis/AF_and_OR_calcs.R +++ b/meta_data_analysis/AF_and_OR_calcs.R @@ -1,9 +1,9 @@ -#============================================ +######################################################### # TASK: To calculate Allele Frequency and # Odds Ratio from master data # and add the calculated params to meta_data extracted from -# pnca_data_extraction.py -#=========================================== +# data_extraction.py +######################################################### getwd() setwd('~/git/LSHTM_analysis/meta_data_analysis') getwd() @@ -14,9 +14,9 @@ gene = 'pncA' gene_match = paste0(gene,'_p.') cat(gene_match) -#======= -# input dir -#======= +#=========== +# input +#=========== # infile1: Raw data #indir = 'git/Data/pyrazinamide/input/original' indir = paste0('~/git/Data') @@ -24,27 +24,26 @@ in_filename = 'original_tanushree_data_v2.csv' infile = paste0(indir, '/', in_filename) cat(paste0('Reading infile1: raw data', ' ', infile) ) -# infile2: gene associated meta data file to extract valid snps and add calcs to -# filename: outfile3 from data_extraction.py +# infile2: gene associated meta data file to extract valid snps and add calcs to. +# This is outfile3 from data_extraction.py indir_metadata = paste0('~/git/Data', '/', drug, '/', 'output') in_filename_metadata = 'pnca_metadata.csv' infile_metadata = paste0(indir_metadata, '/', in_filename_metadata) cat(paste0('Reading infile2: gene associated metadata:', infile_metadata)) -#========= -# output dir -#========= +#=========== +# output +#=========== # outdir = 'git/Data/pyrazinamide/output' outdir = paste0('~/git/Data', '/', drug, '/', 'output') out_filename = paste0(tolower(gene),'_', 'meta_data_with_AFandOR.csv') outfile = paste0(outdir, '/', out_filename) cat(paste0('Output file with full path:', outfile)) - #%% end of variable assignment for input and output files -#=============== -# Step 1: Read master/raw data stored in Data/ -#=============== +######################################################### +# 1: Read master/raw data stored in Data/ +######################################################### raw_data_all = read.csv(infile, stringsAsFactors = F) raw_data = raw_data_all[,c("id" @@ -55,9 +54,9 @@ rm(raw_data_all) rm(indir, in_filename, infile) -##### +#=========== # 1a: exclude na -##### +#=========== raw_data = raw_data[!is.na(raw_data$pyrazinamide),] total_samples = length(unique(raw_data$id)) @@ -66,16 +65,16 @@ cat(paste0('Total samples without NA in', ' ', drug, 'is:', total_samples)) # sanity check: should be true is.numeric(total_samples) -##### +#=========== # 1b: combine the two mutation columns -##### +#=========== raw_data$all_mutations_pyrazinamide = paste(raw_data$dr_mutations_pyrazinamide , raw_data$other_mutations_pyrazinamide) head(raw_data$all_mutations_pyrazinamide) -##### +#=========== # 1c: create yet another column that contains all the mutations but in lower case -##### +#=========== raw_data$all_muts_pnca = tolower(raw_data$all_mutations_pyrazinamide) # sanity checks @@ -104,9 +103,10 @@ table(mut, dst) #fisher.test(table(mut, dst)) #table(mut) -#=============== -# Step 2: Read valid snps for which OR can be calculated (infile_comp_snps.csv) -#=============== +######################################################### +# 2: Read valid snps for which OR +# can be calculated (infile_comp_snps.csv) +######################################################### cat(paste0('Reading metadata infile:', infile_metadata)) pnca_metadata = read.csv(infile_metadata @@ -188,7 +188,7 @@ hist(log(ors) , breaks = 100 ) -# FIXME: could be good to add a sanity check +# sanity check if (table(names(ors) == names(pvals)) & table(names(ors) == names(afs)) & table(names(pvals) == names(afs)) == length(pnca_snps_unique)){ cat('PASS: names of ors, pvals and afs match: proceed with combining into a single df') } else{ @@ -214,9 +214,9 @@ if (table(rownames(comb_AF_and_OR) == comb_AF_and_OR$mutation)){ cat('FAIL: rownames and mutation col values mismatch') } -############ -# Merge 1: combine meta data file with calculated num params -########### +######################################################### +# 3: Merge meta data file + calculated num params +######################################################### df1 = pnca_metadata df2 = comb_AF_and_OR @@ -254,7 +254,7 @@ if ( identical(na_count[[length(na_count)]], na_count[[length(na_count)-1]], na_ cat('PASS: No. of NAs for OR, AF and Pvals are equal as expected', '\nNo. of NA: ', na_count[[length(na_count)]]) } else { - cat('FAIl: No. of NAs for OR, AF and Pvals mismatch') + cat('FAIL: No. of NAs for OR, AF and Pvals mismatch') } # reassign custom colnames @@ -289,6 +289,7 @@ cat(paste0('Added', ' ', ncol_added, ' more cols to merged_df:' , '\nno. of cols in merged_df now: ', ncol(merged_df))) #%% write file out: pnca_meta_data_with_AFandOR +#********************************************* cat(paste0('writing output file: ' , '\nFilename: ', out_filename , '\nPath:', outdir)) @@ -300,6 +301,98 @@ cat(paste0('Finished writing:' , out_filename , '\nNo. of rows: ', nrow(merged_df) , '\nNo. of cols: ', ncol(merged_df))) - +#************************************************ cat('======================================================================') rm(out_filename) +cat('End of script: calculated AF, OR, pvalues and saved file') +# End of script +#%% +# sanity check: Count NA in these four cols. +# However these need to be numeric else these +# will be misleading when counting NAs (i.e retrun 0) +#is.numeric(meta_with_afor$OR) +na_var = c('AF', 'OR', 'pvalue', 'logor', 'neglog10pvalue') + +# loop through these vars and check if these are numeric. +# if not, then convert to numeric +check_all = NULL + +for (i in na_var){ + # cat(i) + check0 = is.numeric(meta_with_afor[,i]) + if (check0) { + check_all = c(check0, check_all) + cat('These are all numeric cols') + } else{ + cat('First converting to numeric') + check0 = as.numeric(meta_with_afor[,i]) + check_all = c(check0, check_all) + } +} + +# count na now that the respective cols are numeric +na_count = sapply(meta_with_afor, function(y) sum(length(which(is.na(y))))); na_count +str(na_count) + +# extract how many NAs: +# should be all TRUE +# should be a single number since +# all the cols should have 'equal' and 'same' no. of NAs +# compare if the No of 'NA' are the same for all these cols +na_len = NULL +for (i in na_var){ + temp = na_count[[i]] + na_len = c(na_len, temp) +} + +cat('Checking how many NAs and if these are identical for the selected cols:') +my_nrows = NULL +for ( i in 1: (length(na_len)-1) ){ +# cat(compare(na_len[i]), na_len[i+1]) + c = compare(na_len[i], na_len[i+1]) + if ( c$result ) { + cat('PASS: No. of NAa in selected cols are identical') + my_nrows = na_len[i] } + else { + cat('FAIL: No. of NAa in selected cols mismatch') + } +} + +cat('No. of NAs in each of the selected cols: ', my_nrows) + +# yet more sanity checks: +cat('Check whether the ', my_nrows, 'indices are indeed the same') + +#which(is.na(meta_with_afor$OR)) + +# initialise an empty df with nrows as extracted above +na_count_df = data.frame(matrix(vector(mode = 'numeric' +# , length = length(na_var) + ) + , nrow = my_nrows +# , ncol = length(na_var) + )) + +# populate the df with the indices of the cols that are NA +for (i in na_var){ + cat(i) + na_i = which(is.na(meta_with_afor[i])) + na_count_df = cbind(na_count_df, na_i) + colnames(na_count_df)[which(na_var == i)] <- i +} + +# Now compare these indices to ensure these are the same +check2 = NULL +for ( i in 1: ( length(na_count_df)-1 ) ) { +# cat(na_count_df[i] == na_count_df[i+1]) + check_all = identical(na_count_df[[i]], na_count_df[[i+1]]) + check2 = c(check_all, check2) + if ( all(check2) ) { + cat('PASS: The indices for AF, OR, etc are all the same\n') + } else { + cat ('FAIL: Please check indices which are NA') + } +} + + +