# Combine per-sample featureCounts tables into one count matrix.
# Column names and order come from the sample column of the design file.
# The name column is the common gene name from the GTF used for counting.

suppressPackageStartupMessages(library(rtracklayer))

read_counts <- function(sample, dir = "counts") {
  path <- file.path(dir, paste0(sample, ".txt"))
  tab <- read.delim(path, comment.char = "#", check.names = FALSE)
  setNames(tab[[ncol(tab)]], tab$Geneid)
}

read_names <- function(gtf = "gtf/chr22.gtf") {
  ann <- import(gtf)
  genes <- ann[ann$type == "gene"]
  setNames(genes$gene_name, genes$gene_id)
}

combine_counts <- function(design = "design.csv", dir = "counts", gtf = "gtf/chr22.gtf") {
  samples <- read.csv(design, stringsAsFactors = FALSE)$sample
  cols <- lapply(samples, read_counts, dir = dir)

  genes <- names(cols[[1]])
  mat <- vapply(cols, function(x) x[genes], numeric(length(genes)))
  colnames(mat) <- samples

  data.frame(
    gene = genes,
    name = unname(read_names(gtf)[genes]),
    mat,
    check.names = FALSE,
    stringsAsFactors = FALSE
  )
}

counts <- combine_counts()
dir.create("results", showWarnings = FALSE)
write.csv(counts, "results/gene_counts.csv", row.names = FALSE, quote = FALSE)