#source("~/git/LSHTM_analysis/config/pnca.R")
#source("~/git/LSHTM_analysis/config/alr.R")
#source("~/git/LSHTM_analysis/config/gid.R")
#source("~/git/LSHTM_analysis/config/embb.R")
#source("~/git/LSHTM_analysis/config/katg.R")
source("~/git/LSHTM_analysis/config/rpob.R")

source("/home/tanu/git/LSHTM_analysis/my_header.R")
#########################################################
# TASK: Generate averaged affinity values 
# across all affinity tools for a given structure
# as applicable...
#########################################################

#=======
# output
#=======
outdir_images = paste0("~/git/Writing/thesis/images/results/", tolower(gene))

#OutFile1
outfile_mean_aff = paste0(outdir_images, "/", tolower(gene)
                                , "_mean_affinity_all.csv")
print(paste0("Output file:", outfile_mean_aff))

#OutFile2
outfile_mean_aff_priorty = paste0(outdir_images, "/", tolower(gene)
                          , "_mean_affinity_priority.csv")
print(paste0("Output file:", outfile_mean_aff_priorty))

#%%===============================================================

#=============
# Input
#=============
df3_filename = paste0("/home/tanu/git/Data/", drug, "/output/", tolower(gene), "_merged_df3.csv")
df3 = read.csv(df3_filename)
length(df3$mutationinformation)

# mut_info checks
table(df3$mutation_info)
table(df3$mutation_info_orig)
table(df3$mutation_info_labels_orig)

# used in plots and analyses
table(df3$mutation_info_labels) # different, and matches dst_mode
table(df3$dst_mode)

# create column based on dst mode with different colname
table(is.na(df3$dst))
table(is.na(df3$dst_mode))

#===============
# Create column: sensitivity mapped to dst_mode
#===============
df3$sensitivity = ifelse(df3$dst_mode == 1, "R", "S")
table(df3$sensitivity)

length(unique((df3$mutationinformation)))
all_colnames = as.data.frame(colnames(df3))
common_cols  = c("mutationinformation"
                 , "position"
                 , "dst_mode"
                 , "mutation_info_labels"
                 , "sensitivity"
                 , "ligand_distance")

all_colnames$`colnames(df3)`[grep("scaled", all_colnames$`colnames(df3)`)]
all_colnames$`colnames(df3)`[grep("outcome", all_colnames$`colnames(df3)`)]

#===================
# stability cols
#===================
raw_cols_stability =  c("duet_stability_change"
                        , "deepddg"
                        , "ddg_dynamut2"
                        , "ddg_foldx")

scaled_cols_stability = c("duet_scaled"       
                          , "deepddg_scaled"   
                          , "ddg_dynamut2_scaled"
                          , "foldx_scaled")

outcome_cols_stability = c("duet_outcome"
                           , "deepddg_outcome"
                           , "ddg_dynamut2_outcome"
                           , "foldx_outcome")

#===================
# affinity cols
#===================
raw_cols_affinity =  c("ligand_affinity_change"
                       , "mmcsm_lig"
                       , "mcsm_ppi2_affinity"
                       , "mcsm_na_affinity")

scaled_cols_affinity = c("affinity_scaled" 
                         , "mmcsm_lig_scaled" 
                         , "mcsm_ppi2_scaled" 
                         , "mcsm_na_scaled" )

outcome_cols_affinity  = c( "ligand_outcome"
                            , "mmcsm_lig_outcome"
                            , "mcsm_ppi2_outcome"
                            , "mcsm_na_outcome")

#===================
# conservation cols
#===================
# raw_cols_conservation =  c("consurf_score"
#                            , "snap2_score"
#                            , "provean_score")
# 
# scaled_cols_conservation = c("consurf_scaled"
#                              , "snap2_scaled"
#                              , "provean_scaled")
# 
# # CANNOT strictly be used, as categories are not identical with conssurf missing altogether
# outcome_cols_conservation = c("provean_outcome"
#                               , "snap2_outcome"
#                               #consurf outcome doesn't exist
# )

######################################################################
cols_to_consider = colnames(df3)[colnames(df3)%in%c(common_cols
                                                    , raw_cols
                                                    , scaled_cols
                                                    , outcome_cols_affinity)]

# cols_to_extract  = cols_to_consider[cols_to_consider%in%c(common_cols
#                                                           , outcome_cols_affinity)]
##############################################################
#####################
# Ensemble affinity: affinity_cols
#####################
# extract outcome cols and map numeric values to the categories
# Destabilising == 0, and stabilising == 1 so rescaling can let -1 be destabilising

# df3_plot = df3[, cols_to_extract]
# 
# df3_plot[, outcome_cols_affinity] <- sapply(df3_plot[, outcome_cols_affinity]
#                              , function(x){ifelse(x == "Destabilising", 0, 1)})

df3_plot = df3[, c(common_cols, scaled_cols)]
#=====================================
# Affintiy (2 cols): average the scores
# across predictors ==> average by
# position ==> scale b/w -1 and 1

# column to average: ens_affinity
#=====================================
cols_to_average_affinity = which(colnames(df3_plot)%in%outcome_cols_affinity)
cols_to_average_affinity

# ensemble average across predictors
df3_plot_affinity$ens_affinity = rowMeans(df3_plot_affinity[,cols_to_average_affinity])

head(df3_plot_affinity$position); head(df3_plot_affinity$mutationinformation)
head(df3_plot_affinity$ens_affinity)
table(df3_plot_affinity$ens_affinity)

# ensemble average of predictors by position
mean_ens_affinity_by_position <- df3_plot_affinity %>%
  dplyr::group_by(position) %>%
  dplyr::summarize(avg_ens_affinity = mean(ens_affinity))

# REscale b/w -1 and 1
#en_aff_min = min(mean_ens_affinity_by_position['ens_affinity'])
#en_aff_max = max(mean_ens_affinity_by_position['ens_affinity']) 

# scale the average affintiy value between -1 and 1
# mean_ens_affinity_by_position['avg_ens_affinity_scaled'] = lapply(mean_ens_affinity_by_position['avg_ens_affinity']
#                                                        , function(x) ifelse(x < 0, x/abs(en_aff_min), x/en_aff_max))


mean_ens_affinity_by_position['avg_ens_affinity_scaled'] = lapply(mean_ens_affinity_by_position['avg_ens_affinity']
                                                                  , function(x) {
                                                                    scales::rescale(x, to  = c(-1,1)
                                                                                    #, from = c(en_aff_min,en_aff_max))
                                                                                    , from = c(0,1))
                                                                  })
cat(paste0('Average affintiy scores:\n'
           , head(mean_ens_affinity_by_position['avg_ens_affinity'])
           , '\n---------------------------------------------------------------'
           , '\nAverage affintiy scaled scores:\n'
           , head(mean_ens_affinity_by_position['avg_ens_affinity_scaled'])))

#convert to a df
mean_ens_affinity_by_position = as.data.frame(mean_ens_affinity_by_position)


#FIXME: sanity checks
# TODO: predetermine the bounds
# l_bound_ens_aff = min(mean_ens_affintiy_by_position['avg_ens_affinity_scaled'])
# u_bound_ens_aff = max(mean_ens_affintiy_by_position['avg_ens_affinity_scaled'])
# 
# if ( (l_bound_ens_aff == -1) && (u_bound_ens_aff == 1) ){
#   cat(paste0("PASS: ensemble affinity scores averaged by position and then scaled"
#         , "\nmin ensemble averaged affinity: ", l_bound_ens_aff
#         , "\nmax ensemble averaged affinity: ", u_bound_ens_aff))
# }else{
#   cat(paste0("FAIL: ensemble affinity scores could not be scaled b/w -1 and 1"
#         , "\nmin ensemble averaged affinity: ", l_bound_ens_aff
#         , "\nmax ensemble averaged affinity: ", u_bound_ens_aff))
#   quit()
# } 


######################################################################
##################
# merge: mean ensemble stability and affinity by_position
####################
# if ( class(mean_ens_stability_by_position) && class(mean_ens_affinity_by_position) != "data.frame"){
#   cat("Y")
# }

common_cols = intersect(colnames(mean_ens_stability_by_position), colnames(mean_ens_affinity_by_position))

if (dim(mean_ens_stability_by_position) && dim(mean_ens_affinity_by_position)){
  print(paste0("PASS: dim's match, mering dfs by column :", common_cols))
  #combined = as.data.frame(cbind(mean_duet_by_position, mean_affinity_by_position ))
  combined_df = as.data.frame(merge(mean_ens_stability_by_position
                                    , mean_ens_affinity_by_position
                                    , by = common_cols
                                    , all = T))
  
  cat(paste0("\nnrows combined_df:", nrow(combined_df)
               , "\nnrows combined_df:", ncol(combined_df)))
}else{
    cat(paste0("FAIL: dim's mismatch, aborting cbind!"
          , "\nnrows df1:", nrow(mean_duet_by_position)
          , "\nnrows df2:", nrow(mean_affinity_by_position)))
    quit()      
}
#%%============================================================
# output
write.csv(combined_df, outfile_mean_ens_st_aff
          , row.names = F)
cat("Finished writing file:\n"
    , outfile_mean_ens_st_aff
    , "\nNo. of rows:", nrow(combined_df)
    , "\nNo. of cols:", ncol(combined_df))

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