# Per-taxon comparisons of WT and KO three ways: relative abundance with mice as n; CLR with mice as n; CLR averaged
# per cage with cages as n. Benjamini-Hochberg false discovery rate across the 30 taxa. R 4.3.3, base R.
d <- read.csv("counts.csv")                          # one row per mouse: mouse, genotype, cage, taxon_01 ... taxon_30
taxa <- grep("^taxon_", names(d), value = TRUE)
x <- as.matrix(d[, taxa])
rel <- x / rowSums(x)                                # relative abundance
lx <- log(x + 0.5); clr <- lx - rowMeans(lx)         # centred log-ratio (pseudocount 0.5 for zeros)
ko <- d$genotype == "KO"
welch <- function(v, g) sapply(seq_len(ncol(v)), function(j) t.test(v[g, j], v[!g, j])$p.value)
nm <- function(m) if (any(m)) paste(sub("taxon_", "", taxa[m]), collapse = ", ") else "none"
f2 <- function(p) formatC(signif(p, 2), format = "g")

p1 <- welch(rel, ko); down <- colMeans(rel[ko, ]) < colMeans(rel[!ko, ])
cat(sprintf("Relative abundance, mice as n (32 mice), uncorrected p < 0.05: %d of 30 taxa, %d of them lower in KO (%s)\n",
            sum(p1 < 0.05), sum(p1 < 0.05 & down), nm(p1 < 0.05)))
q2 <- p.adjust(welch(clr, ko), "BH")
cat(sprintf("CLR, mice as n, Benjamini-Hochberg q < 0.05: %d of 30 taxa (%s)\n", sum(q2 < 0.05), nm(q2 < 0.05)))
uc <- unique(d$cage); cm <- t(sapply(uc, function(c) colMeans(clr[d$cage == c, , drop = FALSE]))); kc <- sapply(uc, function(c) d$genotype[d$cage == c][1] == "KO")
p3 <- welch(cm, kc); q3 <- p.adjust(p3, "BH")
cat(sprintf("CLR averaged per cage, cages as n (4 vs 4), Benjamini-Hochberg q < 0.05: %d of 30 taxa (%s)\n", sum(q3 < 0.05), nm(q3 < 0.05)))
tt <- t.test(cm[kc, 1], cm[!kc, 1])
cat(sprintf("Taxon 01, KO vs WT relative to the geometric mean of all taxa (cages as n): %.2f-fold (95%% CI %.2f to %.2f), p = %s, q = %s\n",
            exp(mean(cm[kc, 1]) - mean(cm[!kc, 1])), exp(tt$conf.int[1]), exp(tt$conf.int[2]), f2(tt$p.value), f2(q3[1])))