# Sholl analysis, knockout vs wild-type: intersections per radius analysed four ways. R 4.3.3, base R.
d <- read.csv("sholl.csv")                           # one row per neuron and radius: genotype, mouse, neuron, radius_um, intersections
f2 <- function(p) formatC(signif(p, 2), format = "g")
d$cell <- paste(d$mouse, d$neuron, sep = "-")        # neuron IDs repeat across mice, so combine them
Y <- tapply(d$intersections, list(d$cell, d$radius_um), mean)
cg <- tapply(d$genotype, d$cell, `[`, 1)[rownames(Y)]; cm <- tapply(d$mouse, d$cell, `[`, 1)[rownames(Y)]
M <- apply(Y, 2, function(v) tapply(v, cm, mean)); mg <- tapply(d$genotype, d$mouse, `[`, 1)[rownames(M)]   # mouse means per radius
ko <- cg == "knockout"; wt <- cg == "wild-type"; mko <- mg == "knockout"; mwt <- mg == "wild-type"; R <- as.numeric(colnames(Y))
pw <- colMeans(Y[wt, ]); pko <- colMeans(Y[ko, ])
cat(sprintf("Peak of the mean profile: wild-type %.1f intersections at %d um, knockout %.1f at %d um\n", max(pw), R[which.max(pw)], max(pko), R[which.max(pko)]))
p <- sapply(seq_along(R), function(j) t.test(Y[ko, j], Y[wt, j])$p.value)
hits <- if (any(p < 0.05)) paste(paste(R[p < 0.05], "um"), collapse = ", ") else "none"
cat(sprintf("A. Welch t-test at each of %d radii, neurons as n (%d vs %d): %d with p < 0.05 (radii: %s); smallest p = %s\n", length(p), sum(ko), sum(wt), sum(p < 0.05), hits, f2(min(p))))
tot <- rowSums(Y); t <- t.test(tot[ko], tot[wt])
cat(sprintf("B. Total intersections per neuron, neurons as n, Welch t-test: knockout %.1f, wild-type %.1f, p = %s\n", mean(tot[ko]), mean(tot[wt]), f2(t$p.value)))
mt <- rowSums(M); t <- t.test(mt[mko], mt[mwt])
cat(sprintf("C. Total intersections per mouse (mean of its neurons), mice as n (%d vs %d), Welch t-test: difference %.1f (95%% CI %.1f to %.1f), p = %s\n",
            sum(mko), sum(mwt), mean(mt[mko]) - mean(mt[mwt]), t$conf.int[1], t$conf.int[2], f2(t$p.value)))
pm <- sapply(seq_along(R), function(j) t.test(M[mko, j], M[mwt, j])$p.value); q <- p.adjust(pm, method = "BH")
cat(sprintf("D. Welch t-test at each radius on mouse means, Benjamini-Hochberg adjusted: %d radii with q < 0.05; smallest q = %s\n", sum(q < 0.05), f2(min(q))))