299 lines
8.2 KiB
R
299 lines
8.2 KiB
R
getwd()
|
|
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 #
|
|
########################################################################
|
|
|
|
source("Header_TT.R")
|
|
#require(data.table)
|
|
#require(arsenal)
|
|
#require(compare)
|
|
#library(tidyverse)
|
|
|
|
#################################
|
|
# Read file: normalised file
|
|
# output of step 4 mcsm_pipeline
|
|
#################################
|
|
|
|
inDir = "~/git/Data/pyrazinamide/input/processed/"
|
|
inFile = paste0(inDir, "mcsm_complex1_normalised.csv"); inFile
|
|
|
|
mcsm_data = read.csv(inFile
|
|
, row.names = 1
|
|
, stringsAsFactors = F
|
|
, header = T)
|
|
rm(inDir, inFile)
|
|
|
|
str(mcsm_data)
|
|
|
|
table(mcsm_data$DUET_outcome); sum(table(mcsm_data$DUET_outcome) )
|
|
|
|
# spelling Correction 1: DUET
|
|
mcsm_data$DUET_outcome[mcsm_data$DUET_outcome=='Stabilizing'] <- 'Stabilising'
|
|
mcsm_data$DUET_outcome[mcsm_data$DUET_outcome=='Destabilizing'] <- 'Destabilising'
|
|
|
|
# checks: should be the same as above
|
|
table(mcsm_data$DUET_outcome); sum(table(mcsm_data$DUET_outcome) )
|
|
head(mcsm_data$DUET_outcome); tail(mcsm_data$DUET_outcome)
|
|
|
|
# spelling Correction 2: Ligand
|
|
table(mcsm_data$Lig_outcome); sum(table(mcsm_data$Lig_outcome) )
|
|
|
|
mcsm_data$Lig_outcome[mcsm_data$Lig_outcome=='Stabilizing'] <- 'Stabilising'
|
|
mcsm_data$Lig_outcome[mcsm_data$Lig_outcome=='Destabilizing'] <- 'Destabilising'
|
|
|
|
# checks: should be the same as above
|
|
table(mcsm_data$Lig_outcome); sum(table(mcsm_data$Lig_outcome) )
|
|
head(mcsm_data$Lig_outcome); tail(mcsm_data$Lig_outcome)
|
|
|
|
# count na in each column
|
|
na_count = sapply(mcsm_data, function(y) sum(length(which(is.na(y))))); na_count
|
|
|
|
# sort by Mutationinformation
|
|
mcsm_data = mcsm_data[order(mcsm_data$Mutationinformation),]
|
|
head(mcsm_data$Mutationinformation)
|
|
|
|
# get freq count of positions and add to the df
|
|
setDT(mcsm_data)[, occurrence := .N, by = .(Position)]
|
|
|
|
pos_count_check = data.frame(mcsm_data$Position, mcsm_data$occurrence)
|
|
|
|
###########################
|
|
# 2: Read file: meta data with AFandOR
|
|
###########################
|
|
|
|
inDir = "~/git/Data/pyrazinamide/input/processed/"
|
|
inFile2 = paste0(inDir, "meta_data_with_AFandOR.csv"); inFile2
|
|
|
|
meta_with_afor <- read.csv(inFile2
|
|
, stringsAsFactors = F
|
|
, header = T)
|
|
|
|
rm(inDir, inFile2)
|
|
|
|
str(meta_with_afor)
|
|
|
|
# sort by Mutationinformation
|
|
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
|
|
###########################
|
|
|
|
# link col name = Mutationinforamtion
|
|
head(mcsm_data$Mutationinformation)
|
|
head(meta_with_afor$Mutationinformation)
|
|
|
|
#########
|
|
# merge 1a: meta data with mcsm
|
|
#########
|
|
merged_df2 = merge(x = meta_with_afor
|
|
,y = mcsm_data
|
|
, by = "Mutationinformation"
|
|
, all.y = T)
|
|
|
|
head(merged_df2$Position)
|
|
|
|
# sort by Position
|
|
head(merged_df2$Position)
|
|
merged_df2 = merged_df2[order(merged_df2$Position),]
|
|
head(merged_df2$Position)
|
|
|
|
merged_df2v2 = merge(x = meta_with_afor
|
|
,y = mcsm_data
|
|
, 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
|
|
#!=!=!=!=!=!=!=!
|
|
|
|
# should be False
|
|
identical(merged_df2, merged_df2v2)
|
|
table(merged_df2$Position%in%merged_df2v2$Position)
|
|
|
|
rm(merged_df2v2)
|
|
|
|
#########
|
|
# merge 1b:remove duplicate mutation information
|
|
#########
|
|
|
|
#==#=#=#=#=#=#
|
|
# 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 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")
|
|
} else {
|
|
print("Error!: check data, nrows is not as expected")
|
|
}
|
|
|
|
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
# 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
|
|
###########################
|
|
|
|
#########
|
|
# merge 2a:same as merge 1 but excluding NA
|
|
#########
|
|
merged_df2_comp = merged_df2[!is.na(merged_df2$AF),]
|
|
|
|
#########
|
|
# merge 2b: remove duplicate mutation information
|
|
#########
|
|
merged_df3_comp = merged_df2_comp[!duplicated(merged_df2_comp$Mutationinformation),]
|
|
|
|
# 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)
|
|
|
|
#=====================
|
|
# write_output files
|
|
#=====================
|
|
# output dir
|
|
outDir = "~/git/Data/pyrazinamide/output/"
|
|
getwd()
|
|
|
|
outFile1 = paste0(outDir, "merged_df3.csv"); outFile1
|
|
write.csv(merged_df3, outFile1)
|
|
|
|
#outFile2 = paste0(outDir, "merged_df3_comp.csv"); outFile2
|
|
#write.csv(merged_df3_comp, outFile2)
|
|
|
|
rm(outDir
|
|
, outFile1
|
|
# , outFile2
|
|
)
|
|
#============================= end of script
|
|
|