tanu/LSHTM_analysis
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

updated combining df scripts for duet & lig

Browse source
This commit is contained in:
Tanushree Tunstall 2020-03-24 18:28:52 +00:00
parent c184841951
commit 96ebb85069
3 changed files with 458 additions and 342 deletions
Download patch file Download diff file Expand all files Collapse all files

320
mcsm_analysis/pyrazinamide/scripts/combining_two_df.R
View file

@ -1,17 +1,17 @@
getwd()
setwd("~/git/LSHTM_analysis/mcsm_analysis/pyrazinamide/scripts/")
getwd()
######################################################### #########################################################
# TASK: To combine mcsm and meta data with af and or # 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/')
getwd()
##########################################################
# Installing and loading required packages
##########################################################
######################################################################## source('Header_TT.R')
# Installing and loading required packages #
########################################################################
#source("Header_TT.R")
#require(data.table) #require(data.table)
#require(arsenal) #require(arsenal)
#require(compare) #require(compare)
@ -22,14 +22,58 @@ getwd()
# output of step 4 mcsm_pipeline # output of step 4 mcsm_pipeline
################################# #################################
inDir = "~/git/Data/pyrazinamide/input/processed/" #%% variable assignment: input and output paths & filenames
inFile = paste0(inDir, "mcsm_complex1_normalised.csv"); inFile 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 , row.names = 1
, stringsAsFactors = F , stringsAsFactors = F
, header = T) , header = T)
rm(inDir, inFile)
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) 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) ) table(mcsm_data$Lig_outcome); sum(table(mcsm_data$Lig_outcome) )
head(mcsm_data$Lig_outcome); tail(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 # count na in each column
na_count = sapply(mcsm_data, function(y) sum(length(which(is.na(y))))); na_count 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),] mcsm_data = mcsm_data[order(mcsm_data$Mutationinformation),]
head(mcsm_data$Mutationinformation) head(mcsm_data$Mutationinformation)
orig_col = ncol(mcsm_data)
# get freq count of positions and add to the df # 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) pos_count_check = data.frame(mcsm_data$Position, mcsm_data$occurrence)
########################### ###########################
# 2: Read file: meta data with AFandOR # 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/" meta_with_afor <- read.csv(infile_comb
inFile2 = paste0(inDir, "meta_data_with_AFandOR.csv"); inFile2
meta_with_afor <- read.csv(inFile2
, stringsAsFactors = F , stringsAsFactors = F
, header = T) , 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) str(meta_with_afor)
@ -85,113 +168,45 @@ head(meta_with_afor$Mutationinformation)
meta_with_afor = meta_with_afor[order(meta_with_afor$Mutationinformation),] meta_with_afor = meta_with_afor[order(meta_with_afor$Mutationinformation),]
head(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(mcsm_data$Mutationinformation)
head(meta_with_afor$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 merged_df2 = merge(x = meta_with_afor
,y = mcsm_data ,y = mcsm_data
, by = "Mutationinformation" , by = 'Mutationinformation'
, all.y = T) , all.y = T)
cat('Dim of merged_df2: '
, '\nNo. of rows: ', nrow(merged_df2)
, '\nNo. of cols: ', ncol(merged_df2))
head(merged_df2$Position) 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 # sort by Position
head(merged_df2$Position) head(merged_df2$Position)
merged_df2 = merged_df2[order(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 merged_df2v2 = merge(x = meta_with_afor
,y = mcsm_data ,y = mcsm_data
, by = "Mutationinformation" , by = 'Mutationinformation'
, all.x = T) , all.x = T)
#!=!=!=!=!=!=!=! #!=!=!=!=!=!=!=!
# COMMENT: used all.y since position 186 is not part of the struc, # COMMENT: used all.y since position 186 is not part of the struc,
# hence doesn't have a mcsm value # hence doesn't have a mcsm value
# but 186 is associated with with mutation # but 186 is associated with mutation
#!=!=!=!=!=!=!=! #!=!=!=!=!=!=!=!
# should be False # should be False
@ -214,8 +229,12 @@ table(merged_df2$Position%in%merged_df2v2$Position)
rm(merged_df2v2) 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 # 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 # sanity checks
# nrows of merged_df3 should be the same as the nrows of mcsm_data # nrows of merged_df3 should be the same as the nrows of mcsm_data
if(nrow(mcsm_data) == nrow(merged_df3)){ 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 { } 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[!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),]
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 # alternate way of deriving merged_df3_comp
foo = merged_df3[!is.na(merged_df3$AF),] foo = merged_df3[!is.na(merged_df3$AF),]
# compare dfs: foo and merged_df3_com # compare dfs: foo and merged_df3_com
@ -271,29 +306,40 @@ all.equal(foo, merged_df3)
summary(comparedf(foo, merged_df3)) summary(comparedf(foo, merged_df3))
#=============== end of combining df #=============== end of combining df
#clear variables #*********************
rm(mcsm_data
, meta_with_afor
, foo)
#rm(diff_n, my_merged, mcsm)
#=====================
# write_output files # write_output files
#=====================
# output dir # output dir
outDir = "~/git/Data/pyrazinamide/output/" #outdir = '~/git/Data/pyrazinamide/output/'
getwd() #uncomment as necessary
#FIXME
#out_filenames = c('merged_df2'
# , 'merged_df3'
# , 'meregd_df2_comp'
# , 'merged_df3_comp'
#)
outFile1 = paste0(outDir, "merged_df3.csv"); outFile1 #cat('Writing output files: '
#write.csv(merged_df3, outFile1) # , '\nPath:', outdir)
#outFile2 = paste0(outDir, "merged_df3_comp.csv"); outFile2 #for (i in out_filenames){
#write.csv(merged_df3_comp, outFile2) # 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 #sapply(out_filenames, function(x) write.csv(x, 'x.csv'))
, outFile1 #*************************
# , outFile2 # 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 #============================= end of script

305
mcsm_analysis/pyrazinamide/scripts/combining_two_df_lig.R
View file

@ -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 # 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(data.table)
#require(arsenal) #require(arsenal)
#require(compare) #require(compare)
@ -22,14 +23,58 @@ getwd()
# output of step 4 mcsm_pipeline # output of step 4 mcsm_pipeline
################################# #################################
inDir = "~/git/Data/pyrazinamide/input/processed/" #%% variable assignment: input and output paths & filenames
inFile = paste0(inDir, "mcsm_complex1_normalised.csv"); inFile 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 , row.names = 1
, stringsAsFactors = F , stringsAsFactors = F
, header = T) , header = T)
rm(inDir, inFile)
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) 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=='Stabilizing'] <- 'Stabilising'
mcsm_data$DUET_outcome[mcsm_data$DUET_outcome=='Destabilizing'] <- 'Destabilising' 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) ) table(mcsm_data$DUET_outcome); sum(table(mcsm_data$DUET_outcome) )
head(mcsm_data$DUET_outcome); tail(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) ) table(mcsm_data$Lig_outcome); sum(table(mcsm_data$Lig_outcome) )
head(mcsm_data$Lig_outcome); tail(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 ########################### !!! only for mcsm_lig
# 4: Filter/subset data # 4: Filter/subset data
# Lig plots < 10Ang # Lig plots < 10Ang
@ -82,57 +131,60 @@ length(unique(mcsm_data2$Mutationinformation))
# count Destabilisinga and stabilising # count Destabilisinga and stabilising
table(mcsm_data2$Lig_outcome) #{RESULT: no of mutations within 10Ang} 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: # Extra sanity check:
# for mcsm_lig ONLY # for mcsm_lig ONLY
# Dis_lig_Ang should be <10 # 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") print ("Sanity check passed: lig data is <10Ang")
}else{ }else{
print ("Error: data should be filtered to be within 10Ang") print ("Error: data should be filtered to be within 10Ang")
} }
#!!!!!!!!!!!!!!!!!!!!!
# REASSIGNMENT: so as not to alter the script
mcsm_data = mcsm_data2
#!!!!!!!!!!!!!!!!!!!!!
# clear variables # clear variables
rm(mcsm_data2) rm(mcsm_data2)
# count na in each column # count na in each column
na_count = sapply(mcsm_data, function(y) sum(length(which(is.na(y))))); na_count 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 # sort by Mutationinformation
mcsm_data = mcsm_data[order(mcsm_data$Mutationinformation),] mcsm_data = mcsm_data[order(mcsm_data$Mutationinformation),]
head(mcsm_data$Mutationinformation) head(mcsm_data$Mutationinformation)
# check orig_col = ncol(mcsm_data)
mcsm_data[grep("Q10P", mcsm_data$Mutationinformation),]
mcsm_data[grep("A102T", mcsm_data$Mutationinformation),]
# get freq count of positions and add to the df # 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) pos_count_check = data.frame(mcsm_data$Position, mcsm_data$occurrence)
########################### ###########################
# 2: Read file: meta data with AFandOR # 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/" meta_with_afor <- read.csv(infile_comb
inFile2 = paste0(inDir, "meta_data_with_AFandOR.csv"); inFile2
meta_with_afor <- read.csv(inFile2
, stringsAsFactors = F , stringsAsFactors = F
, header = T) , 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) str(meta_with_afor)
# sort by Mutationinformation # sort by Mutationinformation
@ -140,114 +192,45 @@ head(meta_with_afor$Mutationinformation)
meta_with_afor = meta_with_afor[order(meta_with_afor$Mutationinformation),] meta_with_afor = meta_with_afor[order(meta_with_afor$Mutationinformation),]
head(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(mcsm_data$Mutationinformation)
head(meta_with_afor$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 merged_df2 = merge(x = meta_with_afor
,y = mcsm_data ,y = mcsm_data
, by = "Mutationinformation" , by = 'Mutationinformation'
, all.y = T) , all.y = T)
cat('Dim of merged_df2: '
, '\nNo. of rows: ', nrow(merged_df2)
, '\nNo. of cols: ', ncol(merged_df2))
head(merged_df2$Position) 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 # sort by Position
head(merged_df2$Position) head(merged_df2$Position)
merged_df2 = merged_df2[order(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 merged_df2v2 = merge(x = meta_with_afor
,y = mcsm_data ,y = mcsm_data
, by = "Mutationinformation" , by = 'Mutationinformation'
, all.x = T) , all.x = T)
#!=!=!=!=!=!=!=! #!=!=!=!=!=!=!=!
# COMMENT: used all.y since position 186 is not part of the struc, # COMMENT: used all.y since position 186 is not part of the struc,
# hence doesn't have a mcsm value # hence doesn't have a mcsm value
# but 186 is associated with with mutation # but 186 is associated with mutation
#!=!=!=!=!=!=!=! #!=!=!=!=!=!=!=!
# should be False # should be False
@ -271,8 +253,12 @@ table(merged_df2$Position%in%merged_df2v2$Position)
rm(merged_df2v2) 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 # Cannot trust lineage, country from this df as the same mutation
@ -285,64 +271,55 @@ head(merged_df3$Position) ; tail(merged_df3$Position) # should be sorted
# sanity checks # sanity checks
# nrows of merged_df3 should be the same as the nrows of mcsm_data # nrows of merged_df3 should be the same as the nrows of mcsm_data
if(nrow(mcsm_data) == nrow(merged_df3)){ 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 { } 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 # alternate way of deriving merged_df3_comp
foo = merged_df3[!is.na(merged_df3$AF),] foo = merged_df3[!is.na(merged_df3$AF),]
# compare dfs: foo and merged_df3_com # compare dfs: foo and merged_df3_com
all.equal(foo, merged_df3) all.equal(foo, merged_df3)
summary(comparedf(foo, merged_df3)) summary(comparedf(foo, merged_df3))
#=============== end of combining df #=============== 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 # write_output files
#===================== # Not required as this is a subset of the combining_two_df.R
#*************************
# Not required as this is a subset of the "combining_two_df.R" script # 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

153
meta_data_analysis/AF_and_OR_calcs.R
View file

@ -1,9 +1,9 @@
#============================================ #########################################################
# TASK: To calculate Allele Frequency and # TASK: To calculate Allele Frequency and
# Odds Ratio from master data # Odds Ratio from master data
# and add the calculated params to meta_data extracted from # and add the calculated params to meta_data extracted from
# pnca_data_extraction.py # data_extraction.py
#=========================================== #########################################################
getwd() getwd()
setwd('~/git/LSHTM_analysis/meta_data_analysis') setwd('~/git/LSHTM_analysis/meta_data_analysis')
getwd() getwd()
@ -14,9 +14,9 @@ gene = 'pncA'
gene_match = paste0(gene,'_p.') gene_match = paste0(gene,'_p.')
cat(gene_match) cat(gene_match)
#======= #===========
# input dir # input
#======= #===========
# infile1: Raw data # infile1: Raw data
#indir = 'git/Data/pyrazinamide/input/original' #indir = 'git/Data/pyrazinamide/input/original'
indir = paste0('~/git/Data') indir = paste0('~/git/Data')
@ -24,27 +24,26 @@ in_filename = 'original_tanushree_data_v2.csv'
infile = paste0(indir, '/', in_filename) infile = paste0(indir, '/', in_filename)
cat(paste0('Reading infile1: raw data', ' ', infile) ) cat(paste0('Reading infile1: raw data', ' ', infile) )
# infile2: gene associated meta data file to extract valid snps and add calcs to # infile2: gene associated meta data file to extract valid snps and add calcs to.
# filename: outfile3 from data_extraction.py # This is outfile3 from data_extraction.py
indir_metadata = paste0('~/git/Data', '/', drug, '/', 'output') indir_metadata = paste0('~/git/Data', '/', drug, '/', 'output')
in_filename_metadata = 'pnca_metadata.csv' in_filename_metadata = 'pnca_metadata.csv'
infile_metadata = paste0(indir_metadata, '/', in_filename_metadata) infile_metadata = paste0(indir_metadata, '/', in_filename_metadata)
cat(paste0('Reading infile2: gene associated metadata:', infile_metadata)) cat(paste0('Reading infile2: gene associated metadata:', infile_metadata))
#========= #===========
# output dir # output
#========= #===========
# outdir = 'git/Data/pyrazinamide/output' # outdir = 'git/Data/pyrazinamide/output'
outdir = paste0('~/git/Data', '/', drug, '/', 'output') outdir = paste0('~/git/Data', '/', drug, '/', 'output')
out_filename = paste0(tolower(gene),'_', 'meta_data_with_AFandOR.csv') out_filename = paste0(tolower(gene),'_', 'meta_data_with_AFandOR.csv')
outfile = paste0(outdir, '/', out_filename) outfile = paste0(outdir, '/', out_filename)
cat(paste0('Output file with full path:', outfile)) cat(paste0('Output file with full path:', outfile))
#%% end of variable assignment for input and output files #%% 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_all = read.csv(infile, stringsAsFactors = F)
raw_data = raw_data_all[,c("id" raw_data = raw_data_all[,c("id"
@ -55,9 +54,9 @@ rm(raw_data_all)
rm(indir, in_filename, infile) rm(indir, in_filename, infile)
##### #===========
# 1a: exclude na # 1a: exclude na
##### #===========
raw_data = raw_data[!is.na(raw_data$pyrazinamide),] raw_data = raw_data[!is.na(raw_data$pyrazinamide),]
total_samples = length(unique(raw_data$id)) 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 # sanity check: should be true
is.numeric(total_samples) is.numeric(total_samples)
##### #===========
# 1b: combine the two mutation columns # 1b: combine the two mutation columns
##### #===========
raw_data$all_mutations_pyrazinamide = paste(raw_data$dr_mutations_pyrazinamide raw_data$all_mutations_pyrazinamide = paste(raw_data$dr_mutations_pyrazinamide
, raw_data$other_mutations_pyrazinamide) , raw_data$other_mutations_pyrazinamide)
head(raw_data$all_mutations_pyrazinamide) head(raw_data$all_mutations_pyrazinamide)
##### #===========
# 1c: create yet another column that contains all the mutations but in lower case # 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) raw_data$all_muts_pnca = tolower(raw_data$all_mutations_pyrazinamide)
# sanity checks # sanity checks
@ -104,9 +103,10 @@ table(mut, dst)
#fisher.test(table(mut, dst)) #fisher.test(table(mut, dst))
#table(mut) #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)) cat(paste0('Reading metadata infile:', infile_metadata))
pnca_metadata = read.csv(infile_metadata pnca_metadata = read.csv(infile_metadata
@ -188,7 +188,7 @@ hist(log(ors)
, breaks = 100 , 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)){ 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') cat('PASS: names of ors, pvals and afs match: proceed with combining into a single df')
} else{ } 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') 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 df1 = pnca_metadata
df2 = comb_AF_and_OR 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', cat('PASS: No. of NAs for OR, AF and Pvals are equal as expected',
'\nNo. of NA: ', na_count[[length(na_count)]]) '\nNo. of NA: ', na_count[[length(na_count)]])
} else { } 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 # 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))) , '\nno. of cols in merged_df now: ', ncol(merged_df)))
#%% write file out: pnca_meta_data_with_AFandOR #%% write file out: pnca_meta_data_with_AFandOR
#*********************************************
cat(paste0('writing output file: ' cat(paste0('writing output file: '
, '\nFilename: ', out_filename , '\nFilename: ', out_filename
, '\nPath:', outdir)) , '\nPath:', outdir))
@ -300,6 +301,98 @@ cat(paste0('Finished writing:'
, out_filename , out_filename
, '\nNo. of rows: ', nrow(merged_df) , '\nNo. of rows: ', nrow(merged_df)
, '\nNo. of cols: ', ncol(merged_df))) , '\nNo. of cols: ', ncol(merged_df)))
#************************************************
cat('======================================================================') cat('======================================================================')
rm(out_filename) 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')
}
}