require(tidyverse)
require(dr4pl)

#tested with R version 3.5.1 (2018-07-02)
#tidyverse 1.2.1
#dr4pl 1.1.8



#please load the functions below into your environment
#wrapper function for fitting a 4 parameter logistic model
fit_dr4pl <- function(df, drc_formula = 'response~dose', assay = 'motility', ini_met = 'logistic', rob_met = NULL) {
  
  fit_formula <- eval(parse(text = drc_formula))
  print(fit_formula)
  
  if(assay == 'motiliy') {
    df_trend <- df$datamode %>% unique()
  } else {
    
    df_trend <- 'auto'
    
  }
  save_fit <- safely(~dr4pl::dr4pl(fit_formula, data = df, trend = df_trend, method.ini = ini_met, method.rob = rob_met))
  
  
  fit_result <- df %>% 
    save_fit
  
  print('summary check:')
  
  summary_check <- try(summary(fit_result$result), silent = TRUE)
  
  if (is.null(fit_result$result) | inherits(summary_check, 'try-error')) {
    
    print('did not work')
    ini_met = 'Mead'
    print(paste('trying to fit with ', ini_met, ' method.ini argument'))
    fit_result <- df %>%
      save_fit
    
  }
  if (is.null(fit_result$result)| inherits(summary_check, 'try-error')) {
    
    print('did not work again')
    rob_met = 'absolute'
    print(paste('trying to fit with ', rob_met, ' method.ini argument'))
    fit_result <- df %>%
      save_fit
    
  } else {
    
    print('worked')
    
  }
  return(fit_result)
}

#wrapper function to extract data from dr4pl object
get_info_dr4pl <- function(dr4pl_obj, df = NULL, perform_fit = FALSE) {
  
  if(!perform_fit){
    
    dr4pl_object <- dr4pl_obj
    
  } else {
    
    dr4pl_object <- fit_dr4pl(df)
    
  }
  
  estimate_values <- summary(dr4pl_object) %>% coef() %>% print()
  ec50 <- estimate_values$Estimate[2] %>% round(digits = 2) %>% print()
  UpLi <- estimate_values$Estimate[1] %>% round() %>% print()
  LoLi <- estimate_values$Estimate[4] %>% round() %>% print()
  max_effect <- max(dr4pl_object$data$Response)
  min_effect <- min(dr4pl_object$data$Response)
  hill <- estimate_values$Estimate[3] %>% round(digits = 3) %>% print()
  
  gof_dr4pl <- try(dr4pl::gof.dr4pl(dr4pl_object), silent = TRUE)
  
  
  if(class(gof_dr4pl) != 'try-error') {
    
    #RESIDUAL squared error, can be decomposed into Lack-of-fit error and pure error sum of squares, which are caluclated by gof_dr4pl function of package
    RSS <- sum(gof_dr4pl[3] %>% as.numeric(), gof_dr4pl[4] %>% as.numeric()) %>% print()
    
    #calculate total error
    TSS <- sum((dr4pl_object$data$Response - mean(dr4pl_object$data$Response))^2) %>% print()
    
    gof <- 1 - (RSS/TSS)
    
  } else {
    
    gof <- NA
    
  }
  
  return_df <- tibble(
    EC50 = ec50,
    hill_slope = hill,
    Efficacy_high = UpLi,
    Max_effect = max_effect,
    Mean_effec = mean(dr4pl_object$data$Response),
    Min_effect = min_effect,
    Efficacy_low = LoLi,
    rsquared = gof
  )
  
  return(return_df)
  
}

#wrapper function to extract fitting values
get_fitting_stats <- function(df) {
  
  return_df <- df %>% 
    mutate(fit_data = map(fit_check, get_info_dr4pl)) %>%
    unnest(fit_data) %>% 
    mutate(
      datamode_new = ifelse(sign(hill_slope)==1 & (ceiling(Max_effect-Min_effect) > 25), 'Increasing',
                            ifelse(sign(hill_slope<0) & (ceiling(Max_effect-Min_effect) >= 15), 'Decreasing', 'Inactive')),
      Efficacy = ifelse(datamode_new == 'Decreasing', ceiling(abs(Max_effect-Min_effect)),
                        ifelse(datamode_new == 'Increasing', ceiling(abs(Min_effect-Max_effect)), 0))) %>%
    print()
  
}

#load data; please setup path to csv file
fit_acr <- read_csv('supplementary_file_4.csv') %>% 
  rename(reframe_id = compound_id) %>% 
  group_by(Population, reframe_id) %>%
  #filter out unstained population
  filter(Population != 'Pi-_PNA-') %>%
  nest() %>%
  mutate(fit = map(data, ~fit_dr4pl(., drc_formula = 'ratio~conc_uM')),
         fit_check = map(fit, 'result'))

#select flow cytometry population
fit_acr_filter <- filter(fit_acr, Population == 'Pi+_PNA+')

for (i in seq_len(nrow(fit_acr_filter))) {

  
  sub_fit_df <- fit_acr_filter %>% 
    filter(dplyr::row_number()==i)

  sub_dataset <- sub_fit_df %>%
    unnest(data) %>% print()

  pop <- sub_fit_df %>% 
    pull(Population) %>% print()
  
  compound <- sub_fit_df %>% 
    pull(reframe_id) %>% print()
  
  #extraction from https://gist.github.com/andrewheiss/a418909911fad3d65668c8c219fe62db
  mo <- sub_fit_df %>% 
    '[['('fit_check') %>%
    .[[1]] %>% print()
  
  
  estimate_values <- summary(mo) %>% coef() %>% print()
  
  EC50 <- estimate_values$Estimate[2] %>% round(digits = 3) %>% { ifelse(.<10, format(., nsmall = 3), '>10') } %>% print()
  UpLi <- estimate_values$Estimate[1] %>% round(digits = 1) %>% { ifelse(.<10,format(., nsmall = 1), '>10') } %>%  print()
  LoLi <- estimate_values$Estimate[4] %>% round(digits = 1) %>% { ifelse(.>0, format(., nsmall = 1), '<0.0' ) } %>%  print()
  hillSlope <- estimate_values$Estimate[3] %>% round(digits = 3) %>% format(nsmall = 3) %>% print()
  #compound <- sub_dataset$Chemical_name %>% unique()
  y_label <- paste('%', gsub('_', '', pop)) %>% print()

  control_stats <- sub_dataset %>%
    group_by(reframe_id) %>% 
    summarize(neg_mean = mean(neg_med),
              neg_sd = sd(neg_mad)) %>%  print()
  
  
  
  ggplot(mo$data, aes(Dose, Response)) +
    stat_summary(fun.data = 'mean_sdl', fun.args = list(mult = 1), geom = 'errorbar', color = 'darkgrey',
                 width = 0.2) +
    stat_summary(fun.y = 'mean', geom = 'point', color = 'black') +
    stat_function(fun = dr4pl::MeanResponse, args = list(theta = mo$parameters)) +
    geom_hline(yintercept = control_stats$neg_mean, size = 0.5, linetype = 'dashed', color = 'red') +
    scale_x_log10(breaks = c(0.001, 0.01, 0.1, 1, 10), limits =  c(0.0007, 15),
                  labels = c('0.001', '0.01', '0.1', '1', '10')) +
    scale_y_continuous(breaks = seq(0, 100, 20), limits = c(-10,110)) +
    labs(x = expression('Concentration ['*mu*'M]'), 
         y = y_label,
         title = compound,
         subtitle = paste('Hill slope:', hillSlope, '\t','EC50:', EC50, '\u03BCM')) +
    theme(
      panel.grid.minor = element_blank(),
      plot.title = element_text(face = 'bold', size = 12),
      plot.subtitle = element_text(face = 'italic', size = 10)
    )
  ggsave(paste('Reframe', pop, compound, 'fit.png', sep = '_'), width = 2.5, height = 2.5, scale = 1.3)
  
  
}


reframe_acr_sum <- fit_acr %>% 
  group_by(Population, reframe_id) %>% 
  nest() %>% 
  mutate(fit_data = map(data, get_fitting_stats)) %>% 
  print()

#use Pi+ PNA+ for tyrothricin
tyrothricin <- reframe_acr_sum %>% 
  filter(reframe_id == 'RFM-011-941-5', Population == 'Pi+_PNA+')

remaining_cpds <- reframe_acr_sum %>% 
  filter(Population == 'Pi-_PNA+', reframe_id != 'RFM-011-941-5')

final_reframe <- bind_rows(tyrothricin, remaining_cpds) %>% 
  unnest(fit_data) %>% print() %>% 
  select(-data, -data1, -fit, -fit_check, -hill_slope:-Efficacy_low) %>%
  mutate(EC50 = ifelse(EC50 > 10, '> 10', EC50)) %>% print()

