#!/usr/bin/env Rscript             
#########################################################
# TASK: formatting data that will be used for various plots

# useful links
#https://stackoverflow.com/questions/38851592/r-append-column-in-a-dataframe-with-frequency-count-based-on-two-columns
#########################################################
# working dir and loading libraries
getwd()
setwd("~/git/LSHTM_analysis/scripts/plotting")
getwd()

#source("Header_TT.R")
library(ggplot2)
library(data.table)
library(dplyr)
require("getopt", quietly = TRUE) #cmd parse arguments
source("dirs.R")

#========================================================
# command line args
spec = matrix(c(
  "drug"   , "d", 1, "character",
  "gene"   , "g", 1, "character"
), byrow = TRUE, ncol = 4)

opt = getopt(spec)

#FIXME: detect if script running from cmd, then set these
#drug = opt$drug
#gene = opt$gene

# hardcoding when not using cmd
drug = "streptomycin"
gene = "gid"

if(is.null(drug)|is.null(gene)) {
  stop("Missing arguments: --drug and --gene must both be specified (case-sensitive)")
}

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


cat('columns based on variables:\n'
      , drug
      , '\n'
      , dr_muts_col
      , '\n'
      , other_muts_col
      , "\n"
      , resistance_col
      , '\n===============================================================')

#%%===============================================================
###########################
# Read file: struct params
###########################
cat("Reading struct params including mcsm:", in_filename_params)
    
my_df = read.csv(infile_params, header = T)

cat("\nInput dimensions:", dim(my_df)) 

###########################
# add foldx outcome category
# and foldx scaled values 

# This will enable to always have these variables available
# when calling for plots
###########################

#------------------------------
# adding foldx scaled values
# scale data b/w -1 and 1
#------------------------------
n = which(colnames(my_df) == "ddg"); n 

my_min = min(my_df[,n]); my_min 
my_max = max(my_df[,n]); my_max 

my_df$foldx_scaled = ifelse(my_df[,n] < 0
                            , my_df[,n]/abs(my_min)
                            , my_df[,n]/my_max) 
# sanity check
my_min = min(my_df$foldx_scaled); my_min 
my_max = max(my_df$foldx_scaled); my_max

if (my_min == -1 && my_max == 1){
  cat("\nPASS: foldx ddg successfully scaled b/w -1 and 1"
      , "\nProceeding with assigning foldx outcome category")
}else{
  cat("\nFAIL: could not scale foldx ddg values"
      , "Aborting!")
}

#------------------------------
# adding foldx outcome category
# ddg<0 = "Stabilising" (-ve)
#------------------------------
c1 = table(my_df$ddg < 0)
my_df$foldx_outcome = ifelse(my_df$ddg < 0, "Stabilising", "Destabilising")
c2 = table(my_df$ddg < 0)

if ( all(c1 == c2) ){
  cat("\nPASS: foldx outcome successfully created")
}else{
  cat("\nFAIL: foldx outcome could not be created. Aborting!")
  exit()
}

###########################
# extract unique mutation entries
###########################

# check for duplicate mutations
if ( length(unique(my_df$mutationinformation)) != length(my_df$mutationinformation)){
  cat(paste0("\nCAUTION:", " Duplicate mutations identified"
             , "\nExtracting these..."))
  dup_muts = my_df[duplicated(my_df$mutationinformation),]
  dup_muts_nu = length(unique(dup_muts$mutationinformation))
  cat(paste0("\nDim of duplicate mutation df:", nrow(dup_muts)
             , "\nNo. of unique duplicate mutations:", dup_muts_nu
             , "\n\nExtracting df with unique mutations only"))
  my_df_u = my_df[!duplicated(my_df$mutationinformation),]
}else{
  cat(paste0("\nNo duplicate mutations detected"))
  my_df_u = my_df
}

upos = unique(my_df_u$position)
cat("\nDim of clean df:"); cat(dim(my_df_u))
cat("\nNo. of unique mutational positions:"); cat(length(upos), "\n")
  

###########################
# extract mutations <10Angstroms and symbols
###########################
table(my_df_u$ligand_distance<10)

my_df_u_lig = my_df_u[my_df_u$ligand_distance <10,]


#==================
# Angstroms symbol
#==================

angstroms_symbol = "\u212b"
cat(paste0("There are ", nrow(my_df_u_lig), " sites lying within 10", angstroms_symbol, " of the ligand\n"))

#==================
# Delta symbol
#==================

delta_symbol = "\u0394"; delta_symbol

###########################
# variables for my cols
###########################

mcsm_red2 =  "#ae301e" # most negative
mcsm_red1 =  "#f8766d"

mcsm_mid = "white" # middle

mcsm_blue1 = "#00bfc4"
mcsm_blue2 = "#007d85" # most positive


########################################################################
#               end of data extraction and cleaning for plots          #
########################################################################
# clear variables
rm(opt, spec)