#!/usr/bin/env Rscript
#########################################################
# TASK: prediction_ps

#=======================================================================
# working dir and loading libraries
getwd()
setwd("~/git/LSHTM_analysis/scripts/")
getwd()

source("plotting/combining_dfs_plotting.R")

#=======
# output
#=======
ps_unadjusted = paste0(outdir, "/results/", tolower(gene), "_unadjusted_logistic_PS.csv")
ps_adjusted = paste0(outdir, "/results/", tolower(gene), "_adjusted_logistic_PS.csv")

####################################################################
#               end of loading libraries and functions              
####################################################################

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# ps
table(merged_df3$mutation_info)
merged_df3$mutation_info_labels = ifelse(merged_df3$mutation_info == dr_muts_col, 1, 0)
table(merged_df3$mutation_info_labels)

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

###############################################################################
#========================================
# merged_df3: UNadjusted,loop
#========================================
my_ivs = c("or_mychisq", "or_kin", "pval_fisher", "af"
           , "ligand_distance"
           , "rsa"
           , "rd_values"
           , "kd_values"
           , "duet_stability_change" 
           , "duet_scaled"
           , "duet_outcome"
           , "ddg"
           , "foldx_scaled"
           , "foldx_outcome")

ps_logistic_df3  = data.frame()

for( i in my_ivs){
  print(i)
  
  df = data.frame(var_name = NA
                  , number_samples = NA
                  , beta = NA
                  , odds_ratio = NA
                  , pvalue = NA
                  , se = NA
                  , zvalue = NA
                  , ci_lower = NA
                  , ci_upper = NA)
  
  model = glm(mutation_info_labels ~ eval(parse(text=i))
              , data = merged_df3
              , family = "binomial")
  
  var_name = i
  
  number_samples = nobs(model)
  
  beta_logistic = summary(model)$coefficients[2,1]; beta_logistic
  
  or_logistic = exp(summary(model)$coefficients[2,1])
  
  pval_logistic = summary(model)$coefficients[2,4]
  
  se_logistic = summary(model)$coefficients[2,2]
  
  zval_logistic = summary(model)$coefficients[2,3]
  
  ci_mod = exp(confint(model))[2,]
  
  ci_lower_logistic = ci_mod[["2.5 %"]]
  ci_upper_logistic = ci_mod[["97.5 %"]]
  
  print(c(var_name, beta_logistic, or_logistic, pval_logistic, se_logistic, zval_logistic, ci_mod))
  
  df$var_name = var_name 
  df$number_samples = number_samples
  df$beta = beta_logistic
  df$odds_ratio = or_logistic
  df$pvalue = pval_logistic
  df$se = se_logistic
  df$zvalue = zval_logistic
  df$ci_lower = ci_lower_logistic 
  df$ci_upper = ci_upper_logistic 
  
  print(df)
  ps_logistic_df3 = rbind(ps_logistic_df3, df)
  
}
#--------------------
# formatting df
#--------------------
ps_logistic_df3$data_source = "df3"
ps_logistic_df3$model = "unadjusted"

ps_logistic_df3$odds_ratio = round(ps_logistic_df3$odds_ratio, 2)
ps_logistic_df3$ci_lower = round(ps_logistic_df3$ci_lower, 2)
ps_logistic_df3$ci_upper = round(ps_logistic_df3$ci_upper, 2)

# adding pvalue_signif
ps_logistic_df3$pvalue_signif = ps_logistic_df3$pvalue
str(ps_logistic_df3$pvalue_signif)

ps_logistic_df3 = dplyr::mutate(ps_logistic_df3
                                , pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
                                                            , pvalue_signif <=0.0001 ~ '****'
                                                            , pvalue_signif <=0.001 ~ '***'
                                                            , pvalue_signif <=0.01 ~ '**'
                                                            , pvalue_signif <0.05 ~ '*'
                                                            , TRUE ~ 'ns'))
# rearranging columns
ps_logistic_df3_o = ps_logistic_df3[c("var_name"
                                      , "number_samples"
                                      , "model"
                                      , "odds_ratio"
                                      , "pvalue"
                                      , "pvalue_signif"
                                      , "beta"
                                      , "se"
                                      , "zvalue"
                                      , "ci_lower"
                                      , "ci_upper"
                                      , "data_source")]
# writing file
write.csv(ps_logistic_df3_o, ps_unadjusted, row.names = F)

#========================================
# merged_df3: adjusted, loop
#========================================
#model_adjusted_df3 = glm(mutation_info_labels ~ or_mychisq + rsa + rd_values + kd_values +
#                       duet_stability_change + duet_scaled + duet_outcome + ddg + foldx_scaled + foldx_outcome
#            , data = merged_df3
#           , family = "binomial")

model_adjusted_df3 = glm(mutation_info_labels ~ rd_values +
                           ligand_distance + duet_stability_change 
                         , data = merged_df3
                         , family = "binomial");summary(model_adjusted_df3)

var_names_df = as.data.frame(names(model_adjusted_df3$coefficients))
names(var_names_df) = c("var_name")

ci_mod = exp(confint(model_adjusted_df3))
ci_mod_df = as.data.frame(ci_mod)
names(ci_mod_df) = c("ci_lower", "ci_upper")
ci_mod_df$ci_lower = round(ci_mod_df$ci_lower, 2)
ci_mod_df$ci_upper = round(ci_mod_df$ci_upper, 2)


estimates_df = as.data.frame(summary(model_adjusted_df3)$coefficients)
names(estimates_df) = c("beta", "se", "zvalue", "pvalue")
estimates_df$odds_ratio = round(exp(estimates_df$beta), 2)
number_samples = nobs(model_adjusted_df3)
estimates_df$number_samples = number_samples
estimates_df$data_source = "df3"
estimates_df$model = "adjusted"

names(ps_logistic_adjusted_df3)

if ( all(rownames(estimates_df) == rownames(ci_mod_df)) ){
  cat("PASS: rownames match. Preparing to merge...")
  ps_logistic_adjusted_df3 = merge(estimates_df, ci_mod_df, by = "row.names", all = T)
}

colnames(ps_logistic_adjusted_df3)[1] <- c("var_name")
d2 = which(ps_logistic_adjusted_df3$var_name ==  "(Intercept)")
ps_logistic_adjusted_df3 = ps_logistic_adjusted_df3[-d2,]

# adding pvalue_signif
ps_logistic_adjusted_df3$pvalue_signif = ps_logistic_adjusted_df3$pvalue
str(ps_logistic_adjusted_df3$pvalue_signif)

ps_logistic_adjusted_df3 = dplyr::mutate(ps_logistic_adjusted_df3
                                         , pvalue_signif = case_when(pvalue_signif == 0.05 ~ "."
                                                                     , pvalue_signif <=0.0001 ~ '****'
                                                                     , pvalue_signif <=0.001 ~ '***'
                                                                     , pvalue_signif <=0.01 ~ '**'
                                                                     , pvalue_signif <0.05 ~ '*'
                                                                     , TRUE ~ 'ns'))
# rearranging columns
colnames(ps_logistic_adjusted_df3)
ps_logistic_adjusted_df3_o = ps_logistic_adjusted_df3[c("var_name"    
                                                        , "number_samples"
                                                        , "model"
                                                        , "odds_ratio"
                                                        , "pvalue"
                                                        , "pvalue_signif"
                                                        , "beta"
                                                        , "se"
                                                        , "zvalue"
                                                        , "ci_lower"
                                                        , "ci_upper"
                                                        ,"data_source")]
# writing file
write.csv(ps_logistic_adjusted_df3_o, ps_adjusted, row.names = F)
###############################################################################