utils.R

################################################################################
# Utility Functions for IC50 Analysis
#
# This file contains reusable functions for IC50 dose-response curve analysis.
# These utilities are used by multiple scripts in this repository to reduce
# code duplication while keeping the main workflow scripts readable.
#
# Functions:
#   - calculate_ic50(): Fit 4-parameter log-logistic model to dose-response data
#   - create_ic50_plot(): Generate interactive plotly visualization
################################################################################

library(drc)
library(plotly)
library(pracma)

################################################################################
# calculate_ic50
#
# Calculates IC50 (half maximal inhibitory concentration) using a 4-parameter
# log-logistic model (LL.4). This is the most common model for dose-response
# curves in pharmacology.
#
# Parameters:
#   df                - Data frame containing the dose-response data
#   concentration_col - Name of the column containing drug concentrations
#   mortality_col     - Name of the column containing mortality/response values
#
# Returns:
#   A list containing:
#     $ic50       - IC50 value (concentration at 50% response)
#     $slope      - Slope of the dose-response curve
#     $std_error  - Standard error of IC50 estimate
#     $t_value    - t-statistic for IC50
#     $p_value    - p-value for IC50 significance
#     $model      - The fitted drc model object (for advanced users)
#
# Model Details:
#   The LL.4 model has 4 parameters:
#     - Slope: Steepness of the curve
#     - Lower limit: Fixed at 0 (minimum response)
#     - Upper limit: Fixed at 100 (maximum response)
#     - IC50: Concentration at 50% response (what we're solving for)
#
# Example:
#   df <- read.csv("examples/sample_ic50_data.csv")
#   result <- calculate_ic50(df, "Cell.Mortality.Concentration", "Cell.Mortality.Mortality.24h")
#   print(paste("IC50:", result$ic50))
#   print(paste("p-value:", result$p_value))
################################################################################

calculate_ic50 <- function(df, concentration_col, mortality_col) {

  # Validate inputs
  if (!concentration_col %in% names(df)) {
    stop(paste("Column not found:", concentration_col,
               "\nAvailable columns:", paste(names(df), collapse = ", ")))
  }

  if (!mortality_col %in% names(df)) {
    stop(paste("Column not found:", mortality_col,
               "\nAvailable columns:", paste(names(df), collapse = ", ")))
  }

  # Create formula for drc model
  formula <- as.formula(paste(mortality_col, "~", concentration_col))

  # Fit 4-parameter log-logistic model
  # fixed = c(NA, 0, 100, NA) means:
  #   Slope: NA (estimate from data)
  #   Lower limit: 0 (fixed)
  #   Upper limit: 100 (fixed)
  #   IC50: NA (estimate from data)
  tryCatch({
    model <- drm(
      formula,
      data = df,
      fct = LL.4(
        fixed = c(NA, 0, 100, NA),
        names = c("Slope", "LS Inferior", "LS Superior", "IC50")
      )
    )

    # Extract model summary
    summary_result <- summary(model)
    coefficients <- summary_result$coefficients

    # Extract key statistics
    # Coefficients matrix has rows: [1]=Slope, [2]=IC50
    # Columns: [1]=Estimate, [2]=Std.Error, [3]=t-value, [4]=p-value
    result <- list(
      ic50 = round(coefficients[2, 1], 2),      # IC50 estimate
      slope = round(coefficients[1, 1], 2),     # Slope estimate
      std_error = round(coefficients[2, 2], 2), # Std error of IC50
      t_value = round(coefficients[2, 3], 2),   # t-statistic
      p_value = signif(coefficients[2, 4], 3),  # p-value (3 sig figs)
      model = model                              # Full model object
    )

    return(result)

  }, error = function(e) {
    stop(paste(
      "IC50 calculation failed:",
      e$message,
      "\nPlease check:",
      "  1. Data has sufficient range of concentrations",
      "  2. Response values are between 0-100",
      "  3. No missing values in concentration or response columns",
      sep = "\n"
    ))
  })
}

################################################################################
# create_ic50_plot
#
# Creates an interactive plotly visualization of dose-response data with
# fitted IC50 curve. The plot shows:
#   - Scatter points: Original data
#   - Smooth curve: Fitted 4-parameter log-logistic model
#   - Log scale on x-axis (standard for dose-response curves)
#
# Parameters:
#   df                - Data frame containing the data
#   concentration_col - Name of column with drug concentrations
#   mortality_col     - Name of column with mortality/response values
#   ic50_value        - IC50 value (from calculate_ic50 function)
#   slope             - Slope value (from calculate_ic50 function)
#   title             - Optional: Custom title for the plot (default: auto-generated)
#
# Returns:
#   A plotly object that can be:
#     - Displayed interactively in RStudio
#     - Saved with htmlwidgets::saveWidget()
#     - Embedded in Shiny apps
#
# The Logistic Function:
#   The 4-parameter logistic equation is:
#     y = D + (A - D) / (1 + (x/C)^B)
#   Where:
#     A = upper limit (100)
#     B = slope
#     C = IC50 (inflection point)
#     D = lower limit (0)
#
# Example:
#   result <- calculate_ic50(df, "Cell.Mortality.Concentration", "Cell.Mortality.Mortality.24h")
#   plot <- create_ic50_plot(df, "Cell.Mortality.Concentration", "Cell.Mortality.Mortality.24h",
#                            result$ic50, result$slope)
#   htmlwidgets::saveWidget(plot, "ic50_plot.html", selfcontained = TRUE)
################################################################################

create_ic50_plot <- function(df, concentration_col, mortality_col,
                             ic50_value, slope, title = NULL) {

  # Validate inputs
  if (!concentration_col %in% names(df)) {
    stop(paste("Column not found:", concentration_col))
  }

  if (!mortality_col %in% names(df)) {
    stop(paste("Column not found:", mortality_col))
  }

  # Define 4-parameter logistic function
  # This is the mathematical model that describes the dose-response curve
  logistic4 <- function(x, A, B, C, D) {
    return((A - D) / (1.0 + ((x / C) ** B)) + D)
  }

  # Extract data columns
  conc_data <- df[[concentration_col]]
  mort_data <- df[[mortality_col]]

  # Create scatter plot of observed data
  plot <- plot_ly(
    data = df,
    x = conc_data,
    y = mort_data,
    type = "scatter",
    mode = "markers",
    name = "Observed Data",
    marker = list(size = 10, color = 'rgba(50, 120, 200, 0.8)')
  )

  # Generate default title if not provided
  if (is.null(title)) {
    title <- paste("IC50 Dose-Response Curve (IC50 =", ic50_value, ")")
  }

  # Configure plot layout
  plot <- layout(
    plot,
    xaxis = list(
      type = "log",  # Log scale for concentrations (standard practice)
      title = "Concentration"
    ),
    yaxis = list(
      title = paste(mortality_col, "(%)")
    ),
    title = title,
    hovermode = "closest"
  )

  # Generate smooth curve using the fitted model
  # linspace creates evenly spaced points for a smooth curve
  xs <- linspace(min(conc_data), max(conc_data), n = 1000)

  # Calculate fitted y-values using the logistic function
  # Parameters: A=100 (upper), B=slope, C=IC50, D=0 (lower)
  ys <- logistic4(xs, 100, slope, ic50_value, 0)

  # Add fitted curve to plot (line only, no markers)
  plot <- add_trace(
    plot,
    x = xs,
    y = ys,
    type = "scatter",
    mode = "lines",
    name = "Fitted Curve",
    line = list(
      color = 'rgba(200, 50, 50, 0.8)',
      width = 2
    ),
    inherit = FALSE  # Don't inherit data from previous trace
  )

  return(plot)
}

################################################################################
# End of utils.R
#
# Usage Notes:
#   - Always source this file before using the functions: source("utils.R")
#   - Both functions include input validation and helpful error messages
#   - The calculate_ic50() function is designed to work with any column names
#   - The create_ic50_plot() function returns a plotly object, not a static image
#
# For more examples, see the scripts in the examples/ directory
################################################################################