# 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. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats
from scipy.optimize import curve_fit

d = pd.read_csv("binding.csv")                       # one row per well: experiment, condition, ligand_nM, well, total, nonspecific
one_site = lambda l, bmax, kd: bmax * l / (kd + l)
def specific(x):                                     # total wells minus the mean nonspecific binding at that concentration
    nsb = x.groupby("ligand_nM").nonspecific.transform("mean"); return x.ligand_nM.to_numpy(), (x.total - nsb).to_numpy()
def fit(x):
    l, s = specific(x); p, _ = curve_fit(one_site, l, s, p0=[s.max(), np.median(l)]); return p
def scatchard_kd(x):                                 # bound/free against bound; slope = -1/Kd
    l, s = specific(x); m = pd.DataFrame({"l": l, "s": s}).groupby("l").s.mean()
    slope, _ = np.polyfit(m.to_numpy(), (m / m.index).to_numpy(), 1); return -1 / slope
res = []
for e in sorted(d.experiment.unique()):
    row = {}
    for c in ("control", "treated"):
        x = d[(d.experiment == e) & (d.condition == c)]; bmax, kd = fit(x); row[c] = (bmax, kd, scatchard_kd(x))
    res.append(row)
    print(f"Experiment {e}: control Bmax {row['control'][0]:.0f} fmol/mg, Kd {row['control'][1]:.2f} nM [Scatchard {row['control'][2]:.2f}]; "
          f"treated Bmax {row['treated'][0]:.0f}, Kd {row['treated'][1]:.2f} [Scatchard {row['treated'][2]:.2f}]")
for i, lab in ((0, "Bmax"), (1, "Kd")):
    r = np.array([np.log(x["treated"][i] / x["control"][i]) for x in res]); t = stats.ttest_1samp(r, 0); ci = t.confidence_interval()
    print(f"{lab}, treated / control, one ratio per experiment (paired, n = {len(r)}): {np.exp(r.mean()):.2f} (95% CI {np.exp(ci.low):.2f} to {np.exp(ci.high):.2f}), p = {t.pvalue:.2g}")
pkd = np.array([-np.log10(x["control"][1] * 1e-9) for x in res])
print(f"Control pKd: mean {pkd.mean():.2f} (SD {pkd.std(ddof=1):.2f}); geometric mean Kd {10 ** (-pkd.mean()) * 1e9:.2f} nM")
# Global fit pooling all wells: a Bmax for each membrane preparation, one Kd per condition (or one shared Kd)
pieces = [(e, c) for e in sorted(d.experiment.unique()) for c in ("control", "treated")]
L_all, S_all, G_all, T_all = [], [], [], []
for gi, (e, c) in enumerate(pieces):
    l, s = specific(d[(d.experiment == e) & (d.condition == c)]); L_all.append(l); S_all.append(s)
    G_all.append(np.full(len(l), gi)); T_all.append(np.full(len(l), int(c == "treated")))
L_all, S_all, G_all, T_all = map(np.concatenate, (L_all, S_all, G_all, T_all)); ng = len(pieces)
full = lambda X, *q: np.array(q)[X[1].astype(int)] * X[0] / (np.array(q)[ng + X[2].astype(int)] + X[0])
red = lambda X, *q: np.array(q)[X[1].astype(int)] * X[0] / (q[ng] + X[0])
X = np.vstack([L_all, G_all, T_all])
pf, _ = curve_fit(full, X, S_all, p0=[500] * ng + [1, 1]); pr, _ = curve_fit(red, X, S_all, p0=[500] * ng + [1])
ssf, ssr = np.sum((S_all - full(X, *pf)) ** 2), np.sum((S_all - red(X, *pr)) ** 2); dfe = len(S_all) - ng - 2
F = (ssr - ssf) / (ssf / dfe)
print(f"Global fit of all {len(S_all)} wells (Kd shared across experiments within each condition): Kd treated / control {pf[ng + 1] / pf[ng]:.2f}; "
      f"extra-sum-of-squares F(1, {dfe}) = {F:.2f}, p = {stats.f.sf(F, 1, dfe):.2g}")