<p>Comparing read counts from an RNA-seq experiment for two select genes <em>before</em> and <em>after</em> using DESeq2's <code>collapseReplicates()</code> and <code>plotCounts()</code> functions yields interesting results:</p>
<p><strong>Before <code>collapseReplicates()</code> and <code>plotCounts()</code>:</strong></p>
<div class="s-table-container">
<table class="s-table">
<thead>
<tr>
<th style="text-align: left;">Geneid</th>
<th style="text-align: right;">foo1.1</th>
<th style="text-align: right;">foo1.2</th>
<th style="text-align: right;">foo2.1</th>
<th style="text-align: right;">foo2.2</th>
<th style="text-align: right;">bar1.1</th>
<th style="text-align: right;">bar1.2</th>
<th style="text-align: right;">bar2.1</th>
<th style="text-align: right;">bar2.2</th>
<th style="text-align: right;">baz1.1</th>
<th style="text-align: right;">baz1.2</th>
<th style="text-align: right;">baz2.1</th>
<th style="text-align: right;">baz2.2</th>
<th style="text-align: right;">baz3.1</th>
<th style="text-align: right;">baz3.2</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;">WASH7P</td>
<td style="text-align: right;">6</td>
<td style="text-align: right;">5</td>
<td style="text-align: right;">0</td>
<td style="text-align: right;">2</td>
<td style="text-align: right;">1</td>
<td style="text-align: right;">1</td>
<td style="text-align: right;">8</td>
<td style="text-align: right;">5</td>
<td style="text-align: right;">0</td>
<td style="text-align: right;">0</td>
<td style="text-align: right;">0</td>
<td style="text-align: right;">0</td>
<td style="text-align: right;">0</td>
<td style="text-align: right;">0</td>
</tr>
<tr>
<td style="text-align: left;">SOX3</td>
<td style="text-align: right;">1880</td>
<td style="text-align: right;">1861</td>
<td style="text-align: right;">1950</td>
<td style="text-align: right;">2055</td>
<td style="text-align: right;">1189</td>
<td style="text-align: right;">1181</td>
<td style="text-align: right;">2415</td>
<td style="text-align: right;">2482</td>
<td style="text-align: right;">3887</td>
<td style="text-align: right;">3810</td>
<td style="text-align: right;">1851</td>
<td style="text-align: right;">1738</td>
<td style="text-align: right;">3217</td>
<td style="text-align: right;">3406</td>
</tr>
</tbody>
</table>
</div>
<p><strong>After <code>collapseReplicates()</code> and <code>plotCounts()</code>:</strong></p>
<div class="s-table-container">
<table class="s-table">
<thead>
<tr>
<th style="text-align: left;">Geneid</th>
<th style="text-align: right;">foo1</th>
<th style="text-align: right;">foo2</th>
<th style="text-align: right;">bar1</th>
<th style="text-align: right;">bar2</th>
<th style="text-align: right;">baz1</th>
<th style="text-align: right;">baz2</th>
<th style="text-align: right;">baz3</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;">WASH7P</td>
<td style="text-align: right;">9.877191</td>
<td style="text-align: right;">2.279384</td>
<td style="text-align: right;">3.478891</td>
<td style="text-align: right;">11.613875</td>
<td style="text-align: right;">0.500000</td>
<td style="text-align: right;">0.500000</td>
<td style="text-align: right;">0.500000</td>
</tr>
<tr>
<td style="text-align: left;">SOX3</td>
<td style="text-align: right;">3189.598</td>
<td style="text-align: right;">3563.717</td>
<td style="text-align: right;">3530.486</td>
<td style="text-align: right;">4187.011</td>
<td style="text-align: right;">6473.122</td>
<td style="text-align: right;">7991.460</td>
<td style="text-align: right;">5390.721</td>
</tr>
</tbody>
</table>
</div><hr />
<p><strong>Note:</strong> 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).</p>
<hr />
<p>Comparing the tables above, it appears as though DESeq2 is <em>NOT</em> taking the average <em>or</em> sum of columns being collapsed.</p>
<p>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'.</p>
<p>So, how does <code>collapseReplicates()</code> &quot;combine counts into single columns of the count matrix&quot; as is described in the <a href="http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#collapsing-technical-replicates" rel="nofollow noreferrer">DESeq2 vignette</a>?</p>
<hr />
<p>Here is the code to collapse replicates and retrieve the number of read counts for a specific gene (e.g., <em>WASH7P</em>) in a dataframe to be used in a <a href="http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#plot-counts" rel="nofollow noreferrer">count plot</a>, getting the &quot;after&quot; results shown above:</p>
<pre><code>dds &lt;- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)
ddsColl &lt;- collapseReplicates(dds, dds<span class="math-container">$sample, dds$</span>run)
keep &lt;- rowSums(counts(ddsColl)) &gt;= 10
dds &lt;- ddsColl[keep,]
dds<span class="math-container">$condition &lt;- relevel(dds$</span>condition, ref = &quot;baz&quot;)
dds &lt;- DESeq(dds)
countsdf &lt;- plotCounts(dds, gene=&quot;WASH7P&quot;, intgroup=&quot;condition&quot;, returnData=TRUE)
</code></pre>

CC BY-SA 4.0
2022-01-17T19:40:11.127Z
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`).
Gawain
CC BY-SA 4.0
2022-01-17T19:51:03.123Z
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?
Gawain
<rna-seq><deseq2>
How does DESeq2 "collapseReplicates()" function work on read counts data?
# How does DESeq2 "collapseReplicates()" function work on read counts data?

(Asked by: username_0 on 2022-01-17)

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)

**Tags:** rna-seq, deseq2
> **Comments:**
> **username_0**: 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`).  
> **username_0**: 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?  

---

## 0 answers

