# Immunofluorescence intensity, knockdown vs control: corrected total cell fluorescence (CTCF) per cell, analysed with
# cells, images or experiments as n, and without background subtraction. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats
d = pd.read_csv("cells.csv")                         # one row per cell: experiment, condition, image, cell, area_px, int_den, bg_mean (per image)
d["ctcf"] = d.int_den - d.area_px * d.bg_mean        # CTCF = integrated density - area x mean background
d["log_ctcf"] = np.log(d.ctcf); d["log_raw"] = np.log(d.int_den)
print(f"Cells: {len(d)}, CTCF median {d.ctcf.median():,.0f}, smallest {d.ctcf.min():,.0f}; background is {100 * (d.area_px * d.bg_mean).sum() / d.int_den.sum():.0f}% of the raw integrated density")
img = d.groupby(["experiment", "condition", "image"])[["log_ctcf", "log_raw"]].mean().reset_index()
ex = img.groupby(["experiment", "condition"])[["log_ctcf", "log_raw"]].mean().unstack()
kd, ct = d.condition == "knockdown", d.condition == "control"
t = stats.ttest_ind(d.log_ctcf[kd], d.log_ctcf[ct], equal_var=False)
print(f"A. Cells as n ({kd.sum()} vs {ct.sum()}), Welch t-test on log CTCF: knockdown / control {np.exp(d.log_ctcf[kd].mean() - d.log_ctcf[ct].mean()):.2f}, p = {t.pvalue:.2g}")
ik, ic = img.condition == "knockdown", img.condition == "control"
t = stats.ttest_ind(img.log_ctcf[ik], img.log_ctcf[ic], equal_var=False)
print(f"B. Images as n ({ik.sum()} vs {ic.sum()}), Welch t-test on mean log CTCF: knockdown / control {np.exp(img.log_ctcf[ik].mean() - img.log_ctcf[ic].mean()):.2f}, p = {t.pvalue:.2g}")
for v, lab in (("log_ctcf", "C. Experiments as n, CTCF"), ("log_raw", "D. Experiments as n, raw integrated density (no background subtraction)")):
    r = ex[(v, "knockdown")] - ex[(v, "control")]; t = stats.ttest_1samp(r, 0); ci = t.confidence_interval()
    print(f"{lab}, log ratio per experiment (n = {len(r)}): knockdown / control {np.exp(r.mean()):.2f} (95% CI {np.exp(ci.low):.2f} to {np.exp(ci.high):.2f}), p = {t.pvalue:.2g}")