3
$\begingroup$

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.

$\endgroup$
0

2 Answers 2

4
$\begingroup$

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,]
$\endgroup$
2
  • 1
    $\begingroup$ 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! $\endgroup$ Dec 29, 2022 at 7:44
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 29, 2022 at 7:58
1
$\begingroup$

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


$\endgroup$
5
  • 2
    $\begingroup$ I feel it's unfair to write that without mentioning the DESeq author response: 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). $\endgroup$
    – Alexlok
    Dec 20, 2022 at 15:00
  • 3
    $\begingroup$ 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. $\endgroup$
    – ATpoint
    Dec 25, 2022 at 8:42
  • $\begingroup$ Hi @ATpoint thanks its a good discussion, albeit maybe not this site. $\endgroup$
    – M__
    Dec 25, 2022 at 15:12
  • 3
    $\begingroup$ @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. $\endgroup$
    – ATpoint
    Dec 25, 2022 at 17:40
  • $\begingroup$ @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). $\endgroup$
    – M__
    Dec 25, 2022 at 18:21

Your Answer

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

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