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

#########################################################
# TASK:

#########################################################
source("Header_TT.R")
require(data.table)
require(dplyr)
#========================================================
#	Read file: call script for combining df for PS
#source("../combining_two_df.R")

#???????????

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

#=======
# output
#=======
out_filename_mean_stability = paste0(tolower(gene), "_mean_stability.csv") 
outfile_mean_stability = paste0(outdir, "/", out_filename_mean_stability)
print(paste0("Output file:", outfile_mean_stability))

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

my_df = read.csv(infile_params
                 #, stringsAsFactors = F
                 , header = T)

cat('Input dimensions:', dim(my_df)) 

# clear variables
rm(in_filename_params, infile_params)

# quick checks
colnames(my_df)
str(my_df)

# check for duplicate mutations
if ( length(unique(my_df$mutationinformation)) != length(my_df$mutationinformation)){
  cat(paste0('CAUTION:', ' 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('No duplicate mutations detected'))
  my_df_u = my_df
}

upos = unique(my_df_u$position)
cat('Dim of clean df:')
cat(dim(my_df_u))
cat('\nNo. of unique mutational positions:'); cat(length(upos))

########################################################################
#               end of data extraction and cleaning for plots          #
########################################################################
#================
# Data for plots
#================
# REASSIGNMENT as necessary
df  = my_df_u
rm(my_df)

###########################
# Data for bfactor figure
# PS average 
# Lig average
###########################
head(df$position); head(df$mutationinformation)
head(df$duet_scaled)

# order data frame 
#df = df[order(df$position),] #already done

head(df$position); head(df$mutationinformation)
head(df$duet_scaled)

#***********
# PS: average by position
#***********
mean_duet_by_position <- df %>%
  group_by(position) %>%
  summarize(averaged.duet = mean(duet_scaled))

#***********
# Lig: average by position
#***********
mean_affinity_by_position <- df %>%
  group_by(position) %>%
  summarize(averaged.affinity = mean(affinity_scaled))

#***********
# cbind:mean_duet_by_position and mean_affinity_by_position
#***********

combined = as.data.frame(cbind(mean_duet_by_position, mean_affinity_by_position ))

# sanity check
# mean_PS_affinity_Bfactor

colnames(combined)

colnames(combined) = c("position"
                       , "average_duet_scaled"
                       , "position2"
                       , "average_affinity_scaled")

colnames(combined)

identical(combined$position, combined$position2)

n = which(colnames(combined) == "position2"); n

combined_df = combined[,-n]

max(combined_df$average_duet_scaled) ; min(combined_df$average_duet_scaled)

max(combined_df$average_affinity_scaled) ; min(combined_df$average_affinity_scaled)

head(combined_df$position); tail(combined_df$position)
#%%============================================================
# output
write.csv(combined_df, outfile_mean_stability
          , row.names = F)
cat("Finished writing file:\n"
    , outfile_mean_stability
    , "\nNo. of rows:", nrow(combined_df)
    , "\nNo. of cols:", ncol(combined_df))

# end of script
#===============================================================