# Summarize expression in the gene count matrix.
# Samples are averaged within each group from the design file.
# A gene is not expressed when that group mean is 0.
# Bins are mean raw counts: 0, <1, 1-9, 10-99, 100-999, 1000+.

group_means <- function(counts, design) {
  groups <- unique(design$group)
  means <- vapply(groups, function(g) {
    samples <- design$sample[design$group == g]
    rowMeans(counts[, samples, drop = FALSE])
  }, numeric(nrow(counts)))
  colnames(means) <- groups
  means
}

expr_bins <- function(x) {
  labs <- c("0", "<1", "1-9", "10-99", "100-999", "1000+")
  idx <- rep(1L, length(x))
  pos <- x > 0
  idx[pos] <- findInterval(x[pos], c(1, 10, 100, 1000)) + 2L
  setNames(tabulate(idx, nbins = length(labs)), labs)
}

print_summary <- function(means) {
  groups <- colnames(means)
  bins <- sapply(groups, function(g) expr_bins(means[, g]))
  silent <- colSums(means == 0)
  expressed <- colSums(means > 0)

  line <- function(label, vals) {
    cat(sprintf("%-16s", label), paste(sprintf("%8s", vals), collapse = ""), "\n", sep = "")
  }

  line("", groups)
  line("not expressed", silent)
  line("expressed", expressed)
  cat("\n")
  line("level", groups)
  for (lvl in rownames(bins)) {
    line(lvl, bins[lvl, ])
  }
}

design <- read.csv("design.csv", stringsAsFactors = FALSE)
counts <- read.csv("results/gene_counts.csv", check.names = FALSE)
print_summary(group_means(counts, design))