Dada2: Silva 132 vs 138 for identifying cyanobacterial ASVs

[pdcountway]

Hi - I've used assignTaxonomy with the Silva v132 reference database to identify cyanobacterial ASVs. A high proportion of ASVs were identified as Cyanobium PCC-6307. I aligned the sequences of the 10 most abundant 'Cyanobium' ASVs in Geneious and discovered a large amount of sequence variability. I then aligned these 10 ASV sequences to Silva v138 using the online Silva/Sina aligner and the results suggested that only 5 of these ASVs were Cyanobium, while 4 were identified as 'chloroplast' and one belonged to another genus of cyanobacteria, Aphanizomenon. This new taxonomic assignment made a lot more sense in terms of the observed sequence diversity.
I've seen a few posts about updating the dada2 reference database to v138 - is the bottleneck still waiting for the folks at mothur to prepare a v138 database? Thanks!

[benjjneb]

I've seen a few posts about updating the dada2 reference database to v138 - is the bottleneck still waiting for the folks at mothur to prepare a v138 database? Thanks!

Yes, that's the starting point for our current automated Silva database curation script: https://github.com/benjjneb/dada2/blob/master/R/taxonomy.R#L488

We probably should work with the SIlva release files directly, but I'm short on bandwidth right now to make that happen. If anyone curates a Silva v138 database for assignTaxonomy, feel free to share!

[cjfields]

@benjjneb would this be using the mothur versions? It looks like these were just updated:

https://blog.mothur.org/2020/03/04/SILVA-v138-reference-files/

I may have some time to dedicate to this, would a Google Drive link (or similar) work?

EDIT: just reread @pdcountway's message and his note about mothur at the end. Will see what I can generate.

[cjfields]

@benjjneb @pdcountway I generated the files for Silva v138 and added them to Google Drive. Completely untested, let me know if these work.

[pdcountway]

@benjjneb @pdcountway I generated the files for Silva v138 and added them to Google Drive. Completely untested, let me know if these work.

@cjfields Thanks for creating the the Silva v138 files, that's great! I'll re-run assignTaxonomy on my cyano sequences and see if there's an improvement in taxonomic identification. Cheers!

[pdcountway]

@benjjneb I'm just getting back to this and have a basic R question. I have a file saved called seqtab.nochim.csv - generated after the original chimera removal step. I'd like to use this file for taxonomic re-assignment with the new Silva v138 that @cjfields provided. What is the correct way to load this file back into R so that it can be understood by dada2? Or is it not that simple? I tried assigning the results of read.csv to the variable: seqtab.nochim <- read.csv("path to filename/seqtab.nochim.csv") and then running assignTaxonomy per usual but received this error message: Error in getUniques(object, collapse = collapse, silence = silence) : Unrecognized format: Requires named integer vector, fastq filename, dada-class, derep-class, sequence matrix, or a data.frame with $sequence and $abundance columns.
Thanks, Pete

[benjjneb]

Thanks @cjfields ! To clarify, are those the files obtained by running the dada2 scripts makeTaxonomyFastaSilva etc. on the new mothur files? I.e. no additional tweaks or alternative processing?

[benjjneb]

@pdcountway I highly recommend using the saveRDS/readRDS functions for saving and loading data from/into R. They will give you back exactly the same object you saved, e.g. `

saveRDS(foo, "foo.rds")
bar <- readRDS("foo.rds")
identical(foo, bar) # TRUE

With the reading/writing to text formats like csv, a lot of details of the R objects can get changed, such as object classes and storage formats.

Is it possible to regenerate seqtab.nochim and save it that way? If so, that would be easiest. If not, that's OK and we can try to troubleshoot the read.csv step.

[pdcountway]

@benjjneb Thanks for the suggestion on save and read RDS. I will do that in the future. It sounds like the easiest path is to regenerate the seqtab.nochim with the pipeline - just thought I might be able to take a shortcut.

[cjfields]

@benjjneb yes, specifically:

library(dada2)
dada2:::makeTaxonomyFasta_Silva('./mothur/silva.nr_v138.align', './mothur/silva.nr_v138.tax', "./silva_nr_v138_train_set.fa.gz")
dada2:::makeSpeciesFasta_Silva("./silva/SILVA_138_SSURef_tax_silva.fasta.gz", 'silva_species_assignment_v138.fa.gz')

Note the makeSpeciesFasta_Silva followed these directions:

https://github.com/benjjneb/dada2/blob/759e4886fc63b7e533c0997d3cfb697d6fd4d745/R/taxonomy.R#L514

The file SILVA_138_SSURef_tax_silva.fasta.gz was retrieved from the main Silva database:

https://www.arb-silva.de/fileadmin/silva_databases/release_138/Exports/SILVA_138_SSURef_tax_silva.fasta.gz

R information:

> sessionInfo()
R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Catalina 10.15.3

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.6/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] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dada2_1.14.1 Rcpp_1.0.3  

loaded via a namespace (and not attached):
 [1] plyr_1.8.6                  RColorBrewer_1.1-2          pillar_1.4.3                compiler_3.6.2              GenomeInfoDb_1.22.0         XVector_0.26.0              bitops_1.0-6                tools_3.6.2                
 [9] zlibbioc_1.32.0             lattice_0.20-40             lifecycle_0.1.0             tibble_2.1.3                gtable_0.3.0                png_0.1-7                   pkgconfig_2.0.3             rlang_0.4.5                
[17] Matrix_1.2-18               DelayedArray_0.12.2         rstudioapi_0.11             parallel_3.6.2              GenomeInfoDbData_1.2.2      stringr_1.4.0               hwriter_1.3.2               dplyr_0.8.4                
[25] Biostrings_2.54.0           S4Vectors_0.24.3            IRanges_2.20.2              stats4_3.6.2                grid_3.6.2                  tidyselect_1.0.0            Biobase_2.46.0              glue_1.3.1                 
[33] R6_2.4.1                    jpeg_0.1-8.1                BiocParallel_1.20.1         latticeExtra_0.6-29         reshape2_1.4.3              ggplot2_3.3.0               purrr_0.3.3                 magrittr_1.5               
[41] Rsamtools_2.2.3             matrixStats_0.55.0          scales_1.1.0                BiocGenerics_0.32.0         GenomicRanges_1.38.0        GenomicAlignments_1.22.1    ShortRead_1.44.3            assertthat_0.2.1           
[49] SummarizedExperiment_1.16.1 colorspace_1.4-1            stringi_1.4.6               RCurl_1.98-1.1              RcppParallel_4.4.4          munsell_0.5.0               crayon_1.3.4 

[pdcountway]

Hi @benjjneb I ran through the dada2 pipeline again and saved seqtab.nochim with the saveRDS function per your suggestion. The reason I'm doing this is that I keep having problems with dada2 quitting during assignTaxonomy. It may be related to the way we multithread on our supercomputer and we're trying to troubleshoot that - so I'll park that issue for now. However, I'm not sure the *.rds file will be useful to restart assignTaxonomy. Here's why...
I ran through the pipeline in interactive mode, generated a fresh copy of seqtab.nochim (saved a copy for possible use - see below) then ran assignTaxonomy with the new silva v138 database provided by @cjfields. The new reference database loaded and the process began running - using all 32 cores @100%. Several hours later the process crashed - potentially due to how we're queuing up jobs (noted above). When I tried to re-run assignTaxonomy using the saved *.rds file, the new silva v138 failed to load and the process eventually crashed.

This is how I saved seqtab.nochim and eventually tried reloading and rerunning assignTaxonomy.
Initial save:

saveRDS(seqtab.nochim, "seqtab.nochim.rds")
Assigned saved *.rds file to the variable seqtab.nochim:
seqtab.nochim <- readRDS("seqtab.nochim.rds")
Re-launch assignTaxonomy:
taxa <- assignTaxonomy(seqtab.nochim, refFasta = "/path_to_ref_database/silva_nr_v138_train_set_dada2.fasta", multithread = 32)

Normally there's an indication that the database has been read, and then you wait until the assignTaxonomy completes (48-72 hours is typical for me - which I still don't understand). This time there was no indication that the database was read.
Do you have any advice on what may have gone wrong in terms of saving the file and reloading it?
Thanks, Pete

[cjfields]

@pdcountway interesting, I haven't run this yet (still completely untested) so I could try a test run locally to see if I can replicate it. I have a V3-V4 that I can try. I know there were some fairly substantial changed between release v132 and v138.

[benjjneb]

@pdcountway We had a very important bug fix recently in release version 1.14.1 and devel version 1.15.2 that resolved a segfault bug in assignTaxonomy(..., multithread=TRUE): https://github.com/benjjneb/dada2/issues/916

If that is the step that is failing on you, and you are using 1.14.0 (or 1.15.0), you should upgrade to 1.14.1 to resolve that bug.

[pdcountway]

@benjjneb Thanks for the heads up. Yes, we are running 1.14.1 though I never saw the segfault error in 1.14.0, but I suspected it could be related. The crashes I've had are with 1.14.1 and typically the process just ends. I've seen the error message: qsub: read error: Connection timed out, which makes us suspect that it has something to do with our job queuing software PBS Pro.

[pdcountway]

@pdcountway interesting, I haven't run this yet (still completely untested) so I could try a test run locally to see if I can replicate it. I have a V3-V4 that I can try. I know there were some fairly substantial changed between release v132 and v138.

@cjfields Thanks, it would be good to see if you can load a *.rds derived file back in to the assignTaxonomy step. I think your new v138 database is probably OK as I did get it to load in the first step of assignTaxonomy before something else (maybe our job queuing software?) caused it to crash.

[pdcountway]

@cjfields @benjjneb I ran through the process one more time last night and could not get the new silva v138 to load at the assignTaxonomy step. Ultimately the process crashed, ending my R session, no error messages. I will go back and run through the dada2 tutorial with the example data to try to isolate the issues I'm having.

[ac-simpson]

@cjfields @pdcountway I (random user) was able to use the silva v138 training set successfully (running R 3.6.3 on jupyter notebooks on my own laptop). A quick look the new tax table for my top ASVs: the taxonomy is mostly the same compared to silva v132, with some changes at the genus level - so it looks like I'm in business. Thanks a lot!!

[cjfields]

@ac-simpson good to hear! I've been a little pre-occuppied with current events and hadn't tried it yet (it's on my todo list).

[pdcountway]

@ac-simpson @cjfields @benjjneb I have used dada2 v1.14.1 with silva v138 recently to complete data processing through assignTaxonomy - using a batch script on our supercomputer provided by @kguay. The entire process took ~15 minutes for ~500,000 raw sequences and ~600 ASVs. It appears that something about the way I was running dada2 in interactive mode, via our PBS Pro job scheduler, was causing intermittent crashes and not actually utilizing all of the assigned processors. I will provide more details on this later.

[benjjneb]

Thanks @cjfields

We now have Silva v138 up at Zenodo and linked from the main taxonomic training fasta page.

[lideyongJNU]

@benjjneb I have processed PacBio sequences with DADA2 v1.2.1 and assigned taxonomy against the train set on Zenodo. The result looks stranged. For example, one of ASVs was assigned to Gammaproteobacteria; Burkholderiales
ASV_106 Bacteria Proteobacteria Gammaproteobacteria Burkholderiales Nitrosomonadaceae Nitrosomonas NA.
But Order Burkholderiales is under Class Betaproteobacteria according to NCBI. When I search for the train set download from Zenodo, I also found that strange taxonomy assignment.
(See line 549 in the train set)
>Bacteria;Proteobacteria;Gammaproteobacteria;Burkholderiales;Nitrosomonadaceae;MND1;
I am wondering that is there something wrong with the train set?

[and3k]

@lideyongJNU The training set is correct.

SILVA r138 now takes the GTDB into account, which classifies the Burkholderiales into the Gammaproteobacteria: https://gtdb.ecogenomic.org/tree?r=g__Nitrosomonas

The Class Betaproteobacteria does not exist any more in GTDB or SILVA.