# Treated vs control for one cell type: cells as replicates (Wilcoxon rank-sum per gene, like Seurat's default) versus
# pseudobulk (counts summed per mouse, log2 counts per million, Welch's t-test per gene). Benjamini-Hochberg across genes.
# 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: cell, mouse, condition, gene_001 ... gene_100
genes = [c for c in d.columns if c.startswith("gene_")]
x = d[genes].to_numpy(float); tr = (d.condition == "treated").to_numpy()
lib = x.sum(1, keepdims=True); ln = np.log1p(x / lib * 1e4)      # per-cell normalised log expression
p_cell = np.array([stats.mannwhitneyu(ln[tr, j], ln[~tr, j], alternative="two-sided").pvalue for j in range(len(genes))])
q_cell = stats.false_discovery_control(p_cell)
pb = d.groupby(["mouse", "condition"])[genes].sum().reset_index()   # pseudobulk: one profile per mouse
cpm = np.log2(pb[genes].to_numpy(float) / pb[genes].to_numpy(float).sum(1, keepdims=True) * 1e6 + 1)
tm = (pb.condition == "treated").to_numpy()
p_pb = stats.ttest_ind(cpm[tm], cpm[~tm], equal_var=False).pvalue
q_pb = stats.false_discovery_control(p_pb)
def summary(q):
    s = q < 0.05; return f"{s.sum()} genes ({s[:5].sum()} of the 5 truly changed, {s[5:].sum()} of the 95 unchanged)"
print(f"Cells as replicates ({tr.sum()} vs {(~tr).sum()} cells), Wilcoxon, BH q < 0.05: {summary(q_cell)}")
print(f"Pseudobulk ({tm.sum()} vs {(~tm).sum()} mice), Welch's t-test on log2 CPM, BH q < 0.05: {summary(q_pb)}")
lfc = cpm[tm].mean(0) - cpm[~tm].mean(0)
print("Pseudobulk log2 fold change, genes 001-005: " + ", ".join(f"{v:.2f}" for v in lfc[:5]) +
      f"; largest absolute value among the 95 unchanged genes {np.abs(lfc[5:]).max():.2f}")
top = np.argsort(p_pb, kind="stable")[:5]
print("Pseudobulk, five smallest p-values: " + ", ".join(f"{genes[j].replace('gene_', '')} (p = {p_pb[j]:.3f})" for j in top))