Dada2: Best way to make dada2 parallel?

[zachcp]

Hi Team Dada,

Is there a best-practices for running dada2 in parallel? I know that its one of the theoretical selling points but I am not sure of the best way to make use of it. In the tutorial, for example, the derep and dada steps are held as a list so it might be possible to parallelize there. However, dada accepts the list as a whole. If given a dataset too big to fit in memory might it be possible to chunk the data and recombine only the output? I'd love to hear your opinions on this.

Thanks,
zach cp

[joey711]

I'm a fan of parallel nodes in a cluster/cloud. This explicitly avoids the "peak RAM" issue you bring up, which absolutely would be an issue if you were attempting to denoise, say, thousands of large-complexity libraries on one node at one time -- requiring that node to have more memory than roughly the sum of the peak mem requirements of each sample.

Commonly encountered sample complexity-sizes can be efficiently analyzed with a surprisingly small peak RAM, so despite my warning above, very big memory nodes might still be able to easily parallelize an experiment with typical R tools like "foreach" package. Moreover, many cheap commodity nodes can be tasked with analyzing one or more samples in serial, and then writing the results to a central location. Aggregating dada results from thousands of samples can be extremely cheap, fast on a single modest node.

Hope that helps!

joey

I'll close for now, but @benjjneb can still comment if he has anything to add.

[zachcp]

thanks @joey711

my first stab at DADA2 resulted in huge memory use and an incredibly long processing time. However, I was running all the samples at once which keeps quite a bit of extra stuff in memory. Chopping up the samples, processing themselves independently, and saving the per-sample DADA2 outputs as .RDS files is working smashingly. And, as you mentioned this can now be parallelized via parallels' mcmapply and could be modified to run across multiple machines, etc.

However, this raises a second issue: logging. It would be great to be able to log the progress of each sample. Did the error model converge? How many cycles did it take? What were the input parameters? It seems to me that Dada2 + phyloseq has the potential to provide a "full-stack" microbiome analytic capability but it would then be very useful to log the DADA2 process - either as a log file or as data attached to the DADA2 object somehow.

Once again, great work with DADA2,
zach cp

[joey711]

@zachcp thanks for the feedback and kind words. Glad you figured out your parallelization problem. As for the log file, I agree with you, and I typically do keep a record of all the "chatter" that comes out of these runs. I typically parallelize even the earlier sequence filter steps, and save the intermediate files and logs. I recommend this kind of approach, but I could see the utility in the "dada" object itself logging things like the number of iterations to convergence. Maybe @benjjneb can comment on that, too, but this would be pretty easy and very low-mem to add without changing anything else. FYI, the error models at convergence are stored in the dada object, and you can inspect those. There are even plot functions for exploring this, including the self-consistency of complementary substitutions (which were estimated independently, so are a useful measure).

I agree that this is getting close to a "full stack". The obvious missing piece is taxonomic classification/annotation, which in best-practice should include a blacklist garbage removal of things like PhiX, known reagent contaminants, etc. There are some very good tools out there for annotating genotype/rep-seqs - like usearch and pplacer - so unclear if this aspect of the "stack" will be re-coded in R. If usearch ever becomes open source it will be an obvious choice, as I understand it to be low level C that an R package can easily wrap.

Thanks again,

joey

[zachcp]

@joey711,

As for taxonomic classification/annotation, I hadn't seen the new Usearch utax features. Thats a nice tip. One thought: what about the RDP Classifier? I haven't wrapped Java programs and haven't even used RDP but it is Java and the R/Java interop is fine (rCDK is perfectly functional interface to CDK). I think all it would need is exposing one or two methods for passing sequences in and getting classifications back out.

It looks like there is a helper package for importing the RDP output into phyloseq which is mentioned in an RDP tutorial. I might poke around on this although, as you suggest, if Robert Edgar has written software to do what you want he has probably done it faster and better.

[benjjneb]

Running DADA2 on large experiments is something that can be done right now, but needs its own tutorial. @zachcp has definitely hit on some of the key points, but this is the ideal large experiment pipeline as I currently see it:

Filter and Trim: No change.

Estimate errors on a subset of the data: The error parameters can generally be well-estimated on tens of thousands of reads, get better up to hundreds of thousands of reads, and then don't change much with more data. Thus, running dada(...) in selfConsist mode on a subset of the samples can provide the error rates estimates needed in reasonable time and memory.

Use estimated error rates to run dada(...) on each sample individually: Now that the error rates are known, each sample can be run individually inside a parallelized loop (eg. mclapply) with the err parameter provided and selfConsist=FALSE. Memory usage can be kept constant by "streaming" data, i.e. derepFastq and dada(..) inside the loop and save the output to an Rdata or Rds file. Chimera checking can also be run at this stage.

Build the final table: At this point all the dada(...) objects need to be read back in to memory, but they can be reduced to just their associated uniques vector if its really a huge experiment and memory becomes limiting.

Putting up a walkthrough of this process is definitely on the todo list (and I'm reopening this issue for that reason). I'm not terribly familiar with R's parallel processing capabilities, so any further suggestions on that welcome (mclapply seems the way to go?).

[zachcp]

i like the way this pipeline sounds, @benjjneb, and I think for lots of sample it could go very quickly since most of the current time is recomputing the error model for each sample (at least I think thats whats happening).

Parallel on R is pretty straightforward. Here's my current pseudocode. The only thing I noted was that I couldn't max out the CPUs on our server becaues it caused a crash. There was a big CPU spike at the beginning of the job that would settle after a minute or so. I ended up using only about 1/2 of the CPUS. Might be some tuning that would improve that.

library(dada2)
library(parallel)

#get Filenames and Samplenames
fnF  <- list of Forward Fastqs
fnR  <- list of Reverse fastqs
samplenames <-  sapply(strsplit(fastqs, "\\_"), `[`, 1)

runDadaAndSave <- function(forwardReadpath,  reverseReadpath, samplename) {
    outfile <- paste0(samplename,"_dada.RDS")
    filterF  <-
    filterR  <-
    derepF <-
    derepR <- 
    dadaF <- 
    dadaR <- 
    bimeraF <-
    bimeraR <- 
    mergers <- ...
   saveRDS(mergers, outfile)
} 

#parallel mapply here
mcmapply(runDadaAndSave, fnFs, fnRs, samplenames, mc.cores = 10000000, SIMPLIFY = TRUE)

[joey711]

You have 10-million cores??
And how many samples??
Whether Ben's suggestion or mine, you only really need one core per sample (or even fewer cores if you can wait for a few samples to process in serial on the same core).

To be honest, I use a queuing engine to direct the processing of my hundreds or thousands of samples across a farm of nodes, where each node is a separate machine with plenty of RAM. I also use selfConsist=TRUE on each node, and don't bother with a preliminary estimate of error profiles. I think Ben and I respectfully disagree on this point. Here is my reasoning: each sample in my larger projects are run in separate PCRs at a minimum, but often separate sequencer runs, or even separate sequencers. I'm sure a global preliminary error matrix is a great initialization, but I do want the more precise sample-wise estimate from individualized self-consist, and I'm willing to pay the few cents (at most) extra cost per sample to get it.

That said, I agree it would be very useful and informative to include a tutorial on the mcmapply() approach. To be more platform/parallelization-scheme agnostic, I prefer the foreach package. Ben, if you haven't already, the documentation for foreach is very good, and explains how they abstract the parallel implementation from the loop code. Very elegant solution.

My two cents. Good discussion. Sorry I closed this prematurely.

[zachcp]

@joey711 haha. No, no. I think you'd have to scale to Amazon for that many cores. I'm in the hundreds-to-thousands of samples and tens-of-cores range. I'll check out foreach - especially if you think it gives you flexibility down the road. If you don't mind me asking: do you use foreach to farm out to your nodes or is that used only on each local node?

@benjjneb @joey711 I can see the value of both approaches. I think in the most common use case (as determined by only the few projects I've been involved in) you get back per-sample fastq data from the same run. In this case I would expect the error models to be the same. Do you guys think that different PCR reactions using the same reagents would make a difference? I would naively assume not.

Any thoughts on the RDP classifier? i.e. is it worth looking at?

[zachcp]

Note: the source of the massive CPU spike i mentioned earlier is due to the ShortRead Package. As documented in this thread, and now as part of the ShortRead documentation:

FastqSampler and FastqStreamer use OpenMP threads (when available) during creation of the return value. This may sometimes create problems when a process is already running on multiple threads, e.g., with an error message like

    libgomp: Thread creation failed: Resource temporarily unavailable

A solution is to precede problematic code with the following code snippet, to disable threading

    nthreads <- .Call(ShortRead:::.set_omp_threads, 1L)
    on.exit(.Call(ShortRead:::.set_omp_threads, nthreads))

An mcmapply-like solution will need to set this parameter otherwise any sample reaching FastqStreamer will suddenly hog all of the available CPUs.

[luhugerth]

Hey all,
Thanks for the enlightening discussion. My dataset includes samples from many different MiSeq runs, so it's great to be able to calculate error models for each run individually and then combine them back. I've managed to go as far as writing out the .RDS files for each file. When attempting to read them all back with readRDS is where I get the problem:

Error in FUN(X[[i]], ...) :
Unrecognized format: Requires named integer vector, dada-class, derep-class, or a data.frame with $sequence and $abundance columns.

I've tried writing it in a few different ways with the same result. The help for makeSequenceTable has the example

makeSequenceTable(list(sample1=dada1, sample2=dada2))

but how does this extend to an arbitrarily long list of samples?

Sorry for bringing down the level of the discussion, but I'm sure @zachcp already has the code for that, and this will make this thread into a complete reference for future users...

[zachcp]

hi @luhugerth,

I ran into the same problem. You need a named list. Try something like this:

#If my dada files are saved as RDS files in the working directory
filelist         <- list.files(".")
dadalist         <-  Map(readRDS, filelist)
names(dadalist)  <-  filelist #adding the filenames as sample names. you probably want to modify a bit
makeSequenceTable(dadalist)

[benjjneb]

Hi Luisa,
makeSequenceTable needs a list in which each element can be acted on by the
getUniques() function. The error you are getting indicates that getUniques
doesn't recognize some of the elements in the list you provided.

My guess is that you read back in your lists of mergers from each run, and
then gave makeSequenceTable a list of those lists. You need to flatten the
list out, which will probably happen naturally if you use c(...) to combine
the lists, rather than list(...).

If that doesn't work it might help if you clarified what commands you are
using after reading the data back into your R session.

-Ben

On Wed, Nov 4, 2015 at 6:14 PM, Luisa Hugerth [email protected]
wrote:

Hey all,
Thanks for the enlightening discussion. My dataset includes samples from
many different MiSeq runs, so it's great to be able to calculate error
models for each run individually and then combine them back. I've managed
to go as far as writing out the .RDS files for each file. When attempting
to read them all back with readRDS is where I get the problem:

Error in FUN(X[[i]], ...) :
Unrecognized format: Requires named integer vector, dada-class,
derep-class, or a data.frame with $sequence and $abundance columns.

I've tried writing it in a few different ways with the same result. The
help for makeSequenceTable has the example

makeSequenceTable(list(sample1=dada1, sample2=dada2))

but how does this extend to an arbitrarily long list of samples?

Sorry for bringing down the level of the discussion, but I'm sure @zachcp
https://github.com/zachcp already has the code for that, and this will
make this thread into a complete reference for future users...

—
Reply to this email directly or view it on GitHub
https://github.com/benjjneb/dada2/issues/41#issuecomment-153931142.

[luhugerth]

Hi,
@zachcp 's suggestion worked great!
I was using c(...), but I supose I was missing the Map step. That solved it.
Thanks for the help!

[zachcp]

For future users I'm linking to a few helper scripts. I use them in a Makefile-based workflow. Theres one for running dada2 and another for aggregating and saving the results.

https://github.com/zachcp/phyloseq-tools/blob/master/R/scripts/runDada2.R
https://github.com/zachcp/phyloseq-tools/blob/master/R/scripts/agreggateRDSfiles.R

[benjjneb]

The new 1.1 version of the dada2 package has implemented multithreading within the main dada(...) method (see 7472c4cce4d12d83f6ce2531fc66870e54028329).

dada(..., multithread=TRUE) will ask the OS how many threads are available and attempt to use them all. dada(..., multithread=5) will run with the specified number of threads.

I've seen performance gains of up to 10x when processing datasets on a 16 core processor.

We're working on extending the multithreading support to other key functions in the package, and will be working on a new multithreaded workflow. Suggestions welcome!

[zachcp]

BOOM! very nice. I'll try to give this a spin sometime this week. How is the memory consumption?

[joey711]

Mem is a good question. Ben, I suspect peak memory is now O ~ N_core * Mem_prev_peak + Mem_para_overhead? Unless there's something else in the multithreading I don't know about?

[benjjneb]

I haven't tested the memory yet, but I think it will be much better than joey is thinking, and only marginally higher than before. The threads have read-access to the original memory space so it isn't necessary to copy everything to each thread.

[joey711]

Ah, that's great!
Makes sense that this should be possible. Very cool to hear that's already what you implemented.

[bashma]

@benjjneb when will this option be available or when will 1.1 version be released? I tried the development version:
dadaRs <- dada(derepRs, err = inflateErr(tperr1,3), selfConsist = TRUE, multithread=5)
but I only see one thread being used when I look at running processes.

thank you, Mina

EDIT:
packageVersion("dada2")
[1] ‘1.1.3’
dadaFs <- dada(derepFs, err = inflateErr(tperr1,3), selfConsist = TRUE, multithread=5)
Warnmeldung:
In dada(derepFs, err = inflateErr(tperr1, 3), selfConsist = TRUE, :
multithread is not a valid DADA option.

[benjjneb]

@bashma

It's available right now, but you'll need to pull the current devel-version (1.1.1) from github. The Bioconductor devel-branch (the 1.1.3 version you have) is not up to date.

Apologies for the confusion, I should probably resync the github branch and bioconductor-devel. The most current version will always be found at github though.

[elsherbini]

The multithreading option would be great to have on the bioconductor release. I'm running dada on a cluster with an ooooold gcc version and it fails when I try to use devtools to pull from github.

[spholmes]

Joseph,
When you say it fails from github, could you tell us a little more about
your
system and the error message, hard to debug with no symptoms,
Thanks
Susan

On Fri, Jul 29, 2016 at 1:47 PM, Joseph Elsherbini <[email protected]

wrote:

The multithreading option would be great to have on the bioconductor
release. I'm running dada on a cluster with an ooooold gcc version and it
fails when I try to use devtools to pull from github.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
https://github.com/benjjneb/dada2/issues/41#issuecomment-236289244, or mute
the thread
https://github.com/notifications/unsubscribe-auth/ABJcvbkmsBfa8gm4d6oM8gw7puJyaCLOks5qambLgaJpZM4GGK2U
.

Susan Holmes
Professor, Statistics and BioX
John Henry Samter Fellow in Undergraduate Education
Sequoia Hall,
390 Serra Mall
Stanford, CA 94305
http://www-stat.stanford.edu/~susan/

[benjjneb]

The changes in the devel-branch are far too major to move into Bioconductor release. I will work on moving the current devel version on github over to the Bioconductor-devel branch though, so that more up-to-date binaries are available (see Using bioc-devel.

Will work on that today. If all goes well, binaries will be available next week. If the BioC build system balks at the changes (eg. C++11 and BiocParallel) then it could take longer to sort out.

[benjjneb]

@elsherbini Version 1.1.3 of the dada2 package, which includes multithreading, successfully built on all three Bioconductor platforms. You can now install multithreading-enabled binaries through bioc-devel.

[elsherbini]

Thanks for taking interest @spholmes. The issue seems to be that the gcc version I have access to on the cluster is from 2012 (!) and doesn't support the std=c++11 option. Here is a link to the full log from trying to install dada2 (it fails at compiling RcppParallel): gist

Below is the system information.

OS

>uname -a
Linux eofe4 2.6.32-358.el6.x86_64 #1 SMP Tue Jan 29 11:47:41 EST 2013 x86_64 x86_64 x86_64 GNU/Linux

> lsb_release -a
LSB Version:    :base-4.0-amd64:base-4.0-noarch:core-4.0-amd64:core-4.0-noarch:graphics-4.0-amd64:graphics-4.0-noarch:printing-4.0-amd64:printing-4.0-noarch
Distributor ID: RedHatEnterpriseServer
Description:    Red Hat Enterprise Linux Server release 6.6 (Santiago)
Release:        6.6
Codename:       Santiago

gcc

> gcc -version
gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-11)
Copyright (C) 2010 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

R

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

[elsherbini]

@benjjneb I tried installing dada2 with biocLite and couldn't figure out how to get any version other than 1.0.3

Here is the output from doing that:
https://gist.github.com/elsherbini/c882ef536ca1f3a00364a6335729b7a9

[benjjneb]

@elsherbini See the instructions on installing from bioc-devel here: https://www.bioconductor.org/developers/how-to/useDevel/

Basically you need to call the useDevel() function, and then install.

[elsherbini]

Thanks a ton, that worked for me.

[benjjneb]

Also, the below is my current parallelized workflow:

# Read in names of filtered files
fns <- list.files()
fns <- fns[grepl("_F.fq.gz", fns)]
sams <- sapply(strsplit(fns, "_"), `[`, 1)
names(fns) <- sams

# Learn error rates (usually aim for subset of about 5M total reads)
ii <- sample(length(fns), 30)
drp.learn <- derepFastq(fns[ii])
dd.learn <- dada(drp.learn, err=NULL, selfConsist=TRUE, multithread=TRUE)
err <- dd.learn[[1]]$err_out
rm(drp.learn)
rm(dd.learn)

# Infer sequence variants
dds <- vector("list", length(sams))
names(dds) <- sams
for(sam in sams) {
    derep <- derepFastq(fns[[sam]], verbose=TRUE)
    dds[[sam]] <- dada(derep, err=err, multithread=TRUE)
}

# Make table and remove chimeras
seqtab <- makeSequenceTable(dds)
seqtab <- removeBimeraDenovo(seqtab, multithread=TRUE, verbose=TRUE)
saveRDS(seqtab, "seqtab.rds")

On a fairly typical compute node (16 cores, 64GB of memory), this processed an entire Hiseq run of forward reads (~120M 250bp reads, ~750 samples) in 2 hours.

[elsherbini]

I was able to install version 1.1.3 of dada2, but it says that multithread is not a valid option:

> library(dada2)
Loading required package: Rcpp
> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS release 6.8 (Final)

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dada2_1.1.3 Rcpp_0.12.6

loaded via a namespace (and not attached):
 [1] magrittr_1.5               XVector_0.12.1            
 [3] GenomicRanges_1.24.2       BiocGenerics_0.18.0       
 [5] zlibbioc_1.18.0            GenomicAlignments_1.8.4   
 [7] IRanges_2.6.1              munsell_0.4.3             
 [9] BiocParallel_1.6.3         colorspace_1.2-6          
[11] lattice_0.20-33            stringr_1.0.0             
[13] hwriter_1.3.2              plyr_1.8.4                
[15] GenomeInfoDb_1.8.3         tools_3.3.1               
[17] SummarizedExperiment_1.2.3 parallel_3.3.1            
[19] grid_3.3.1                 data.table_1.9.6          
[21] Biobase_2.32.0             gtable_0.2.0              
[23] latticeExtra_0.6-28        reshape2_1.4.1            
[25] RColorBrewer_1.1-2         ggplot2_2.1.0             
[27] S4Vectors_0.10.2           bitops_1.0-6              
[29] stringi_1.1.1              scales_0.4.0              
[31] Biostrings_2.40.2          Rsamtools_1.24.0          
[33] ShortRead_1.30.0           stats4_3.3.1              
[35] chron_2.3-47              
> dada(multithread=TRUE)                                                                                                                                                
Error in dada(multithread = TRUE) : 
  argument "derep" is missing, with no default
In addition: Warning message:
In dada(multithread = TRUE) : multithread is not a valid DADA option.

[elsherbini]

Unrelated to the above issue: It'd be nice to calculate the error matrix with a random subsampling of reads, rather than a subset of the samples. Is this easily accomplished in R? I was thinking about doing it separately and writing out the subsample as a fastq file using an approach like this from the shell.

[benjjneb]

@elsherbini I did not check closely enough before: The devel-branch has not yet produced an updated build since I checked in the multithreading functionality. The current devel build is still the pre-multithreading version.

Not sure why this is, but will look into this and let you know when it does become available.

As for a random subsampling of reads, there is currently no convenience function in dada2 for doing this, but it can be accomplished in R by subsampling reads from various fastq files using the ShortRead package. Also by some command line utilities (eg. usearch/vsearch).

[joey711]

"It'd be nice to calculate the error matrix with a random subsampling of reads, rather than a subset of the samples."

Procedurally, this should definitely be a different issue post.

Theoretically, I think this warrants some caution. A poorly-defined subsampling of reads across samples might not provide a representative set for estimating the error matrix. There are at least a few criteria that must be satisfied, and @benjjneb will likely enumerate these better than me on the new thread once it is posted O:-).

Quick answer: using just a few samples from the same run and read direction does a pretty good job ensuring these criteria are met. Hence, my preliminary question would be, why is this sample-wise approach not sufficient? What problem are you trying to solve?

[benjjneb]

@elsherbini Builds of the latest version of dada2 are now installable from bioc-devel. Go ahead and reinstall through biocLite (after making sure you are on the devel branch: useDevel). Check the version, it should be 1.1.5 after reinstallation.

Contains multithreading and a couple other new features (eg. assignSpecies: http://benjjneb.github.io/dada2/species.html).

[elsherbini]

Thanks. Unfortunately that new version is not compatible with my gcc, whereas the old version is. So for now I'm back to dada2 version 1.0.3

[benjjneb]

Just pushed the DADA2 "Big Data" workflow today: http://benjjneb.github.io/dada2/bigdata.html

For now I consider that the answer to parallelization best practices, so closing this issue.

[elsherbini]

From the docs:

"Finally, a powerful computing approach that is enabled by the parallelism in this workflow is farming out chunks of the computationally intensive sample inference stage to different nodes, each of which has low resource requirements. 10 Amazon 8GB instances can do the job as well as one larger (and more costly) compute node!"

If I have access to hundreds/thousands of nodes, is it possible to do the inference step for each sample independently and then average the $err_out's? Is there a function in the library for averaging the $err_out's?

[elsherbini]

Going further, since I do have lots of nodes, the best workflow for me would be to keep every sample separate for as long as possible. Besides needing to average together the inferred error rates from many samples, is there any reason you couldn't keep the samples processing independently of one another all the way from quality filtering -> derep -> dada2 denoising (given the error matrix) -> merge -> make sequence table, then just merge the sequence tables at the end?

[joey711]

@elsherbini

Yes, you can do this. I even created a function for merging results from read-pairs that were processed separately. For some reason, @benjjneb hasn't exposed this as an exported function in the package, but it is there, as an internal function. Use the dada2::: prefix to the function to use it.

There's some debate about the utility of pooling. I think the answer depends a lot on how similar the error-rates and real sequences are among the set of samples being pooled. If either of these things is quite different among your set, and each sample is large enough to converge on good estimates of the error rates, then analyzing each sample independently in distributed independent operations is an effective approach. I've done this.

[benjjneb]

@elsherbini

is there any reason you couldn't keep the samples processing independently of one another all the way from quality filtering -> derep -> dada2 denoising (given the error matrix) -> merge -> make sequence table, then just merge the sequence tables at the end?

You can definitely do that.

To take it all the way, do your chimera removal on a per-sample basis as well.