diff --git a/R/ExperimentalDesignSimulation.R b/R/ExperimentalDesignSimulation.R index 3043c7d..a18d01c 100644 --- a/R/ExperimentalDesignSimulation.R +++ b/R/ExperimentalDesignSimulation.R @@ -134,7 +134,7 @@ futureExperimentSimulation = function(N_proteins = 300, } template = readRDS(template1_path) # Use template1 as default - # Note: You could also merge template1 and template3 or choose based on parameters + } # Run simulation diff --git a/R/IC50_prediction.R b/R/IC50_prediction.R index 7189410..2a4859e 100644 --- a/R/IC50_prediction.R +++ b/R/IC50_prediction.R @@ -93,7 +93,11 @@ IC50_upper_bound = NA ) } else { - ic50 = -log10(ic50_est) + if (transform_dose) { + ic50 = -log10(ic50_est) + } else { + ic50 = ic50_est + } if (bootstrap) { if(ratio_response){ @@ -101,14 +105,16 @@ dose = x, response = y, n_samples = n_samples, alpha = alpha, increasing = increasing, - target_response = target_response + target_response = target_response, + transform_x = transform_dose ) } else { bootstrap_res = bootstrapIC50LogScale( x = x, y = y, n_samples = n_samples, alpha = alpha, increasing = increasing, - target_response = target_response + target_response = target_response, + transform_x = transform_dose ) } lower = as.numeric(bootstrap_res$ci_lower_transform) @@ -172,7 +178,7 @@ PredictIC50 = function(fit, target_response = 0.5) { #' @param alpha Significance level for confidence interval (default = 0.10). #' @param increasing Logical. Fit non-decreasing if TRUE. #' @param target_response Numeric value for response level (default = 0.5). -#' @param transform_x Logical. If TRUE, applies log10(dose + 1) transformation. Default = TRUE. +#' @param transform_x Logical. If TRUE, applies log10(x + 1) transformation. Default = TRUE. #' #' @return List with mean IC50, CI bounds, and transformed estimates. #' @importFrom stats quantile @@ -211,14 +217,24 @@ bootstrapIC50 = function(dose, response, n_samples = 1000, alpha = 0.10, ci_upper = stats::quantile(ic50_values_clean, 1 - alpha / 2, na.rm = TRUE) mean_ic50 = mean(ic50_values_clean, na.rm = TRUE) + if (transform_x) { + mean_ic50_transform = -log10(mean_ic50) + ci_lower_transform = -log10(ci_lower) + ci_upper_transform = -log10(ci_upper) + } else { + mean_ic50_transform = mean_ic50 + ci_lower_transform = ci_lower + ci_upper_transform = ci_upper + } + return(list( ic50_values = ic50_vals, mean_ic50 = mean_ic50, ci_lower = ci_lower, ci_upper = ci_upper, - mean_ic50_transform = -log10(mean_ic50), - ci_lower_transform = -log10(ci_lower), - ci_upper_transform = -log10(ci_upper) + mean_ic50_transform = mean_ic50_transform, + ci_lower_transform = ci_lower_transform, + ci_upper_transform = ci_upper_transform )) } @@ -231,11 +247,13 @@ bootstrapIC50 = function(dose, response, n_samples = 1000, alpha = 0.10, #' @param increasing Logical. Fit non-decreasing if TRUE. #' @param target_response Numeric value for response level (default = 0.5). #' +#' @param transform_x Logical. If TRUE, applies log10(x + 1) transformation. Default = TRUE. #' @return List with mean IC50, CI bounds, and transformed estimates. #' @importFrom stats quantile #' @importFrom dplyr filter select mutate group_by summarise arrange distinct bootstrapIC50LogScale = function(x, y, n_samples = 1000, alpha = 0.05, - increasing = FALSE, target_response = 0.5) { + increasing = FALSE, target_response = 0.5, + transform_x = TRUE) { ic50_vals = numeric(n_samples) df = data.frame(dose = x, response = y) @@ -260,7 +278,7 @@ bootstrapIC50LogScale = function(x, y, n_samples = 1000, alpha = 0.05, tryCatch({ fit_sample = fitIsotonicRegression(x_sample, y_sample, increasing = increasing, - transform_x = TRUE, + transform_x = transform_x, ratio_y = FALSE, test_significance = FALSE) ic50_est = PredictIC50(fit_sample, target_response = adjusted_target_response) @@ -275,13 +293,23 @@ bootstrapIC50LogScale = function(x, y, n_samples = 1000, alpha = 0.05, ci_upper = stats::quantile(ic50_values_clean, 1 - alpha / 2, na.rm = TRUE) mean_ic50 = mean(ic50_values_clean, na.rm = TRUE) + if (transform_x) { + mean_ic50_transform = -log10(mean_ic50) + ci_lower_transform = -log10(ci_lower) + ci_upper_transform = -log10(ci_upper) + } else { + mean_ic50_transform = mean_ic50 + ci_lower_transform = ci_lower + ci_upper_transform = ci_upper + } + return(list( ic50_values = ic50_vals, mean_ic50 = mean_ic50, ci_lower = ci_lower, ci_upper = ci_upper, - mean_ic50_transform = -log10(mean_ic50), - ci_lower_transform = -log10(ci_lower), - ci_upper_transform = -log10(ci_upper) + mean_ic50_transform = mean_ic50_transform, + ci_lower_transform = ci_lower_transform, + ci_upper_transform = ci_upper_transform )) } diff --git a/R/Visualize_Isotonic_Fit.R b/R/Visualize_Isotonic_Fit.R index 2c56be7..05b790c 100644 --- a/R/Visualize_Isotonic_Fit.R +++ b/R/Visualize_Isotonic_Fit.R @@ -63,6 +63,7 @@ visualizeResponseProtein = function(data, ratio_response = TRUE, transform_dose = TRUE, show_ic50 = TRUE, + target_response = 0.5, add_ci = FALSE, n_samples = 1000, alpha = 0.10, @@ -96,6 +97,7 @@ visualizeResponseProtein = function(data, ratio_response = ratio_response, transform_dose = transform_dose, increasing = increasing, + target_response = target_response, n_samples = n_samples, alpha = alpha) }) @@ -109,10 +111,12 @@ visualizeResponseProtein = function(data, fit = fit, ratio = ratio_response, show_ic50 = show_ic50, + target_response = target_response, ci = ci_bounds, drug_name = drug_name, protein_name = protein_name, - y_lab = y_lab + y_lab = y_lab, + transform_x = transform_dose ) return(p) @@ -142,6 +146,7 @@ visualizeResponseProtein = function(data, plotIsotonic = function(fit, ratio = TRUE, show_ic50 = FALSE, + target_response = target_response, drug_name = NULL, protein_name = NULL, x_lab = expression(Log[10] ~ "[drug (M)]"), @@ -150,7 +155,8 @@ plotIsotonic = function(fit, ci = NULL, legend = FALSE, theme_style = "classic", - original_label = FALSE) { + original_label = FALSE, + transform_x = TRUE) { # Construct title if not provided if (is.null(title) && !is.null(drug_name) && !is.null(protein_name)) { @@ -165,9 +171,26 @@ plotIsotonic = function(fit, y_lab = expression(Log[2] ~ "Intensity") } + if (transform_x) { + x_lab = expression(Log[10] ~ "[drug (M)]") + } else { + x_lab = "x" + } + # Prepare data - fit_df = data.frame(dose = log10(fit$x), y_pred = fit$y_pred) - orig_df = data.frame(x = log10(fit$x), y = fit$original_y) + #fit_df = data.frame(dose = log10(fit$x), y_pred = fit$y_pred) + #orig_df = data.frame(x = log10(fit$x), y = fit$original_y) + + # Prepare data - fit$x is in the correct scale based on transform_x + fit_df = data.frame(dose = fit$x, y_pred = fit$y_pred) + orig_df = data.frame(x = fit$x, y = fit$original_y) + if (transform_x) { + fit_df = data.frame(dose = log10(fit$x), y_pred = fit$y_pred) + orig_df = data.frame(x = log10(fit$x), y = fit$original_y) + } else { + fit_df = data.frame(dose = fit$x, y_pred = fit$y_pred) + orig_df = data.frame(x = fit$x, y = fit$original_y) + } # Base plot model_fit_plot = ggplot2::ggplot() + @@ -197,10 +220,19 @@ plotIsotonic = function(fit, # Optional IC50 if (show_ic50) { if (ratio) { - target_response = 0.5 + target_response = target_response } else { - dmso_idx = which(orig_df$x == -Inf) - target_response = mean(orig_df$y[dmso_idx], na.rm = TRUE) - 1 + if (transform_x) { + #dmso_idx = which(orig_df$x == 0) # log10(0+1) = 0 + dmso_idx = which(is.infinite(orig_df$x) & orig_df$x < 0) + target_response = mean(orig_df$y[dmso_idx], na.rm = TRUE) + log2(target_response) + + } else { + dmso_idx = which(fit$original_x == 0) # Raw dose = 0 + target_response = mean(orig_df$y[dmso_idx], na.rm = TRUE) + log2(target_response) + } + #dmso_idx = which(orig_df$x == -Inf) + #target_response = mean(orig_df$y[dmso_idx], na.rm = TRUE) - 1 } ic50_pred = PredictIC50(fit, target_response = target_response) @@ -220,29 +252,55 @@ plotIsotonic = function(fit, } #y_ic50 = stats::approx(fit$x, fit$y_pred, xout = ic50_pred)$y - ic50_pred_transform = -log10(ic50_pred) #pIC50 + + if (transform_x) { + # ic50_pred is on log10(dose+1) scale, convert to pIC50 for label + ic50_pred_transform = -log10(ic50_pred) #pIC50 + ic50_label = paste("pIC50 =", round(ic50_pred_transform, 2)) + x_pos = log10(ic50_pred) # Already on correct scale for plotting + } else { + # ic50_pred is on raw scale, use as-is + ic50_label = paste("IC50 =", round(ic50_pred, 2)) + x_pos = ic50_pred + } + + # model_fit_plot = model_fit_plot + + # ggplot2::geom_point(ggplot2::aes(x = log10(ic50_pred), y = y_ic50), + # shape = 23, size = 3.5, fill = "red", color = "black") + + # ggplot2::annotate("text", x = log10(ic50_pred) + 0.35, y = y_ic50, + # label = paste("pIC50 =", round(ic50_pred_transform, 2)), + # vjust = -0.5, hjust = 0.5, size = 4, fontface = "bold") model_fit_plot = model_fit_plot + - ggplot2::geom_point(ggplot2::aes(x = log10(ic50_pred), y = y_ic50), + ggplot2::geom_point(ggplot2::aes(x = x_pos, y = y_ic50), shape = 23, size = 3.5, fill = "red", color = "black") + - ggplot2::annotate("text", x = log10(ic50_pred) + 0.35, y = y_ic50, - label = paste("pIC50 =", round(ic50_pred_transform, 2)), + ggplot2::annotate("text", x = x_pos + 0.35, y = y_ic50, + label = ic50_label, vjust = -0.5, hjust = 0.5, size = 4, fontface = "bold") - # Show IC50 confidence interval if provided + # optional display IC50 confidence interval if (show_ic50 && !is.null(ci)) { - ci_lower_M = 10^(-ci$ci_lower) - ci_upper_M = 10^(-ci$ci_upper) + if (transform_x) { + # CI bounds are in pIC50, convert to log10(dose+1) scale for plotting + ci_lower_dose = 10^(-ci$ci_lower) # pIC50 to dose + ci_upper_dose = 10^(-ci$ci_upper) + ci_lower_x = log10(ci_lower_dose) # dose to log10(dose+1) + ci_upper_x = log10(ci_upper_dose) + } else { + # CI bounds are already in raw scale + ci_lower_x = ci$ci_lower + ci_upper_x = ci$ci_upper + } model_fit_plot = model_fit_plot + - ggplot2::geom_segment(ggplot2::aes(x = log10(ci_lower_M), xend = log10(ci_lower_M), + ggplot2::geom_segment(ggplot2::aes(x = ci_lower_x, xend = ci_lower_x, y = 0, yend = target_response), linetype = "dotted", color = "gray40", linewidth = 0.8) + - ggplot2::geom_segment(ggplot2::aes(x = log10(ci_upper_M), xend = log10(ci_upper_M), + ggplot2::geom_segment(ggplot2::aes(x = ci_upper_x, xend = ci_upper_x, y = 0, yend = target_response), linetype = "dotted", color = "gray40", linewidth = 0.8) + ggplot2::annotate("rect", - xmin = log10(ci_lower_M), xmax = log10(ci_upper_M), + xmin = ci_lower_x, xmax = ci_upper_x, ymin = 0, ymax = target_response, alpha = 0.1, fill = "red") } diff --git a/man/VisualizeResponseProtein.Rd b/man/VisualizeResponseProtein.Rd index a877347..7d89203 100644 --- a/man/VisualizeResponseProtein.Rd +++ b/man/VisualizeResponseProtein.Rd @@ -11,6 +11,7 @@ visualizeResponseProtein( ratio_response = TRUE, transform_dose = TRUE, show_ic50 = TRUE, + target_response = 0.5, add_ci = FALSE, n_samples = 1000, alpha = 0.1, diff --git a/man/bootstrapIC50.Rd b/man/bootstrapIC50.Rd index 8b0069f..7d0a535 100644 --- a/man/bootstrapIC50.Rd +++ b/man/bootstrapIC50.Rd @@ -27,7 +27,7 @@ bootstrapIC50( \item{target_response}{Numeric value for response level (default = 0.5).} -\item{transform_x}{Logical. If TRUE, applies log10(dose + 1) transformation. Default = TRUE.} +\item{transform_x}{Logical. If TRUE, applies log10(x + 1) transformation. Default = TRUE.} } \value{ List with mean IC50, CI bounds, and transformed estimates. diff --git a/man/bootstrapIC50LogScale.Rd b/man/bootstrapIC50LogScale.Rd index a3b0248..1ad9ebf 100644 --- a/man/bootstrapIC50LogScale.Rd +++ b/man/bootstrapIC50LogScale.Rd @@ -10,7 +10,8 @@ bootstrapIC50LogScale( n_samples = 1000, alpha = 0.05, increasing = FALSE, - target_response = 0.5 + target_response = 0.5, + transform_x = TRUE ) } \arguments{ @@ -25,6 +26,8 @@ bootstrapIC50LogScale( \item{increasing}{Logical. Fit non-decreasing if TRUE.} \item{target_response}{Numeric value for response level (default = 0.5).} + +\item{transform_x}{Logical. If TRUE, applies log10(x + 1) transformation. Default = TRUE.} } \value{ List with mean IC50, CI bounds, and transformed estimates. diff --git a/man/plotIsotonic.Rd b/man/plotIsotonic.Rd index dd46149..e50d99f 100644 --- a/man/plotIsotonic.Rd +++ b/man/plotIsotonic.Rd @@ -8,6 +8,7 @@ plotIsotonic( fit, ratio = TRUE, show_ic50 = FALSE, + target_response = target_response, drug_name = NULL, protein_name = NULL, x_lab = expression(Log[10] ~ "[drug (M)]"), @@ -16,7 +17,8 @@ plotIsotonic( ci = NULL, legend = FALSE, theme_style = "classic", - original_label = FALSE + original_label = FALSE, + transform_x = TRUE ) } \arguments{ diff --git a/vignettes/MSstatsResponse.Rmd b/vignettes/MSstatsResponse.Rmd index 7420be9..1490d6f 100644 --- a/vignettes/MSstatsResponse.Rmd +++ b/vignettes/MSstatsResponse.Rmd @@ -179,7 +179,7 @@ dia_prepared %>% ## Drug-protein interaction detection -### Understanding Isotonic Regression +### Isotonic Regression Isotonic regression fits a monotonic function to the data without assuming a parametric form. For this example, we have drug inhibitor dose-response data. When treated with inhibitors we typically expect: @@ -362,8 +362,23 @@ p2 <- visualizeResponseProtein( # Combine plots (requires gridExtra) gridExtra::grid.arrange(p1, p2, ncol = 2) ``` +### Optional: use target_response argument to make predictions for any level - +In this example we predict IC25 and generate confidence intervals. Default is 0.50. + +```{r visualize_raw_scale, fig.height=5, fig.width=7} +# Visualize another target +visualizeResponseProtein( + data = dia_prepared, + protein_name = "PROTEIN_B", + drug_name = "Drug1", + ratio_response = TRUE, + show_ic50 = TRUE, + add_ci = TRUE, + transform_dose = TRUE, + target_response = 0.25 +) +``` ## Experimental design planning