# Saturation binding, control vs treated: specific binding (total minus mean nonspecific at each concentration) fitted
# with a one-site hyperbola in each experiment; Kd and Bmax compared across experiments on the log scale (paired by
# experiment). Also the Scatchard estimate for comparison. R 4.3.3, base R.
d <- read.csv("binding.csv")                         # one row per well: experiment, condition, ligand_nM, well, total, nonspecific
f2 <- function(p) formatC(signif(p, 2), format = "g")
specific <- function(x) { nsb <- ave(x$nonspecific, x$ligand_nM, FUN = mean); data.frame(l = x$ligand_nM, s = x$total - nsb) }
fit <- function(x) { y <- specific(x); coef(nls(s ~ bmax * l / (kd + l), data = y, start = list(bmax = max(y$s), kd = median(y$l)))) }
scatchard_kd <- function(x) { y <- specific(x); m <- tapply(y$s, y$l, mean); l <- as.numeric(names(m)); -1 / coef(lm(I(m / l) ~ m))[[2]] }
res <- list()
for (e in sort(unique(d$experiment))) {
  row <- list()
  for (cnd in c("control", "treated")) { x <- d[d$experiment == e & d$condition == cnd, ]; p <- fit(x); row[[cnd]] <- c(p[["bmax"]], p[["kd"]], scatchard_kd(x)) }
  res[[length(res) + 1]] <- row
  cat(sprintf("Experiment %d: control Bmax %.0f fmol/mg, Kd %.2f nM [Scatchard %.2f]; treated Bmax %.0f, Kd %.2f [Scatchard %.2f]\n",
              e, row$control[1], row$control[2], row$control[3], row$treated[1], row$treated[2], row$treated[3]))
}
for (i in 1:2) {
  r <- sapply(res, function(x) log(x$treated[i] / x$control[i])); tt <- t.test(r)
  cat(sprintf("%s, treated / control, one ratio per experiment (paired, n = %d): %.2f (95%% CI %.2f to %.2f), p = %s\n",
              c("Bmax", "Kd")[i], length(r), exp(mean(r)), exp(tt$conf.int[1]), exp(tt$conf.int[2]), f2(tt$p.value)))
}
pkd <- sapply(res, function(x) -log10(x$control[2] * 1e-9))
cat(sprintf("Control pKd: mean %.2f (SD %.2f); geometric mean Kd %.2f nM\n", mean(pkd), sd(pkd), 10^(-mean(pkd)) * 1e9))
# Global fit pooling all wells: a Bmax for each membrane preparation, one Kd per condition (or one shared Kd)
pieces <- expand.grid(condition = c("control", "treated"), experiment = sort(unique(d$experiment)), stringsAsFactors = FALSE)[, 2:1]
all <- do.call(rbind, lapply(seq_len(nrow(pieces)), function(i) {
  y <- specific(d[d$experiment == pieces$experiment[i] & d$condition == pieces$condition[i], ])
  data.frame(y, grp = i, trt = as.integer(pieces$condition[i] == "treated") + 1) }))
ng <- nrow(pieces)
full <- nls(s ~ bmax[grp] * l / (kd[trt] + l), data = all, start = list(bmax = rep(500, ng), kd = c(1, 1)))
red <- nls(s ~ bmax[grp] * l / (kd + l), data = all, start = list(bmax = rep(500, ng), kd = 1))
ssf <- sum(resid(full)^2); ssr <- sum(resid(red)^2); dfe <- nrow(all) - ng - 2; F <- (ssr - ssf) / (ssf / dfe)
kdf <- coef(full)[c("kd1", "kd2")]
cat(sprintf("Global fit of all %d wells (Kd shared across experiments within each condition): Kd treated / control %.2f; extra-sum-of-squares F(1, %d) = %.2f, p = %s\n",
            nrow(all), kdf[[2]] / kdf[[1]], dfe, F, f2(pf(F, 1, dfe, lower.tail = FALSE))))