# Growth rate and doubling time from cell counts; n = independent experiments. R 4.3.3, base R.
d <- read.csv("counts.csv"); d$logc <- log(d$cells)
g <- aggregate(logc ~ experiment + line + hour, data = d, mean)
rate <- function(x, hours) { x <- x[x$hour %in% hours, ]; unname(coef(lm(logc ~ hour, data = x))[2]) }
keys <- unique(g[, c("experiment", "line")]); keys <- keys[order(keys$experiment, keys$line), ]
cat("Growth rate per hour (doubling time, h), fitted over 0-72 h:\n")
r <- data.frame()
for (i in seq_len(nrow(keys))) {
  x <- g[g$experiment == keys$experiment[i] & g$line == keys$line[i], ]
  re <- rate(x, c(0, 24, 48, 72)); ra <- rate(x, c(0, 24, 48, 72, 96))
  r <- rbind(r, data.frame(experiment = keys$experiment[i], line = keys$line[i], r_exp = re))
  cat(sprintf("  experiment %d %s: %.4f (%.1f h); including 96 h: %.4f (%.1f h)\n", keys$experiment[i], keys$line[i], re, log(2) / re, ra, log(2) / ra))
}
w <- r$r_exp[r$line == "WT"]; k <- r$r_exp[r$line == "KO"]; diff <- k - w
t1 <- t.test(k, w, paired = TRUE); hw <- qt(0.975, length(diff) - 1) * sd(diff) / sqrt(length(diff))
cat(sprintf("Mean doubling time: WT %.1f h, KO %.1f h\n", log(2) / mean(w), log(2) / mean(k)))
cat(sprintf("KO - WT growth rate, paired over 3 experiments: %.4f per hour (95%% CI %.4f to %.4f), p = %.4f\n", mean(diff), mean(diff) - hw, mean(diff) + hw, t1$p.value))
c96 <- d[d$hour == 96, ]
cat(sprintf("For comparison, Welch's t-test on 96-h counts with the 9 wells per line as n: p = %s\n",
            formatC(signif(t.test(c96$cells[c96$line == "KO"], c96$cells[c96$line == "WT"])$p.value, 2), format = "g")))