Comparing read counts from an RNA-seq experiment for two select genes before and after using DESeq2's collapseReplicates()
and plotCounts()
functions yields interesting results:
Before collapseReplicates()
and plotCounts()
:
Geneid | foo1.1 | foo1.2 | foo2.1 | foo2.2 | bar1.1 | bar1.2 | bar2.1 | bar2.2 | baz1.1 | baz1.2 | baz2.1 | baz2.2 | baz3.1 | baz3.2 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
WASH7P | 6 | 5 | 0 | 2 | 1 | 1 | 8 | 5 | 0 | 0 | 0 | 0 | 0 | 0 |
SOX3 | 1880 | 1861 | 1950 | 2055 | 1189 | 1181 | 2415 | 2482 | 3887 | 3810 | 1851 | 1738 | 3217 | 3406 |
After collapseReplicates()
and plotCounts()
:
Geneid | foo1 | foo2 | bar1 | bar2 | baz1 | baz2 | baz3 |
---|---|---|---|---|---|---|---|
WASH7P | 9.877191 | 2.279384 | 3.478891 | 11.613875 | 0.500000 | 0.500000 | 0.500000 |
SOX3 | 3189.598 | 3563.717 | 3530.486 | 4187.011 | 6473.122 | 7991.460 | 5390.721 |
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)
cts
) or annotation file (coldata
). $\endgroup$plotCounts
"normalizes counts by the estimated size factors (or normalization factors if these were used) and adds a pseudocount of 1/2 to allow for log scale plotting", which would explain the "0.5" I observe for the genes with "0" reads mapped to those genes. This is good to know, but I now wonder if this data can be used to accurately illustrate the expression of these genes in a plot? $\endgroup$