Dada2: Improve and document memory scaling of assignTaxonomy

[griffinp]

Hello,

I am wondering whether there's a limit on the size of reference database that can be used, or if perhaps there's a problem with multithreading in assignTaxonomy.

I have made my own custom reference database. It's properly formatted (have tested with subsets of 100,000 reference sequences and it works fine). However, when I try to run assignTaxonomy on 24 samples against the entire reference database (589,276 reference sequences) the job gets killed by SLURM and I get an error:

error: Exceeded step memory limit at some point.

I am running this job on a compute cluster, using 16 cores with 32 GB memory per core. (It would be good to get away 16 GB per core if possible!) I'm requesting multithreaded mode: my assignTaxonomy command is:

taxa_assignment <- assignTaxonomy(seqtab.nochim, "path_to/file.fasta.gz",
multithread=TRUE, taxLevels = c("Kingdom", "Phylum", "Class",
"Order", "Family", "Genus", "Species", "ID"))

Using dada2 v1.4 in R 3.4.0.

Any ideas would be very welcome!

[benjjneb]

Unfortunately memory scaling with large datasets (100s of thousands of sequences+) is an issue with the current implementation of assignTaxonomy.

Guesstimating based on the numbers you've provided, I think it might work with 64GB of memory, but I can't guarantee it. The other workaround you can try is to trim your reference database to the region you sequenced, although that would require an external program as we don't have tools w/in the package to do that.

Apologies for the lack of an easy answer. Scaling of assignTaxonomy is something we should consider revisiting for our next release.

[griffinp]

Thanks for getting back to me so quickly @benjjneb. Indeed it would be fantastic if the scaling of assignTaxonomy can be improved in future releases - in the meantime it would be useful if you can indicate somewhere in the documentation these kinds of limitations (it takes quite a lot of effort to create this kind of custom database, as I'm sure you can imagine!)

I am trying to find a way to get around this problem for my case - I have some "less important" taxonomic groups in the database that can be represented by a subset of sequences (rather than all sequences available in GenBank/BOLD, which I had originally included). Will let you know how it goes.

If you can suggest any alternative approaches for assigning taxonomy in the meantime, that'd be super useful. Ideally without having to make too many changes to the reference database format! I looked at the QIIME assign_taxonomy.py requirements and the formatting looks like a bit of a nightmare...

Thanks again!

[benjjneb]

Unfortunately I'm just not familiar enough with various taxonomic assignment implementations to know how their memory scales. Best I can do is offer two suggestions on places to look:

Robert Edgar usually writes fairly performant software, so the sintax command (which is essentially identical to the RDP classifier for common use-cases) might do the job. However, it probably also requires money, as the usearch software is on a freemium model and the free-version is limited to ~3Gb of memory.

The new microclass R package was written specifically to scale taxonomic classification out to large datasets, both on the sequence side and on the reference side. I can confirm it is very fast. Unfortunately they have not implemented the full RDP classifier (no confidences calculated), so you have to rely on their own k-mer classifer (which is probably fine, but non-standard).

[benjjneb]

Should revisit the C-side implementation to improve memory scaling with reference sequence number, see also #390

[benjjneb]

Comment from another thread:

First thing I can say is that this update [Silva 132] makes the whole thing a lot more computationally strenuous.
addSpecies() crashed the R session after a while on a dataset with ~70k SVs, 120 samples.

[GeoMicroSoares]

Just to further comment on the SILVA_132 issue, addSpecies() on my machine used 33GB RAM + 66GB swap before it crashed the session. assignTaxonomy() takes a while longer with this update, ~2h until completion.

[benjjneb]

addSpecies() on my machine used 33GB RAM + 66GB

Woah, really? I have not seen that kind of memory usage from assignSpecies

Can you give me a bit more info: packageVersion("dada2"), R.version, dim(seqtab) and table(nchar(getSequences(seqtab)))

[GeoMicroSoares]

@benjjneb here's the info you requested:

> ps
phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 69959 taxa and 117 samples ]
sample_data() Sample Data:       [ 117 samples by 80 sample variables ]
tax_table()   Taxonomy Table:    [ 69959 taxa by 6 taxonomic ranks ]

> packageVersion("dada2")
[1] ‘1.6.0’

> R.version
               _                           
platform       x86_64-pc-linux-gnu         
arch           x86_64                      
os             linux-gnu                   
system         x86_64, linux-gnu           
status                                     
major          3                           
minor          4.3                         
year           2017                        
month          11                          
day            30                          
svn rev        73796                       
language       R                           
version.string R version 3.4.3 (2017-11-30)
nickname       Kite-Eating Tree        

> dim(seqtab)
[1]   128 70937

> table(nchar(getSequences(seqtab)))

  260   261   262   263   264   266   267   268   269   270   272   273   274   275   276   277   278 
  357    74     3     4     1     1     6     1     2     1     5     1     2     3     6     9     3 
  279   281   284   285   286   287   288   290   291   292   293   294   295   296   297   298   299 
    2     2     2     2     2     2     2     2     4    11     7     2     6     1     3     2     2 
  302   303   304   306   307   308   310   312   313   314   315   316   317   318   319   320   322 
    2     3     1     3     3     4     1     2     2     2     5     4     4     4     1     4     1 
  324   325   326   327   328   329   330   331   332   333   334   335   336   340   341   342   343 
    1     3     1     2     2     2     5     2     3     3     2     4     1     6     4     4    17 
  344   345   346   347   348   349   351   352   356   357   358   359   360   361   362   363   364 
    4     3     4     1     2     1     2     3    12     7     3     7     6     1     5     4     2 
  365   366   367   368   369   370   371   372   373   374   375   376   377   378   379   380   381 
   14     5     4     3     4     7     7    59    75   166    83     5     5     6     4    23   185 
  382   383   384   385   386   387   388   389   390   391   392   393   394   395   396   397   398 
   66   692   744   823   521   480   108   171    56    12    52    31    51   169   155   187   573 
  399   400   401   402   403   404   405   406   407   408   409   410   411   412   413   414   415 
  158    58   881  7338  6590  4127  3015  1961  1513   469   446   309   504   137   110   501   151 
  416   417   418   419   420   421   422   423   424   425   426   427   428   429   430   431   432 
  442   763   589   755   646   730  4072  1423   707   378  2883 10956  9862  1879   177    44     9 
  433   434   435   436   437   438   439 
   17     7    14    19    16    10    17 

[vhertzb]

Hi Benjjneb, I am trying to run the tutorial, just to make sure everything is in working order, before I let it loose on a different set of fastq's. When I tried to execute the addSpecies step, I got an error as follows:

Error in .Call2("new_XStringSet_from_CHARACTER", ans_class, ans_elementType, : key 70 (char 'F') not in lookup table

and the traceback table is (not numbered right in the preview, but they go from 10 down to 1):

1. | .Call2("new_XStringSet_from_CHARACTER", ans_class, ans_elementType, x, start(solved_SEW), width(solved_SEW), get_seqtype_conversion_lookup("B", seqtype), PACKAGE = "Biostrings")
   -- | --

2. | .charToXStringSet(seqtype, x, start, end, width, use.names)
   -- | --

3. | XStringSet("DNA", x, start = start, end = end, width = width, use.names = use.names)
   -- | --

4. | XStringSet("DNA", x, start = start, end = end, width = width, use.names = use.names)
   -- | --

6.
DNAStringSet(x)

5.
.PDict(x, max.mismatch, tb.start, tb.end, tb.width, algorithm, skip.invalid.patterns)

4.
PDict(seqs[lens == len])

3.
PDict(seqs[lens == len])

1. | assignSpecies(seqs, refFasta = refFasta, allowMultiple = allowMultiple, tryRC = tryRC, verbose = verbose)
   -- | --

2. | addSpecies(seqtab.nochim, "~/Documents/NRSG_741/MiSeqData/MiSeq_SOP/silva_species_assignment_v132.fa.gz")
   -- | --

The session info is as follows:

R version 3.4.3 Patched (2018-02-03 r74240)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] phyloseq_1.22.3 ggplot2_2.2.1
[3] ShortRead_1.36.0 GenomicAlignments_1.14.1
[5] SummarizedExperiment_1.8.1 DelayedArray_0.4.1
[7] matrixStats_0.53.1 Biobase_2.38.0
[9] Rsamtools_1.30.0 GenomicRanges_1.30.1
[11] GenomeInfoDb_1.14.0 Biostrings_2.46.0
[13] XVector_0.18.0 IRanges_2.12.0
[15] S4Vectors_0.16.0 BiocParallel_1.12.0
[17] BiocGenerics_0.24.0 dada2_1.6.0
[19] Rcpp_0.12.15

loaded via a namespace (and not attached):
[1] reshape2_1.4.3 splines_3.4.3 lattice_0.20-35
[4] rhdf5_2.22.0 colorspace_1.3-2 mgcv_1.8-23
[7] yaml_2.1.16 survival_2.41-3 rlang_0.1.6
[10] pillar_1.1.0 RColorBrewer_1.1-2 GenomeInfoDbData_1.0.0
[13] foreach_1.4.4 plyr_1.8.4 stringr_1.3.0
[16] zlibbioc_1.24.0 munsell_0.4.3 gtable_0.2.0
[19] hwriter_1.3.2 codetools_0.2-15 labeling_0.3
[22] latticeExtra_0.6-28 knitr_1.19 permute_0.9-4
[25] biomformat_1.6.0 scales_0.5.0 vegan_2.4-6
[28] RcppParallel_4.3.20.2 jsonlite_1.5 stringi_1.1.6
[31] grid_3.4.3 ade4_1.7-10 tools_3.4.3
[34] bitops_1.0-6 magrittr_1.5 lazyeval_0.2.1
[37] RCurl_1.95-4.10 tibble_1.4.2 cluster_2.0.6
[40] ape_5.0 pkgconfig_2.0.1 MASS_7.3-48
[43] Matrix_1.2-12 data.table_1.10.4-3 iterators_1.0.9
[46] multtest_2.34.0 igraph_1.1.2 nlme_3.1-131.1
[49] compiler_3.4.3

For what it's worth, I have pretty much the same setup as you recommend, I just don't have the silva files nested in a training folder within the same folder as the fastq files.

Any advice on whatever stupid thing I have done is appreciated.

Thanks, Vicki

[vhertzb]

BTW, I just updated everything, ran the addSpecies command again, go the same failure. I then tried assignSpecies and got no failure. I will try to dplyr the two objects taxa, resulting from assignTaxonomy, and genus.species, resulting from assignSpecies until somebody figures out what is wrong with addSpecies.

[benjjneb]

On the tax-table that is failing addSpecies can you give the output of:

dim(taxtab)
head(colnames(taxtab))
head(rownames(taxtab))

[vhertzb]

Yesterday I just used cbind to put the species column onto taxa, and all was ok.

Today I had to update to do something else. I've been in meetings all day, so I just ran this again to get the outputs you requested above, and now it works. HOWEVER, I'm not getting the Species header for that column. I know I can add it, but that is not very elegant.

So I have somehow solved the issue I started with. I still have this other issue (header for species column), but in the grand scheme of things, it is very minor. See below.

```{r}

Assign species

taxa <- addSpecies(taxa, "~/Documents/NRSG_741/MiSeqData/MiSeq_SOP/silva_species_assignment_v132.fa.gz")

species <- as.matrix(genus.species[,2])

taxa <- cbind(taxa1, species)

Inspect the taxonomic assignments:

```{r}

taxa.print <- taxa #Removing sequence rownames for display only
rownames (taxa.print) <- NULL
head(taxa.print)

 Kingdom    Phylum          Class         Order

[1,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"
[2,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"
[3,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"
[4,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"
[5,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"
[6,] "Bacteria" "Bacteroidetes" "Bacteroidia" "Bacteroidales"
Family Genus
[1,] "Muribaculaceae" NA NA
[2,] "Muribaculaceae" NA NA
[3,] "Muribaculaceae" NA NA
[4,] "Muribaculaceae" NA NA
[5,] "Bacteroidaceae" "Bacteroides" NA
[6,] "Muribaculaceae" NA NA

[benjjneb]

Glad it works now.

still have this other issue (header for species column), but in the grand scheme of things, it is very minor.

I have seen the species column being wonky myself as well. I'll add a bug tag to this issue and take a closer look when I get time.

[benjjneb]

Just to further comment on the SILVA_132 issue, addSpecies() on my machine used 33GB RAM + 66GB swap before it crashed the session. assignTaxonomy() takes a while longer with this update, ~2h until completion.

Thanks for this report @GeoMicroSoares it should be much better in 82bfb9c0613d11619ef0c32466618368ecb04ca8

I hadn't really profiled the memory usage of the new assignSpecies, and it really gets out of hand as you get into the many thousands of sequences. Now there is a parameter n=2000 (by default) that controls the number of sequences assigned at a time. Speeds things up in practice, as much less memory swapping for those w/o infinite memory.

Let me know if it works for you though!

[ricmedveterinario]

Hello benjjneb,

In order for me to be able to use the new parameter in assignSpecies of n = 2000 is it necessary to reinstall the dada2 package? my package "dada2" is the most current and does not recognize this parameter as valid.

What I have noticed is that the multithreading processors are not working at this stage, with only RAM being triggered exponentially which generates the locking in a lopping way of the R. system.

It would not be necessary to trigger "multithread = TRUE", as in "assignTaxonomy", for "assignSpecies" and "addSpecies"?

I am trying to run assignSpecies with 20 GB RAM and 8 "cores" of processors and it is not possible to finish the task by the system crash. In "assignSpecies" and "addSpecies" there is 1 "core" and all RAM.

Thank you.

[benjjneb]

In order for me to be able to use the new parameter in assignSpecies of n = 2000 is it necessary to reinstall the dada2 package? my package "dada2" is the most current and does not recognize this parameter as valid.

Yes, you'll need to install the current devel version from github, see instructions on that here (install from source, using devtools): https://benjjneb.github.io/dada2/dada-installation.html

Our next release (1.8) will be at the end of April, and so this feature will be available in the release version then.

What I have noticed is that the multithreading processors are not working at this stage, with only RAM being triggered exponentially which generates the locking in a lopping way of the R. system.

Yes, assignSpecies uses a very fast algorithm, but the memory scaling is the problem. The new parameter "chunks" the processing into pieces that should keep memory requirements reasonable. No multithreading, but it isn't necessary and wouldn't really help anyway as its memory that is the limiting problem right now.

[benjjneb]

The big fix to assignSpecies memory scaling is in for 1.8. There are improvements that could be made to assignTaxonomy as well in the longer term.

[ricmedveterinario]

Hello benjjneb,

Thanks for the help first,

I am able to run my data here normally,

DADA2 is getting every update better,

[FabianRoger]

Hej benjjneb.

I run assignTaxonomy on 16 GB machine with DADA2 version ‘1.10.1’ and a table with 13140 ASVs using the silva_species_assignment_v132.fa as reference DB.

I always get Error: vector memory exhausted (limit reached?) error, even with the n parameter set to as low as n = 100.

any thoughts?

[benjjneb]

silva_species_assignment_v132.fa is the reference database for the assignSpecies command, not the assignTaoxnomy command. You'll want to use the silva_nr_v132_train_set.fa.gz reference db for assignTaxonomy.

[FabianRoger]

Sorry, my bad. I am running the add_species command and experience the behaviour above.

[benjjneb]

Can you try the following and see if you still get the memory error?

assignSpecies(getSequences(my_seqtab)[1:100], "silva_species_assignment_v132.fa", ...)
# ... -> whatever other flags you were using

[alexisrapin]

Thank you benjjneb for your support and this awesome tool.

Here is an insight into the peak memory usage of assignSpecies():

R version: 3.5.1 (2018-07-02)
platform: x86_64-pc-linux-gnu
packageVersion("dada2"): 1.10.0

assignSpecies(getSequences(seqtab)[1:75], "silva_species_assignment_v132.fa.gz")    2.576GiB
assignSpecies(getSequences(seqtab)[1:100], "silva_species_assignment_v132.fa.gz")   2.89GiB
assignSpecies(getSequences(seqtab)[1:1000], "silva_species_assignment_v132.fa.gz")  5.3GiB
assignSpecies(getSequences(seqtab)[1:10000], "silva_species_assignment_v132.fa.gz") 19.62GiB
assignSpecies(getSequences(seqtab)[1:20000], "silva_species_assignment_v132.fa.gz") >30.09GiB (reached memlimit)

This looks like the default n=2000 does not actually impact on peak memory usage, does it?

In the meantime, what about a simple incremental workaround like:

# Define seqtab chunks of size 10000
chunk.size <- 10000
chunks <- split(c(1:ncol(seqtab)),
                sort(c(1:ncol(seqtab))%%ceiling(ncol(seqtab)/chunk.size)))

# Run assignSpecies() on each chunk incrementally
chunks.species <- lapply(chunks,
                         function(x){
                           return(assignSpecies(getSequences(seqtab)[x],
                                                "silva_species_assignment_v132.fa.gz"))
                         })

# Combine the results
do.call(rbind, chunks.species)

Or, in a more condensed fashion:

chunk.size <- 10000
do.call(rbind,
        lapply(split(c(1:ncol(seqtab)),
                     sort(c(1:ncol(seqtab))%%ceiling(ncol(seqtab)/chunk.size))),
               function(x){
                 return(assignSpecies(getSequences(seqtab)[x],
                                      "silva_species_assignment_v132.fa.gz"))
                 }))

Peak memory usage was 28.73GiB with ncol(seqtab) = 46776.

Havn't been looking in depth at the algorithm used by assignSpecies() yet...any obvious issue with this approach?

[benjjneb]

Havn't been looking in depth at the algorithm used by assignSpecies() yet...any obvious issue with this approach?

No that should work fine. What is frustrating (on my end) is that is what the method should be doing already based on the provided n. Hm... I wonder if there is something going on with R's garbage collection I don't understand here. Marking this as a bug.

[benjjneb]

@CarlyMuletzWolz Have you tried asking this question on the GitHub repository that code is coming from already? If not can you try that first?

If no response, please come back here (often folks will be happy to know someone is using their code, which is another reason to ask at the source first).

[benjjneb]

To follow up on this, the memory scaling of assignTaxonomy is most importantly dictated by the number of unique taxonomies in the reference database. So, if there are more unique taxonomies, the memory usage will scale roughly linearly with that increase. In this case a unique taxonomy is a unique set of taxonomic levels... this number tends to increase with larger dataset sizes (but less than linearly), but increases super-linearly with the taxonomic level included (e.g. databases at that go to species-level will have substantially more memory requirements than those that stop at the genus-level).