2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ What errors are you seeing? How are you trying to do this? Please edit your question to show us the commands you tried. $\endgroup$ – terdon Oct 13 '18 at 13:08
  • $\begingroup$ Probably you should join all the count tables into a single one (If they share the same rows) $\endgroup$ – llrs Oct 16 '18 at 7:16
1
$\begingroup$

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 = ...)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.