combining quantification (featureCounts) result files into a single dataset

Below is the function I tend to use to read in multiple featureCounts outputs (one per sample):

DESeqDataSetFromFeatureCounts <- function (sampleTable, directory = ".", design, ignoreRank = FALSE, ...) 
{
    if (missing(design)) 
        stop("design is missing")
    l <- lapply(as.character(sampleTable[, 2]), function(fn) read.table(file.path(directory, fn), skip=2))
    if (!all(sapply(l, function(a) all(a$V1 == l[[1]]$V1)))) 
        stop("Gene IDs (first column) differ between files.")
    tbl <- sapply(l, function(a) a$V7)
    colnames(tbl) <- sampleTable[, 1]
    rownames(tbl) <- l[[1]]$V1
    rownames(sampleTable) <- sampleTable[, 1]
    dds <- DESeqDataSetFromMatrix(countData = tbl, colData = sampleTable[, 
        -(1:2), drop = FALSE], design = design, ignoreRank, ...)
    return(dds)
}

You’ll note that it’s essentially identical to DESeqDataSetFromHTSeqCount() or whatever that’s called, just tweaked slightly to work with featureCounts output. You’ll need a sampleTable, just as with the standard functions built into DESeq2.