DESeq2 for large number of samples takes too much RAM

Question

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.

Answer 1

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:

1. 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.
2. 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,]

Answer 2

Absolutely, look at Wilcox-rank, not purely because of the computational overhead, but the parametric assumption that DESeq2 uses are not correct 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

---