# Reading a differential expression table: raw p, BH-adjusted p, a fold-change filter applied after testing, and a test
# against a fold-change threshold (TREAT, McCarthy & Smyth 2009). Student's t-test per gene on log2 expression, 4 vs 4.
# Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats

d = pd.read_csv("log2_expression.csv"); truth = pd.read_csv("truth.csv").true_log2fc.to_numpy()
ctl = d.filter(like="control_").to_numpy(); trt = d.filter(like="treated_").to_numpy()
n1, n2 = ctl.shape[1], trt.shape[1]; df = n1 + n2 - 2
lfc = trt.mean(1) - ctl.mean(1)
s2 = ((n1 - 1) * ctl.var(1, ddof=1) + (n2 - 1) * trt.var(1, ddof=1)) / df; se = np.sqrt(s2 * (1 / n1 + 1 / n2))
p = 2 * stats.t.sf(np.abs(lfc) / se, df); q = stats.false_discovery_control(p)
tau = 1.0                                             # threshold: a 2-fold change
p_treat = stats.t.sf((np.abs(lfc) - tau) / se, df) + stats.t.sf((np.abs(lfc) + tau) / se, df)
q_treat = stats.false_discovery_control(p_treat)
def show(lab, called):
    print(f"{lab}: {called.sum()} genes; {np.sum(called & (truth == 0))} truly unchanged; "
          f"{np.sum(called & (np.abs(truth) < tau))} with a true change below 2-fold (including unchanged)")
show("Raw p < 0.05", p < 0.05)
show("BH-adjusted p (padj) < 0.05", q < 0.05)
show("padj < 0.05 and |log2 fold change| > 1, filtered after testing", (q < 0.05) & (np.abs(lfc) > tau))
show("Test against |log2 fold change| > 1 (TREAT), adjusted p < 0.05", q_treat < 0.05)
print(f"Truly changed genes: {np.sum(truth != 0)}, of which {np.sum(np.abs(truth) > tau)} by more than 2-fold")