Dada2: Formatting the entire Silva database

[valdeanda]

Hi!
Excellent pipeline!
Im dealing with some troubles formating the entire silva database. This is what I did.

From what I could guess from the _assignTaxonomy_ and _assignSpecies functions_
I've downloaded these two files from the current Silva database :

1. SILVA_132_SSURef_Nr99_tax_silva_full_align_trunc.fasta.gz : I renamed as silva.nr_v132.align

2. SILVA_132_SSURef_tax_silva.fasta.gz : For me this file looks like the required fasta format file to _assign taxonomy_ just by doing by doing this: cut -f 2 -d " " SILVA_132_SSURef_tax_silva.fasta

   I use the silva.nr_v132.align file to use the makeTaxonomyFasta, to format the silva database in DADA2 format and It gives me the following error.

path <- "../Databases/Silva.nr_v132"
dada2:::makeTaxonomyFasta_Silva(file.path(path, "silva.nr_v132.align"), file.path(path, "silva.nr_v132.tax"), "../Databases/Silva.nr_v132/silva_nr_v132_formated.fa.gz")

Error: cannot allocate vector of size 256.0 Mb
In addition: Warning message:
  In .Call2("fasta_index", filexp_list, nrec, skip, seek.first.rec,  :
              reading FASTA file ../Databases/Silva.nr_v132/silva.nr_v132.align: ignored 211741253 invalid one-letter sequence codes
path <- "../Databases/Silva.nr_v132"
dada2:::makeTaxonomyFasta_Silva(file.path(path, "silva.nr_v132.align"), file.path(path, "silva.nr_v132.tax"), "../Databases/Silva.nr_v132/silva_nr_v132_formated.fa.gz")

Error: cannot allocate vector of size 256.0 Mb
In addition: Warning message:
  In .Call2("fasta_index", filexp_list, nrec, skip, seek.first.rec,  :
              reading FASTA file ../Databases/Silva.nr_v132/silva.nr_v132.align: ignored 211741253 invalid one-letter sequence codes

So Im not sure if im using the correct files to format the database, and what database do you currently recomend (rdp, greengenes or silva) ?
Thanks in advance

Best
Val

[valdeanda]

To format the SILVA_132_SSURef_tax_silva.fasta I tried this (it was not a simple cut)

sed 's/>/> /g' SILVA_132_SSURef_tax_silva.fasta | cut -f 1,3 -d " " | sed 's/ //g'  | awk '/^>/ {$0=$0 ";"}1'   > silva.132.formated.fasta

To get this

>Bacteria;Epsilonbacteraeota;Campylobacteria;Campylobacterales;Helicobacteraceae;Helicobacter;unidentified;
GCAAGUCGAACGAUGAAGCCUAGCUUGCUAGGUUGAUUAGUGGCGCACGGGUGAGUAAUGCAUAGAUGACAUGCCCUUUA

Then I use the assignTaxonomy function

> taxatest  <- assignTaxonomy(seqtab.nochim, "../Databases/Silva.nr_v132/silva.132.formated.fasta", multithread=TRUE)

And it gives me this error..

Warning message:
In .Call2("fasta_index", filexp_list, nrec, skip, seek.first.rec, :
reading FASTA file ../Databases/Silva.nr_v132/silva.132.formated.fasta: ignored 500379136 invalid one-letter sequence codes

[benjjneb]

Silva uses U instead of T, and DADA2 is expecting the DNA code (T).

If you translate the Us to Ts that should fix that problem.

[arthurdegrandpre]

Hi, what would be the easiest way of doing this U to T translation?
Thank you.

[benjjneb]

Hi, what would be the easiest way of doing this U to T translation?

I'm sure there is some sed or awk 1-liner that could do it, but I'm not much with shell scripting.
Alternatively you can read the file into R and translate there. This is one way using ShortRead/Biostrings:

library(ShortRead)
foo <- readFasta("path/to/silva_file_rna.fasta")
seqs.dna <- as(sread(foo), "DNAStringSet") # Translates from RNA code w/ Us to DNA code w/ Ts
writeFasta(ShortRead(seqs.dna, id(foo)), "path/to/silva_file_dna.fasta", width=20000L)