# 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. R 4.3.3, base R.
d <- read.csv("intensities.csv"); truth <- read.csv("truth.csv")$true_log2fc
x <- as.matrix(d[, -1])                              # proteins x samples, log2 intensity, NA = missing
ctl <- x[, 1:5]; trt <- x[, 6:10]
cat(sprintf("Missing values: %.1f%% of all values; proteins with any missing value: %d of %d\n", 100 * mean(is.na(x)), sum(rowSums(is.na(x)) > 0), nrow(x)))
report <- function(lab, p, keep) {
  q <- rep(NA, nrow(x)); q[keep] <- p.adjust(p, "BH"); called <- keep & !is.na(q) & q < 0.05
  fdp <- sum(called & truth == 0) / max(sum(called), 1)
  cat(sprintf("%s: %d proteins tested; %d called (q < 0.05), %d of them unchanged (false discovery proportion %.1f%%); %d of the 100 changed proteins found\n",
              lab, sum(keep), sum(called), sum(called & truth == 0), 100 * fdp, sum(called & truth != 0)))
}
welch <- function(m) sapply(seq_len(nrow(m)), function(i) t.test(m[i, 6:10], m[i, 1:5])$p.value)
cmin <- apply(x, 2, min, na.rm = TRUE); xm <- x; for (j in seq_len(ncol(x))) xm[is.na(x[, j]), j] <- cmin[j]
report("Impute each sample's minimum", welch(xm), rep(TRUE, nrow(x)))
shift <- apply(x, 2, function(v) mean(v, na.rm = TRUE) - 1.8 * sd(v, na.rm = TRUE)); xs <- x
for (j in seq_len(ncol(x))) xs[is.na(x[, j]), j] <- shift[j]
report("Impute a down-shifted value (sample mean - 1.8 SD)", welch(xs), rep(TRUE, nrow(x)))
keep <- rowSums(!is.na(ctl)) >= 3 & rowSums(!is.na(trt)) >= 3
pk <- sapply(which(keep), function(i) t.test(trt[i, !is.na(trt[i, ])], ctl[i, !is.na(ctl[i, ])])$p.value)
report("Keep proteins with at least 3 valid values per group, no imputation", pk, keep)