# Colocalisation: Pearson's and Manders' coefficients, and comparing conditions with the right n. R 4.3.3, base R.
im <- read.csv("images.csv")
cat("Coefficients for the three example images (threshold = background mean + 3 SD = 16 in both channels):\n")
for (name in sort(unique(im$image))) {
  g <- im[im$image == name, ]; a <- g$ch1; b <- g$ch2; th <- 16
  cat(sprintf("  image %s: Pearson %.3f; Manders M1 %.3f, M2 %.3f\n", name, cor(a, b), sum(a[a > th & b > th]) / sum(a[a > th]), sum(b[a > th & b > th]) / sum(b[b > th])))
}
d <- read.csv("cells.csv"); d$z <- atanh(d$pearson)
t1 <- t.test(d$z[d$condition == "treated"], d$z[d$condition == "control"])
cat(sprintf("Per-cell Welch's t-test (150 cells treated as independent): t = %.2f, p = %.4f\n", t1$statistic, t1$p.value))
m <- tapply(d$z, list(d$experiment, d$condition), mean)
diff <- m[, "treated"] - m[, "control"]; t2 <- t.test(m[, "treated"], m[, "control"], paired = TRUE)
cat("Mean Pearson per experiment (back-transformed): ", paste(sprintf("experiment %s: control %.3f, treated %.3f", rownames(m), tanh(m[, "control"]), tanh(m[, "treated"])), collapse = "; "), "\n", sep = "")
hw <- qt(0.975, length(diff) - 1) * sd(diff) / sqrt(length(diff))
cat(sprintf("Paired t-test on experiment means (n = 3 experiments): difference in z %.3f (95%% CI %.3f to %.3f), t = %.2f, p = %.4f\n",
            mean(diff), mean(diff) - hw, mean(diff) + hw, t2$statistic, t2$p.value))