# Mendelian genotype ratios from heterozygous x heterozygous crosses: chi-square goodness of fit, exact tests and
# knockout viability relative to wild-type littermates with an exact 95% CI. R 4.3.3, base R.
d <- read.csv("pups.csv")                            # one row per genotyped pup: litter, pup, sex, genotype
f2 <- function(p) formatC(signif(p, 2), format = "g")
n_wt <- sum(d$genotype == "wild-type"); n_het <- sum(d$genotype == "heterozygous"); n_ko <- sum(d$genotype == "knockout"); n <- n_wt + n_het + n_ko
obs <- c(n_wt, n_het, n_ko); p0 <- c(0.25, 0.5, 0.25); ex <- n * p0
cat(sprintf("%d litters, %d pups genotyped at weaning\n", length(unique(d$litter)), n))
cat(sprintf("observed wild-type %d, heterozygous %d, knockout %d; expected %.2f, %.2f, %.2f\n", n_wt, n_het, n_ko, ex[1], ex[2], ex[3]))
cs <- chisq.test(obs, p = p0)
cat(sprintf("A. Chi-square goodness of fit to 1:2:1: chi2(2) = %.2f, p = %s\n", cs$statistic, f2(cs$p.value)))
multinom_exact <- function(x, p) {                   # sum of the probabilities of all outcomes no more likely than x
  n <- sum(x); px <- dmultinom(x, n, p); tot <- 0
  for (a in 0:n) for (b in 0:(n - a)) { pr <- dmultinom(c(a, b, n - a - b), n, p); if (pr <= px * (1 + 1e-7)) tot <- tot + pr }
  min(tot, 1) }
cat(sprintf("B. Exact multinomial test of 1:2:1: p = %s\n", f2(multinom_exact(obs, p0))))
b <- binom.test(n_ko, n, 0.25)
cat(sprintf("C. Exact binomial test, knockouts against 1/4 of pups: %d of %d (%.1f%%), p = %s\n", n_ko, n, 100 * n_ko / n, f2(b$p.value)))
r <- binom.test(n_ko, n_ko + n_wt, 0.5); v <- function(q) q / (1 - q)   # knockout share of homozygotes -> knockout / wild-type ratio
cat(sprintf("D. Knockout viability relative to wild-type littermates: %.2f (exact 95%% CI %.2f to %.2f), p = %s\n", v(n_ko / (n_ko + n_wt)), v(r$conf.int[1]), v(r$conf.int[2]), f2(r$p.value)))
fe <- sum(d$sex == "female"); s <- binom.test(fe, n, 0.5)
cat(sprintf("Sex ratio: %d females of %d, exact binomial p = %s\n", fe, n, f2(s$p.value)))