# Organoid size, drug vs control: the same organoids analysed with organoids, wells or batches as n.
# Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats
d = pd.read_csv("organoids.csv")                     # one row per organoid: batch, condition, well, area_um2
d["log_area"] = np.log(d.area_um2)
w = d.groupby(["batch", "condition", "well"]).log_area.agg(["mean", "size"]).reset_index()      # mean log area per well
b = w.groupby(["batch", "condition"])["mean"].mean().unstack()                                   # mean of the wells per batch
for k, r in b.iterrows():
    print(f"Batch {k}: geometric mean area control {np.exp(r.control):,.0f} um2, drug {np.exp(r.drug):,.0f} um2 (ratio {np.exp(r.drug - r.control):.2f})")
print(f"Organoids measured: {len(d)} ({w['size'].min()} to {w['size'].max()} per well)")
def ratio_test(x, y, label):
    t = stats.ttest_ind(x, y, equal_var=False)
    print(f"{label} ({len(x)} vs {len(y)}), Welch t-test on log area: drug / control {np.exp(x.mean() - y.mean()):.2f}, p = {t.pvalue:.2g}")
ratio_test(d.log_area[d.condition == "drug"], d.log_area[d.condition == "control"], "A. Organoids as n")
ratio_test(w["mean"][w.condition == "drug"], w["mean"][w.condition == "control"], "B. Wells as n")
r = b.drug - b.control; t = stats.ttest_1samp(r, 0); ci = t.confidence_interval()
print(f"C. Batches as n, log ratio per batch (paired, n = {len(r)}): drug / control {np.exp(r.mean()):.2f} (95% CI {np.exp(ci.low):.2f} to {np.exp(ci.high):.2f}), p = {t.pvalue:.2g}")