0
$\begingroup$

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)
$\endgroup$
3
  • 1
    $\begingroup$ It’s the sum. Please show reproducible code. github.com/mikelove/DESeq2/blob/master/R/helper.R#L186 $\endgroup$
    – user3051
    Jan 17, 2022 at 19:09
  • $\begingroup$ Hmmm... alright. I have updated the question with the code (at the bottom). Let me know if you need more, like the input count matrix (cts) or annotation file (coldata). $\endgroup$
    – Gawain
    Jan 17, 2022 at 19:40
  • 1
    $\begingroup$ I just realized reading the "Plot counts" section of the DESeq2 vignette that the function 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$
    – Gawain
    Jan 17, 2022 at 19:51

1 Answer 1

1
$\begingroup$

The DESeq2 function collapseReplicates sums the counts for the technical replicates.
Here is the code reference: github.com/mikelove/DESeq2/blob/master/R/helper.R#L186

OPs actual confusion was with the DESeq2 function plotCounts which by default normalizes count data and adds a pseudocount of 0.5 for plotting on log2 scale.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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