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

# FIXME: ADD distance to NA when SP replies
dist_columns = c("ligand_distance", "interface_dist")
DistCutOff = 10
common_cols  = c("mutationinformation"
                 , "X5uhc_position"
                 , "X5uhc_offset"
                 , "position"
                 , "dst_mode"
                 , "mutation_info_labels"
                 , "sensitivity", dist_columns )

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

gene_aff_cols    = colnames(df3)[colnames(df3)%in%scaled_cols_affinity]
gene_stab_cols   = colnames(df3)[colnames(df3)%in%scaled_cols_stability]
gene_common_cols = colnames(df3)[colnames(df3)%in%common_cols]

sel_cols = c(gene_common_cols
             , gene_stab_cols
             ,  gene_aff_cols)

#########################################
#df3_plot = df3[, cols_to_extract]
df3_plot = df3[, sel_cols]

######################
#FILTERING HMMMM!
#all dist <10
#for embb this results in 2 muts
######################
df3_affinity_filtered = df3_plot[ (df3_plot$ligand_distance<10 | df3_plot$interface_dist <10),]
df3_affinity_filtered = df3_plot[ (df3_plot$ligand_distance<10 & df3_plot$interface_dist <10),]

c0u = unique(df3_affinity_filtered$position)
length(c0u)

#df = df3_affinity_filtered
##########################################
#NO FILTERING: annotate the whole df
df = df3_plot
sum(is.na(df))
df2 = na.omit(df)   
c0u = unique(df2$position)
length(c0u)

# reassign orig
my_df_raw = df3

#  now subset
df3 = df2
#######################################################
#=================
# affinity effect
#=================
give_col=function(x,y,df=df3){
  df[df$position==x,y]
}

for (i in unique(df3$position) ){
  #print(i)
  biggest     = max(abs(give_col(i,gene_aff_cols)))
  
  df3[df3$position==i,'abs_max_effect'] = biggest
  df3[df3$position==i,'effect_type']= names(
    give_col(i,gene_aff_cols)[which(
      abs(
        give_col(i,gene_aff_cols)
      ) == biggest, arr.ind=T
    )[, "col"]])
  
  # effect_name = unique(df3[df3$position==i,'effect_type'])
  effect_name = df3[df3$position==i,'effect_type'][1] # pick first one in case we have multiple exact values
  
  ind = rownames(which(abs(df3[df3$position==i,c('position',effect_name)][effect_name])== biggest, arr.ind=T))
  df3[df3$position==i,'effect_sign'] = sign(df3[effect_name][ind,])
}

df3$effect_type = sub("\\.[0-9]+", "", df3$effect_type) # cull duplicate effect types that happen when there are exact duplicate values
df3U = df3[!duplicated(df3[c('position')]), ]
table(df3U$effect_type)
#########################################################
#%% consider stability as well
df4 = df2

#=================
# stability + affinity effect
#=================
effect_cols = c(gene_aff_cols, gene_stab_cols) 

give_col=function(x,y,df=df4){
  df[df$position==x,y]
}

for (i in unique(df4$position) ){
  #print(i)
  biggest     = max(abs(give_col(i,effect_cols)))
  
  df4[df4$position==i,'abs_max_effect'] = biggest
  df4[df4$position==i,'effect_type']= names(
    give_col(i,effect_cols)[which(
      abs(
        give_col(i,effect_cols)
      ) == biggest, arr.ind=T
    )[, "col"]])
  
  # effect_name = unique(df4[df4$position==i,'effect_type'])
  effect_name = df4[df4$position==i,'effect_type'][1] # pick first one in case we have multiple exact values
  
  ind = rownames(which(abs(df4[df4$position==i,c('position',effect_name)][effect_name])== biggest, arr.ind=T))
  df4[df4$position==i,'effect_sign'] = sign(df4[effect_name][ind,])
}

df4$effect_type = sub("\\.[0-9]+", "", df4$effect_type) # cull duplicate effect types that happen when there are exact duplicate values
df4U = df4[!duplicated(df4[c('position')]), ]
table(df4U$effect_type)

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