Skip to content

Main #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
sszvetecz merged 2 commits into from
Feb 10, 2026
Merged

Main #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
eb1be23
added custom target response to visualization
sszvetecz Nov 3, 2025
9e2a4e4
fix x-axis scale and add custom target_response for plotting
sszvetecz Feb 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion R/ExperimentalDesignSimulation.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
52 changes: 40 additions & 12 deletions R/IC50_prediction.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,22 +93,28 @@
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){
bootstrap_res = bootstrapIC50(
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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
))
}

Expand All @@ -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)

Expand All @@ -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)
Expand All @@ -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
))
}
92 changes: 75 additions & 17 deletions R/Visualize_Isotonic_Fit.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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)
})
Expand All @@ -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)
Expand Down Expand Up @@ -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)]"),
Expand All @@ -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)) {
Expand All @@ -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() +
Expand Down Expand Up @@ -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)
Expand All @@ -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")
}
Expand Down
1 change: 1 addition & 0 deletions man/VisualizeResponseProtein.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

2 changes: 1 addition & 1 deletion man/bootstrapIC50.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

5 changes: 4 additions & 1 deletion man/bootstrapIC50LogScale.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

4 changes: 3 additions & 1 deletion man/plotIsotonic.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

19 changes: 17 additions & 2 deletions vignettes/MSstatsResponse.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand Down Expand Up @@ -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

Expand Down