Comparing read counts from an RNA-seq experiment for two select genes before and after using DESeq2's
plotCounts() functions yields interesting results:
Note: In the above comparison, there are two (2) technical replicates (1.1 + 1.2, 2.1 + 2.2, etc.) for each biological replicate (foo1, foo2, etc.) for each of three (3) conditions (foo, bar, and baz).
Comparing the tables above, it appears as though DESeq2 is NOT taking the average or sum of columns being collapsed.
It is also curious - and mildly concerning - that some very low expression (i.e., 0.5 counts) is reported for genes in the matrix of collapsed replicates when, in the original count matrix, zero (0) reads were counted as 'mapped to that gene'.
So, how does
collapseReplicates() "combine counts into single columns of the count matrix" as is described in the DESeq2 vignette?
Here is the code to collapse replicates and retrieve the number of read counts for a specific gene (e.g., WASH7P) in a dataframe to be used in a count plot, getting the "after" results shown above:
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition) ddsColl <- collapseReplicates(dds, dds$sample, dds$run) keep <- rowSums(counts(ddsColl)) >= 10 dds <- ddsColl[keep,] dds$condition <- relevel(dds$condition, ref = "baz") dds <- DESeq(dds) countsdf <- plotCounts(dds, gene="WASH7P", intgroup="condition", returnData=TRUE)