getwd()
setwd("~/git/LSHTM_analysis/scripts/plotting/")
getwd()

#########################################################
# TASK: To combine struct params and meta data for plotting
# Input csv files:
# 1) <gene>_all_params.csv
# 2) <gene>_meta_data.csv

# Output: 
# 1) muts with opposite effects on stability
# 2) large combined df including NAs for AF, OR,etc
# 		Dim: same no. of rows as gene associated meta_data_with_AFandOR
# 3) small combined df including NAs for AF, OR, etc.
# 		Dim: same as mcsm data
# 4) large combined df excluding NAs 
# 		Dim: dim(#1) - no. of NAs(AF|OR) + 1
# 5) small combined df excluding NAs
# 		Dim: dim(#2) - no. of unique NAs - 1
# This script is sourced from other .R scripts for plotting
#########################################################

##########################################################
# 				Installing and loading required packages 			 
##########################################################
source("Header_TT.R")
#require(data.table)
#require(arsenal)
#require(compare)
#library(tidyverse)


#%% variable assignment: input and output paths & filenames
drug = "pyrazinamide"
gene = "pncA"
gene_match = paste0(gene,"_p.")
cat(gene_match)

#=============
# directories
#=============
datadir = paste0("~/git/Data")
indir = paste0(datadir, "/", drug, "/input")
outdir = paste0("~/git/Data", "/", drug, "/output")

#===========
# input
#===========
#in_filename = "mcsm_complex1_normalised.csv"
in_filename_params = paste0(tolower(gene), "_all_params.csv") 
infile_params = paste0(outdir, "/", in_filename_params)
cat(paste0("Input file 1:", infile_params) )

# infile 2: gene associated meta data
#in_filename_gene_metadata = paste0(tolower(gene),  "_meta_data_with_AFandOR.csv")
in_filename_gene_metadata = paste0(tolower(gene),  "_metadata.csv")
infile_gene_metadata = paste0(outdir, "/", in_filename_gene_metadata)
cat(paste0("Input infile 2:", infile_gene_metadata))

#===========
# output
#===========
# mutations with opposite effects
out_filename_opp_muts = paste0(tolower(gene), "_muts_opp_effects.csv")
outfile_opp_muts = paste0(outdir, "/", out_filename_opp_muts)


#%%===============================================================
###########################
# Read file: struct params
###########################
cat("Reading struct params including mcsm:"
    , in_filename_params)

mcsm_data = read.csv(infile_params
                     #, row.names = 1
                     , stringsAsFactors = F
                     , header = T) 

cat("Input dimensions:", dim(mcsm_data)) #416, 86

# clear variables
rm(in_filename_params, infile_params)

str(mcsm_data)

table(mcsm_data$duet_outcome); sum(table(mcsm_data$duet_outcome) )

# spelling Correction 1: DUET incase American spelling needed!
#mcsm_data$duet_outcome[mcsm_data$duet_outcome=="Stabilising"] <- "Stabilizing"
#mcsm_data$duet_outcome[mcsm_data$duet_outcome=="Destabilising"] <- "Destabilizing"

# 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 incase American spelling needed!
table(mcsm_data$ligand_outcome); sum(table(mcsm_data$ligand_outcome) )
#mcsm_data$ligand_outcome[mcsm_data$ligand_outcome=="Stabilising"] <- "Stabilizing"
#mcsm_data$ligand_outcome[mcsm_data$ligand_outcome=="Destabilising"] <- "Destabilizing"

# checks: should be the same as above
table(mcsm_data$ligand_outcome); sum(table(mcsm_data$ligand_outcome) )
head(mcsm_data$ligand_outcome); tail(mcsm_data$ligand_outcome)

# muts with opposing effects on protomer and ligand stability
table(mcsm_data$duet_outcome != mcsm_data$ligand_outcome)
changes = mcsm_data[which(mcsm_data$duet_outcome != mcsm_data$ligand_outcome),]

# sanity check: redundant, but uber cautious!
dl_i = which(mcsm_data$duet_outcome != mcsm_data$ligand_outcome)
ld_i = which(mcsm_data$ligand_outcome != mcsm_data$duet_outcome)

cat("Identifying muts with opposite stability effects")
if(nrow(changes) == (table(mcsm_data$duet_outcome != mcsm_data$ligand_outcome)[[2]]) & identical(dl_i,ld_i)) {
  cat("PASS: muts with opposite effects on stability and affinity correctly identified"
        , "\nNo. of such muts: ", nrow(changes))
}else {
  cat("FAIL: unsuccessful in extracting muts with changed stability effects")
}

#***************************
# write file: changed muts
write.csv(changes, outfile_opp_muts)

cat("Finished writing file for muts with opp effects:"
    , "\nFilename: ", outfile_opp_muts
    , "\nDim:", dim(changes))

# clear variables
rm(out_filename_opp_muts, outfile_opp_muts)
rm(changes, dl_i, ld_i)

#***************************
# count na in each column
na_count = sapply(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)

df_ncols = ncol(mcsm_data)

# REMOVE as this is dangerous due to dup muts
# 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: ", df_ncols)

#pos_count_check = data.frame(mcsm_data$position, mcsm_data$occurrence)

# check duplicate muts
if (length(unique(mcsm_data$mutationinformation)) == length(mcsm_data$mutationinformation)){
  cat("No duplicate mutations in mcsm data")
}else{
  dup_muts = mcsm_data[duplicated(mcsm_data$mutationinformation),]
  dup_muts_nu = length(unique(dup_muts$mutationinformation))
  cat(paste0("CAUTION:", nrow(dup_muts), " Duplicate mutations identified"
             , "\nOf these, no. of unique mutations are:", dup_muts_nu
             , "\nExtracting df with unique mutations only"))
  mcsm_data_u = mcsm_data[!duplicated(mcsm_data$mutationinformation),]
}

if (nrow(mcsm_data_u) == length(unique(mcsm_data$mutationinformation))){
  cat("Df without duplicate mutations successfully extracted")
} else{
  cat("FAIL: could not extract clean df!")
  quit()
}

###########################
# 2: Read file: <gene>_meta data.csv
###########################
cat("Reading meta data file:", infile_gene_metadata)

gene_metadata <- read.csv(infile_gene_metadata
                      , stringsAsFactors = F
                      , header = T)
cat("Dim:", dim(gene_metadata))

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# FIXME: remove
# counting NAs in AF, OR cols: 
if (identical(sum(is.na(gene_metadata$OR))
              , sum(is.na(gene_metadata$pvalue))
              , sum(is.na(gene_metadata$AF)))){
  cat("PASS: NA count match for OR, pvalue and AF\n")
  na_count = sum(is.na(gene_metadata$AF))
  cat("No. of NAs: ", sum(is.na(gene_metadata$OR)))
} else{
  cat("FAIL: NA count mismatch"
      , "\nNA in OR: ", sum(is.na(gene_metadata$OR))
      , "\nNA in pvalue: ", sum(is.na(gene_metadata$pvalue))
      , "\nNA in AF:", sum(is.na(gene_metadata$AF)))
}
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# clear variables
rm(in_filename_gene_metadata, infile_gene_metadata)

str(gene_metadata)

# sort by position (same as mcsm_data)
# earlier it was mutationinformation
#head(gene_metadata$mutationinformation)
#gene_metadata = gene_metadata[order(gene_metadata$mutationinformation),]
##head(gene_metadata$mutationinformation)

head(gene_metadata$position)
gene_metadata = gene_metadata[order(gene_metadata$position),]
head(gene_metadata$position)

###########################
# Merge 1: two dfs with NA
# merged_df2
###########################
head(mcsm_data$mutationinformation)
head(gene_metadata$mutationinformation)

# Find common columns b/w two df
merging_cols = intersect(colnames(mcsm_data), colnames(gene_metadata))

cat(paste0("Merging dfs with NAs: big df (1-many relationship b/w id & mut)"
           , "\nNo. of merging cols:", length(merging_cols)
           , "\nMerging columns identified:"))
print(merging_cols)

#=============
# merged_df2): gene_metadata + mcsm_data
#==============
merged_df2 = merge(x = gene_metadata
                  , y = mcsm_data
                  , by = merging_cols
                  , all.y = T)

cat("Dim of merged_df2: ", dim(merged_df2) #4520, 11
    )
head(merged_df2$position)

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# FIXME: count how many unique muts have entries in meta data
# sanity check
cat("Checking nrows in merged_df2")
if(nrow(gene_metadata) == nrow(merged_df2)){
  cat("nrow(merged_df2) = nrow (gene associated gene_metadata)"
      ,"\nExpected no. of rows: ",nrow(gene_metadata) 
      ,"\nGot no. of rows: ", nrow(merged_df2))
} else{
  cat("nrow(merged_df2)!= nrow(gene associated gene_metadata)"
      , "\nExpected no. of rows after merge: ", nrow(gene_metadata)
      , "\nGot no. of rows: ", nrow(merged_df2)
      , "\nFinding discrepancy")
  merged_muts_u = unique(merged_df2$mutationinformation)
  meta_muts_u = unique(gene_metadata$mutationinformation)
    # find the index where it differs
  unique(meta_muts_u[! meta_muts_u %in% merged_muts_u])
}

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

# sort by position
head(merged_df2$position)
merged_df2 = merged_df2[order(merged_df2$position),]
head(merged_df2$position)

merged_df2v3 = merge(x = gene_metadata
                   , y = mcsm_data
                   , by = merging_cols
                   , all = T) 

merged_df2v2 = merge(x = gene_metadata
                     , y = mcsm_data
                     , by = merging_cols
                     , 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 mutation
#!=!=!=!=!=!=!=!

# should  be False
identical(merged_df2, merged_df2v2)
table(merged_df2$position%in%merged_df2v2$position)

rm(merged_df2v2)

#!!!!!!!!!!! check why these differ

#########
# merge 3b (merged_df3):remove duplicate mutation information
#########
cat("Merging dfs without NAs: small df (removing muts with no AF|OR associated)"
    ,"\nCannot trust lineage info from this"
    ,"\nlinking col: Mutationinforamtion"
    ,"\nfilename: merged_df3")

#==#=#=#=#=#=#
# Cannot trust lineage, country from this df as the same mutation
# can have many different lineages
# but this should be good for the numerical corr plots
#=#=#=#=#=#=#=
merged_df3 = merged_df2[!duplicated(merged_df2$mutationinformation),] 
head(merged_df3$position); tail(merged_df3$position) # should be sorted

# sanity check
cat("Checking nrows in merged_df3")
if(nrow(mcsm_data) == nrow(merged_df3)){
  cat("PASS: No. of rows match with mcsm_data"
      ,"\nExpected no. of rows: ", nrow(mcsm_data)
      ,"\nGot no. of rows: ", nrow(merged_df3))
} else {
  cat("FAIL: No. of rows mismatch"
      , "\nNo. of rows mcsm_data: ", nrow(mcsm_data)
      , "\nNo. of rows merged_df3: ", nrow(merged_df3))
}

# 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")

merged_df2_comp = merged_df2[!is.na(merged_df2$AF),] 
#merged_df2_comp = merged_df2[!duplicated(merged_df2$mutationinformation),] 

# 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 (merged_df3_comp): 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
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
#*********************
# 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 file: "
#    ,"\nFilename: ", out_filename
#    ,"\nPath: ", outdir)

#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)))

#%% 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, gene_metadata, foo, drug, gene, gene_match, indir, merged_muts_u, meta_muts_u, na_count, df_ncols, outdir)
rm(pos_count_check)
#============================= end of script