How to input data for DESeq2 from individual HTSeq count?

Question

I am comparing the gene expression of 2 bacteria under 1 condition. I have now the count tables for 3 tech. replicates for each bacteria.

Bacteria1_1.count
Bacteria1_2.count
Bacteria1_3.count

...same for the other bacteria.

These files look like this:

gene1 10000
gene2 500
gene3 0
gene4 5000

I want to use DESeq2 for differential gene expression analysis. But I cannot figure out how to properly execute the DESeqDataSetFromHTSeqCount() command with this type of data.

Is there another intermediate step to add?

Answer 1

In case still needed, here is a piece of code that I have used for the same purpose (includes a lot of sanity checks). This is highly likely an overkill as each count file would probably have the same rows (as a result of being generated from the same workflow) and thus can easily be merged column-wise.

count_file_dir <- "../input/count_files/"
count_file_names <- list.files(count_file_dir)
count_files <- paste0(count_file_dir, count_file_names)

# read and merge count files, while removing last 5 lines that correspond to summary stats
count_table_list <- mclapply(count_files, function(x) {fread(x)[1:(.N-5)]})
names(count_table_list) <- sub(".count", "", count_file_names)

# check the dimensions of count files
t(sapply(count_table_list, dim))

# check if the count files have the same gene ordering
# compare gene names of the first file with all
gene_names <- sapply(count_table_list, function(x) x$V1)
as.data.frame(apply(gene_names, 2, function(x) all(x == gene_names[,1])))

# add file name as a unique column name for gene counts
for(i in names(count_table_list)) {
  names(count_table_list[[i]])[2] <- i
}

# merge read counts and convert to a matrix
count_matrix <- Reduce(merge, count_table_list)
count_matrix_rownames <- count_matrix$V1
count_matrix$V1 <- NULL
count_matrix <- as.matrix(count_matrix)
rownames(count_matrix) <- count_matrix_rownames
colnames(count_matrix) <- sub(".count", "", colnames(count_matrix))

dim(count_matrix)

dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                              colData = ...,
                              design = ...)