Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-01-17T19:40:11.127Z
Comments.Text: 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`).
Comments.UserDisplayName: Gawain
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-01-17T19:51:03.123Z
Comments.Text: 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?
Comments.UserDisplayName: Gawain
Tags: <rna-seq><deseq2>
Title: How does DESeq2 "collapseReplicates()" function work on read counts data?
ThreadText: # 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](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#collapsing-technical-replicates)?
  
  ---
  
  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](http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#plot-counts), 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
  
  