# Per-taxon comparisons of WT and KO three ways: relative abundance with mice as n; CLR with mice as n; CLR averaged
# per cage with cages as n. Benjamini-Hochberg false discovery rate across the 30 taxa. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats

d = pd.read_csv("counts.csv")                        # one row per mouse: mouse, genotype, cage, taxon_01 ... taxon_30
taxa = [c for c in d.columns if c.startswith("taxon_")]
x = d[taxa].to_numpy(float)
rel = x / x.sum(axis=1, keepdims=True)               # relative abundance
lx = np.log(x + 0.5); clr = lx - lx.mean(axis=1, keepdims=True)   # centred log-ratio (pseudocount 0.5 for zeros)
ko = (d.genotype == "KO").to_numpy()

def welch(v, g):                                     # per-taxon Welch's t-test, KO vs WT
    return np.array([stats.ttest_ind(v[g, j], v[~g, j], equal_var=False).pvalue for j in range(v.shape[1])])
def names(mask): return ", ".join(t.replace("taxon_", "") for t, m in zip(taxa, mask) if m) or "none"

p1 = welch(rel, ko)
down = (rel[ko].mean(0) < rel[~ko].mean(0))
print(f"Relative abundance, mice as n (32 mice), uncorrected p < 0.05: {np.sum(p1 < 0.05)} of 30 taxa, {np.sum((p1 < 0.05) & down)} of them lower in KO ({names(p1 < 0.05)})")
q2 = stats.false_discovery_control(welch(clr, ko))
print(f"CLR, mice as n, Benjamini-Hochberg q < 0.05: {np.sum(q2 < 0.05)} of 30 taxa ({names(q2 < 0.05)})")
cages = d.cage.to_numpy(); uc = list(dict.fromkeys(cages))
cm = np.array([clr[cages == c].mean(axis=0) for c in uc]); kc = np.array([d.genotype[cages == c].iloc[0] == "KO" for c in uc])
p3 = welch(cm, kc); q3 = stats.false_discovery_control(p3)
print(f"CLR averaged per cage, cages as n (4 vs 4), Benjamini-Hochberg q < 0.05: {np.sum(q3 < 0.05)} of 30 taxa ({names(q3 < 0.05)})")
t = stats.ttest_ind(cm[kc, 0], cm[~kc, 0], equal_var=False); ci = t.confidence_interval()
print(f"Taxon 01, KO vs WT relative to the geometric mean of all taxa (cages as n): {np.exp(cm[kc, 0].mean() - cm[~kc, 0].mean()):.2f}-fold "
      f"(95% CI {np.exp(ci.low):.2f} to {np.exp(ci.high):.2f}), p = {t.pvalue:.2g}, q = {q3[0]:.2g}")