# Three ways to handle missing values before a per-protein Welch's t-test with Benjamini-Hochberg correction:
# impute each sample's minimum, impute a down-shifted value (sample mean - 1.8 SD), or keep proteins with at least
# 3 valid values per group and test the observed values only. Python 3.13, NumPy, pandas, SciPy.
import numpy as np, pandas as pd
from scipy import stats

d = pd.read_csv("intensities.csv"); truth = pd.read_csv("truth.csv").true_log2fc.to_numpy()
x = d.drop(columns="protein").to_numpy(float)       # proteins x samples, log2 intensity, NaN = missing
ctl, trt = x[:, :5], x[:, 5:]
print(f"Missing values: {np.isnan(x).mean():.1%} of all values; proteins with any missing value: {np.isnan(x).any(1).sum()} of {len(x)}")
def report(lab, p, keep):
    q = np.full(len(x), np.nan); q[keep] = stats.false_discovery_control(p)
    called = keep & (q < 0.05)
    fdp = np.sum(called & (truth == 0)) / max(called.sum(), 1)
    print(f"{lab}: {keep.sum()} proteins tested; {called.sum()} called (q < 0.05), {np.sum(called & (truth == 0))} of them unchanged "
          f"(false discovery proportion {fdp:.1%}); {np.sum(called & (truth != 0))} of the 100 changed proteins found")
def welch(a, b): return stats.ttest_ind(a, b, axis=1, equal_var=False).pvalue
cmin = np.nanmin(x, axis=0); xm = np.where(np.isnan(x), cmin, x)
report("Impute each sample's minimum", welch(xm[:, 5:], xm[:, :5]), np.ones(len(x), bool))
shift = np.nanmean(x, axis=0) - 1.8 * np.nanstd(x, axis=0, ddof=1); xs = np.where(np.isnan(x), shift, x)
report("Impute a down-shifted value (sample mean - 1.8 SD)", welch(xs[:, 5:], xs[:, :5]), np.ones(len(x), bool))
keep = (np.sum(~np.isnan(ctl), 1) >= 3) & (np.sum(~np.isnan(trt), 1) >= 3)
pk = np.array([stats.ttest_ind(trt[i][~np.isnan(trt[i])], ctl[i][~np.isnan(ctl[i])], equal_var=False).pvalue for i in np.where(keep)[0]])
report("Keep proteins with at least 3 valid values per group, no imputation", pk, keep)