# 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.
# R 4.3.3, base R.
d <- read.csv("cells.csv")                           # one row per cell: cell, mouse, condition, gene_001 ... gene_100
genes <- grep("^gene_", names(d), value = TRUE)
x <- as.matrix(d[, genes]); tr <- d$condition == "treated"
ln <- log1p(x / rowSums(x) * 1e4)                    # per-cell normalised log expression
p_cell <- apply(ln, 2, function(v) wilcox.test(v[tr], v[!tr], exact = FALSE)$p.value)
q_cell <- p.adjust(p_cell, "BH")
pb <- aggregate(d[, genes], by = list(mouse = d$mouse, condition = d$condition), FUN = sum)   # one profile per mouse
pb <- pb[order(pb$mouse), ]
a <- as.matrix(pb[, genes]); cpm <- log2(a / rowSums(a) * 1e6 + 1); tm <- pb$condition == "treated"
p_pb <- apply(cpm, 2, function(v) t.test(v[tm], v[!tm])$p.value); q_pb <- p.adjust(p_pb, "BH")
summ <- function(q) { s <- q < 0.05; sprintf("%d genes (%d of the 5 truly changed, %d of the 95 unchanged)", sum(s), sum(s[1:5]), sum(s[-(1:5)])) }
cat(sprintf("Cells as replicates (%d vs %d cells), Wilcoxon, BH q < 0.05: %s\n", sum(tr), sum(!tr), summ(q_cell)))
cat(sprintf("Pseudobulk (%d vs %d mice), Welch's t-test on log2 CPM, BH q < 0.05: %s\n", sum(tm), sum(!tm), summ(q_pb)))
lfc <- colMeans(cpm[tm, ]) - colMeans(cpm[!tm, ])
cat(sprintf("Pseudobulk log2 fold change, genes 001-005: %s; largest absolute value among the 95 unchanged genes %.2f\n",
            paste(sprintf("%.2f", lfc[1:5]), collapse = ", "), max(abs(lfc[-(1:5)]))))
top <- order(p_pb)[1:5]
cat(sprintf("Pseudobulk, five smallest p-values: %s\n", paste(sprintf("%s (p = %.3f)", sub("gene_", "", genes[top]), p_pb[top]), collapse = ", ")))