# install.packages(c("drc", "ggplot2", "dplyr"))
library(drc)
## Loading required package: MASS
##
## 'drc' has been loaded.
## Please cite R and 'drc' if used for a publication,
## for references type 'citation()' and 'citation('drc')'.
##
## Attaching package: 'drc'
## The following objects are masked from 'package:stats':
##
## gaussian, getInitial
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
##
## select
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# 1. Concentration series and example data
conc_vec <- c(0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30, 60)
# Agonist: response goes up with concentration
agonist_df <- data.frame(
assay_type = "Agonist: upward curve",
conc_uM = rep(conc_vec, each = 3),
response = c(2,4,3, 5,4,6, 8,7,9, 12,15,13, 20,22,19, 35,38,36,
55,58,52, 78,81,76, 92,95,90, 98,101,96, 99,102,97)
)
# Antagonist: response goes down with concentration
antagonist_df <- data.frame(
assay_type = "Antagonist: downward curve",
conc_uM = rep(conc_vec, each = 3),
response = c(100,98,101, 96,97,95, 94,92,93, 88,86,87, 76,73,75, 58,55,57,
38,35,36, 20,18,19, 9,8,10, 4,3,5, 2,1,3)
)
all_df <- bind_rows(agonist_df, antagonist_df)
# 2. Fit a four-parameter Hill curve (LL.4)
# b = Hill slope, c = bottom, d = top, e = EC50/IC50
fit_hill_curve <- function(dat, endpoint_name = "EC50") {
dat <- filter(dat, !is.na(conc_uM), !is.na(response), conc_uM > 0)
fit <- drm(response ~ conc_uM, data = dat,
fct = LL.4(names = c("b", "c", "d", endpoint_name)))
co <- coef(fit)
# R-squared
pred <- predict(fit, newdata = dat)
r2 <- 1 - sum((dat$response - pred)^2) / sum((dat$response - mean(dat$response))^2)
# Smooth curve for plotting
pred_df <- data.frame(
conc_uM = 10^seq(log10(min(dat$conc_uM)), log10(max(dat$conc_uM)), length.out = 300)
)
pred_df$response <- predict(fit, newdata = pred_df)
list(
pred_df = pred_df,
hill_signed = unname(co[1]), # b (upward often negative, downward often positive)
hill_abs = abs(unname(co[1])),
endpoint_val = unname(co[4]), # EC50 / IC50
r_squared = r2
)
}
agonist_fit <- fit_hill_curve(agonist_df, "EC50")
antagonist_fit <- fit_hill_curve(antagonist_df, "IC50")
# 3. Fitted curves for plotting
pred_all <- bind_rows(
mutate(agonist_fit$pred_df, assay_type = "Agonist: upward curve"),
mutate(antagonist_fit$pred_df, assay_type = "Antagonist: downward curve")
)
# 4. Annotation table: one column per metric + combined label column
annotation_df <- data.frame(
assay_type = c("Agonist: upward curve", "Antagonist: downward curve"),
endpoint_name = c("EC50", "IC50"),
endpoint_val = c(agonist_fit$endpoint_val, antagonist_fit$endpoint_val),
hill_signed = c(agonist_fit$hill_signed, antagonist_fit$hill_signed),
hill_abs = c(agonist_fit$hill_abs, antagonist_fit$hill_abs),
r_squared = c(agonist_fit$r_squared, antagonist_fit$r_squared),
x = 0.0015,
y = 95
)
# Build combined label from the split columns
annotation_df$label <- sprintf(
"%s = %.2f uM\nHill slope = %.2f\n|Hill slope| = %.2f\nR2 = %.2f",
annotation_df$endpoint_name, annotation_df$endpoint_val,
annotation_df$hill_signed, annotation_df$hill_abs,
annotation_df$r_squared
)
print(annotation_df)
## assay_type endpoint_name endpoint_val hill_signed hill_abs
## 1 Agonist: upward curve EC50 0.8359492 -0.7552557 0.7552557
## 2 Antagonist: downward curve IC50 0.4595101 0.7136681 0.7136681
## r_squared x y
## 1 0.9968642 0.0015 95
## 2 0.9990373 0.0015 95
## label
## 1 EC50 = 0.84 uM\nHill slope = -0.76\n|Hill slope| = 0.76\nR2 = 1.00
## 2 IC50 = 0.46 uM\nHill slope = 0.71\n|Hill slope| = 0.71\nR2 = 1.00
write.csv(annotation_df, "annotation_df.csv", row.names = FALSE)
# 5. Plot data points + fitted curves
p <- ggplot() +
geom_point(data = all_df, aes(conc_uM, response), size = 2, alpha = 0.7) +
geom_line(data = pred_all, aes(conc_uM, response), linewidth = 1) +
geom_label(data = annotation_df, aes(x, y, label = label),
hjust = 0, vjust = 1, size = 3.5) +
scale_x_log10() +
facet_wrap(~ assay_type) +
labs(x = "Concentration (uM, log scale)", y = "Response (%)",
title = "Four-parameter Hill curve fitting") +
theme_bw()
print(p)
# ggsave("Hill_slope_example.png", p, width = 8, height = 4, dpi = 300)