<p>Absolutely, look at Wilcox-rank, not purely because of the computational overhead, but the parametric assumption that DESeq2 uses  <a href="https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02648-4" rel="nofollow noreferrer">are not correct</a> on large data sets. This might describe your data sets too:</p>
<blockquote>
<p><strong>Abstract</strong>
DESeq2 and edgeR, have unexpectedly high false discovery rates.</p>
<p>Expanding the analysis to limma-voom, NOISeq, dearseq, and <strong>Wilcoxon rank-sum test</strong>, we found that [false discovery rates] FDR control is often failed <strong>except for the Wilcoxon rank-sum test</strong></p>
<p>Particularly, the actual <strong>FDRs of DESeq2</strong> and edgeR sometimes exceed <strong>20%</strong> when the target FDR is 5%. Based on these results, for population-level RNA-seq studies with large sample sizes, <strong>we recommend the Wilcoxon rank-sum test</strong>.</p>
<p><strong>Results and Discussion</strong>
DESeq2 and edgeR had only an 8% overlap in the DEGs they identified</p>
</blockquote>
<hr />

CC BY-SA 4.0
2022-12-20T15:00:39.307Z
I feel it's unfair to write that without mentioning the DESeq author [response](https://twitter.com/mikelove/status/1513468614585765889): in short, there are technical confounders that need to be controlled for. Wilcoxon gets fewer False Positives because of its lower power (while DESeq and edgeR are designed to maximize power), but the right thing to do is really to correct for confounding, whatever test you do next. That being said for this question Wilcoxon is probably reasonable (though a correction with RUV/SVA/... would be appropriate).
Alexlok
CC BY-SA 4.0
2022-12-25T08:42:11.837Z
It is not a matter of parametric assumptions. The linked paper, as said already did not address the fact that in larger cohorts you often find samples that form clusters by some unwanted confounder nested with the variable you're testing. That is not unusual especially in human samples where you have all sorts of confounders, be it eating habits, drug consumption, dietary status etc that one often does not know of. Once this was addressed (RUVseq by the DESeq2 author) the main claim of that paper collapsed. In any case, at larger n I would use linear approaches such as limma-voom.
ATpoint
CC BY-SA 4.0
2022-12-25T15:12:21.587Z
Hi @ATpoint thanks its a good discussion, albeit maybe not this site.
M__
CC BY-SA 4.0
2022-12-25T17:40:20.453Z
@M__, I think these comments are necessary to put answers into context. The paper you reference is basically an advertisement of the non-parametric tool the authors created. It's in the nature of these sort of papers that advertised tool is best and all others are not as good. That goes for any benchmark paper written by method developers, be it aligners, variant callers or stats frameworks. I am not judging this, it is just an observation. Like any of these papers, results should be read with a grain of salt and I think it is important to notice that this was not an unbiased benchmark.
ATpoint
CC BY-SA 4.0
2022-12-25T18:21:25.413Z
@ATpoint we can have this discussion via the site  that you are famous on. There been an objection raised on this site and we must be mindful of that. I will raise this as a question on *your site* not quite now because it's Christmas Day (quite literally), but shortly. I agree 100% is an important topic and I'll explain my interest in the subject (I don't really do eukaryotes BTW).
M__
CC BY-SA 4.0
2022-12-20T00:22:38.917Z
2022-12-24T07:53:34.673Z
2022-12-24T07:53:34.673Z
M__
<p>I'll preface this by saying that I don't think DESeq2 is the right tool to use for ATAC-Seq data. My own study of ATAC-Seq patterns [admittedly only a couple of runs that were our first exploration of the technique] indicates that it has a fractal expression which is difficult to properly summarise at any scale, especially at a bulk sample level. I found <a href="https://twitter.com/ryanjames722/status/1545261940406116352" rel="nofollow noreferrer">this poster</a> from a quick search on Twitter, which recommends <a href="https://github.com/Zhang-lab/ATAC-seq_QC_analysis" rel="nofollow noreferrer">AIAP</a> for ATAC-Seq analysis (but don't read too much into that; I have no experience with doing ATAC-Seq analysis). What I <em>do</em> have a lot of experience with is DESeq2.</p>
<p>DESeq2 has certain assumptions about the expression patterns of genes which are built into its differential expression model, and the further you get from gene-level analysis, the less appropriate DESeq2 is for data analysis (for example, <a href="https://doi.org/10.1093/nar/gkz622" rel="nofollow noreferrer">Swish</a> is a better approach for isoform level and allelic analysis).</p>
<p>However, if you're absolutely set on using DESeq2...  the <a href="https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#speed-up-and-parallelization-thoughts" rel="nofollow noreferrer">DESeq2 data analysis documentation</a> has a few different suggestions for speeding up analysis:</p>
<blockquote>
<p>The above steps should take less than 30 seconds for most analyses. For experiments with complex designs and many samples (e.g. dozens of coefficients, ~100s of samples), one may want to have faster computation than provided by the default run of DESeq. We have two recommendations:</p>
<ol>
<li>By using the argument fitType=&quot;glmGamPoi&quot;, one can leverage the faster NB GLM engine written by Constantin Ahlmann-Eltze. Note that glmGamPoi’s interface in DESeq2 requires use of test=&quot;LRT&quot; and specification of a reduced design.</li>
<li>One can take advantage of parallelized computation. Parallelizing DESeq, results, and lfcShrink can be easily accomplished by loading the BiocParallel package, and then setting the following arguments: parallel=TRUE and BPPARAM=MulticoreParam(4), for example, splitting the job over 4 cores. However, some words of advice on parallelization: first, it is recommend to filter genes where all samples have low counts, to avoid sending data unnecessarily to child processes, when those genes have low power and will be independently filtered anyway; secondly, there is often diminishing returns for adding more cores due to overhead of sending data to child processes, therefore I recommend first starting with small number of additional cores. Note that obtaining results for coefficients or contrasts listed in resultsNames(dds) is fast and will not need parallelization. As an alternative to BPPARAM, one can register cores at the beginning of an analysis, and then just specify parallel=TRUE to the functions when called.</li>
</ol>
</blockquote>
<p>Also, applying a low count filter can help a lot.</p>
<blockquote>
<p><strong>Pre-filtering</strong></p>
<p>While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, <em><strong>we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2.</strong></em> It can also improve visualizations, as features with no information for differential expression are not plotted.</p>
<p>Here we perform a minimal pre-filtering to keep only rows that have at least 10 reads total. Note that more strict filtering to increase power is automatically applied via independent filtering on the mean of normalized counts within the results function.</p>
</blockquote>
<pre><code>keep &lt;- rowSums(counts(dds)) &gt;= 10
dds &lt;- dds[keep,]
</code></pre>
<blockquote>
<p>Alternatively, a popular filter is to ensure at least X samples with a count of 10 or more, where X can be chosen as the sample size of the smallest group of samples:</p>
</blockquote>
<pre><code># unevaluated chunk...
keep &lt;- rowSums(counts(dds) &gt;= 10) &gt;= X
dds &lt;- dds[keep,]
</code></pre>

CC BY-SA 4.0
2022-12-29T07:44:53.733Z
Awesome, thanks for your input! Pre-filtering peaks that have a minimum count significantly reduces the space complexity. Of course, the million dollar question is: what minimum number of counts should each peak have... 5? 10? 100? Anyway, I think that's a question I should try and answer on my own. Thank you very much!
Shashank Nagaraja
CC BY-SA 4.0
2022-12-29T07:58:19.370Z
Adding to the previous comment, yes I am using parallel compute. I was not, however, using the NB GLM model; I am trying it right now and will report the speedup. w.r.t. using DESeq2 for ATAC-seq samples... I totally agree that the implicit assumptions make it a far-from-perfect algorithm to model the distribution in question. The only reason to use it is that other publications use it, which is a **terrible** reason. I think the Wilcoxon rank-sum test might be more appropriate here.
Shashank Nagaraja
CC BY-SA 4.0
2022-12-20T19:30:19.047Z
2022-12-20T19:30:19.047Z
gringer
<p>I am trying to run a very large number of transposase-accessible chromatin (ATAC)-seq samples through DESeq2 to find peaks of differential chromatin accessible across my study genome.</p>
<p>I have about 100 samples, and 1,000,000 &quot;genes&quot; which are really peaks.</p>
<p>Doing this takes a significant amount of RAM; right now I am running on a 1 TB bigmem node and still getting an out-of-memory error. The code works perfectly fine for a small subset of samples.</p>
<p>Does anyone have an alternative strategy to use DESeq2 on a very large dataset? Should I look into doing a Wilcox-rank test myself? Thanks.</p>

<r><deseq2><differential-expression><edger><atac-seq>
DESeq2 for large number of samples takes too much RAM
# DESeq2 for large number of samples takes too much RAM

(Asked by: username_0 on 2022-12-19)

I am trying to run a very large number of transposase-accessible chromatin (ATAC)-seq samples through DESeq2 to find peaks of differential chromatin accessible across my study genome.

I have about 100 samples, and 1,000,000 "genes" which are really peaks.

Doing this takes a significant amount of RAM; right now I am running on a 1 TB bigmem node and still getting an out-of-memory error. The code works perfectly fine for a small subset of samples.

Does anyone have an alternative strategy to use DESeq2 on a very large dataset? Should I look into doing a Wilcox-rank test myself? Thanks.

**Tags:** r, deseq2, differential-expression, edger, atac-seq

---

## 2 answers

### Answer by username_4 on 2022-12-20. Score: 4 | Accepted answer

I'll preface this by saying that I don't think DESeq2 is the right tool to use for ATAC-Seq data. My own study of ATAC-Seq patterns [admittedly only a couple of runs that were our first exploration of the technique] indicates that it has a fractal expression which is difficult to properly summarise at any scale, especially at a bulk sample level. I found this poster from a quick search on Twitter, which recommends AIAP for ATAC-Seq analysis (but don't read too much into that; I have no experience with doing ATAC-Seq analysis). What I do have a lot of experience with is DESeq2.
DESeq2 has certain assumptions about the expression patterns of genes which are built into its differential expression model, and the further you get from gene-level analysis, the less appropriate DESeq2 is for data analysis (for example, Swish is a better approach for isoform level and allelic analysis).
However, if you're absolutely set on using DESeq2...  the DESeq2 data analysis documentation has a few different suggestions for speeding up analysis:

The above steps should take less than 30 seconds for most analyses. For experiments with complex designs and many samples (e.g. dozens of coefficients, ~100s of samples), one may want to have faster computation than provided by the default run of DESeq. We have two recommendations:

By using the argument fitType="glmGamPoi", one can leverage the faster NB GLM engine written by Constantin Ahlmann-Eltze. Note that glmGamPoi’s interface in DESeq2 requires use of test="LRT" and specification of a reduced design.
One can take advantage of parallelized computation. Parallelizing DESeq, results, and lfcShrink can be easily accomplished by loading the BiocParallel package, and then setting the following arguments: parallel=TRUE and BPPARAM=MulticoreParam(4), for example, splitting the job over 4 cores. However, some words of advice on parallelization: first, it is recommend to filter genes where all samples have low counts, to avoid sending data unnecessarily to child processes, when those genes have low power and will be independently filtered anyway; secondly, there is often diminishing returns for adding more cores due to overhead of sending data to child processes, therefore I recommend first starting with small number of additional cores. Note that obtaining results for coefficients or contrasts listed in resultsNames(dds) is fast and will not need parallelization. As an alternative to BPPARAM, one can register cores at the beginning of an analysis, and then just specify parallel=TRUE to the functions when called.


Also, applying a low count filter can help a lot.

Pre-filtering
While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2. It can also improve visualizations, as features with no information for differential expression are not plotted.
Here we perform a minimal pre-filtering to keep only rows that have at least 10 reads total. Note that more strict filtering to increase power is automatically applied via independent filtering on the mean of normalized counts within the results function.

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]


Alternatively, a popular filter is to ensure at least X samples with a count of 10 or more, where X can be chosen as the sample size of the smallest group of samples:

# unevaluated chunk...
keep <- rowSums(counts(dds) >= 10) >= X
dds <- dds[keep,]

> **Comments:**
> **username_0**: Awesome, thanks for your input! Pre-filtering peaks that have a minimum count significantly reduces the space complexity. Of course, the million dollar question is: what minimum number of counts should each peak have... 5? 10? 100? Anyway, I think that's a question I should try and answer on my own. Thank you very much!  
> **username_0**: Adding to the previous comment, yes I am using parallel compute. I was not, however, using the NB GLM model; I am trying it right now and will report the speedup. w.r.t. using DESeq2 for ATAC-seq samples... I totally agree that the implicit assumptions make it a far-from-perfect algorithm to model the distribution in question. The only reason to use it is that other publications use it, which is a **terrible** reason. I think the Wilcoxon rank-sum test might be more appropriate here.  

---

### Answer by username_1 on 2022-12-20. Score: 1

Absolutely, look at Wilcox-rank, not purely because of the computational overhead, but the parametric assumption that DESeq2 uses [are not correct](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-022-02648-4) on large data sets. This might describe your data sets too:

> **Abstract**
> DESeq2 and edgeR, have unexpectedly high false discovery rates.
>
> Expanding the analysis to limma-voom, NOISeq, dearseq, and **Wilcoxon rank-sum test**, we found that [false discovery rates] FDR control is often failed **except for the Wilcoxon rank-sum test**
>
> Particularly, the actual **FDRs of DESeq2** and edgeR sometimes exceed **20%** when the target FDR is 5%. Based on these results, for population-level RNA-seq studies with large sample sizes, **we recommend the Wilcoxon rank-sum test**.
>
> **Results and Discussion**
> DESeq2 and edgeR had only an 8% overlap in the DEGs they identified

---

> **Comments:**
> **username_2**: I feel it's unfair to write that without mentioning the DESeq author [response](https://twitter.com/mikelove/status/1513468614585765889): in short, there are technical confounders that need to be controlled for. Wilcoxon gets fewer False Positives because of its lower power (while DESeq and edgeR are designed to maximize power), but the right thing to do is really to correct for confounding, whatever test you do next. That being said for this question Wilcoxon is probably reasonable (though a correction with RUV/SVA/... would be appropriate).  
> **username_3**: It is not a matter of parametric assumptions. The linked paper, as said already did not address the fact that in larger cohorts you often find samples that form clusters by some unwanted confounder nested with the variable you're testing. That is not unusual especially in human samples where you have all sorts of confounders, be it eating habits, drug consumption, dietary status etc that one often does not know of. Once this was addressed (RUVseq by the DESeq2 author) the main claim of that paper collapsed. In any case, at larger n I would use linear approaches such as limma-voom.  
> **username_1**: Hi @ATpoint thanks its a good discussion, albeit maybe not this site.  
> **username_3**: @M__, I think these comments are necessary to put answers into context. The paper you reference is basically an advertisement of the non-parametric tool the authors created. It's in the nature of these sort of papers that advertised tool is best and all others are not as good. That goes for any benchmark paper written by method developers, be it aligners, variant callers or stats frameworks. I am not judging this, it is just an observation. Like any of these papers, results should be read with a grain of salt and I think it is important to notice that this was not an unbiased benchmark.  
> **username_1**: @ATpoint we can have this discussion via the site  that you are famous on. There been an objection raised on this site and we must be mindful of that. I will raise this as a question on *your site* not quite now because it's Christmas Day (quite literally), but shortly. I agree 100% is an important topic and I'll explain my interest in the subject (I don't really do eukaryotes BTW).  

---

