Dada2: Running assignTaxonomy against multiple references

[NBBense]

I seem to be having an issue where I am unable to get the assignTaxonomy() function working properly with multiple reference files (UNITE full fasta release and Silva v132).

I have aggregated Fungi, Archaeon, and Bacteria reads into a single ASV. If I run the function with the UNITE preceding Silva as follows:

taxa <- assignTaxonomy(seqtab.nochim, c("sh_general_release_dynamic_01.12.2017.fasta","silva_nr_v132_train_set.fa.gz"), multithread=TRUE)

It appears all of my prokaryotic reads come up NA_DADA2_UNSPECIFIED, and I receive the following warning message:

Warning message:
In matrix(unlist(strsplit(genus.unq, ";")), ncol = td, byrow = TRUE) :
data length [108646] is not a sub-multiple or multiple of the number of rows [15521]

If I try running them in the reverse order as follows:
taxa <- assignTaxonomy(seqtab.nochim, c("silva_nr_v132_train_set.fa.gz","sh_general_release_dynamic_01.12.2017.fasta"), multithread=TRUE)

Now it appears my taxa table results contain the fasta header in the kingdom column:
Claussenomyces_sp|KY689630|SH643228.07FU|reps|k__Fungi p__Ascomycota c__Leotiomycetes o__Helotiales f__Helotiaceae

It shouldn't be a huge issue to clean up the files from there, but I am curious if there are any best practices for checking against multiple references within Dada2 or possibly any recommended upstream processing? Our desired outputs are aggregate ASV/taxa tables that include every domain involved in the study.

[benjjneb]

Providing multiple reference files to assignTaxonomy is not supported, and will definitely cause problems when combining UNITE w/ other dbs as there is special sniffing/parsing that is done w/ the UNITE format.

What I'd recommend is to perform taxonomic assignment on each db individually, with outputBootstraps=TRUE. That gives you a return value with $tax equal to the normal taxonomic table, and $boot equal to the bootstrap confidence values at each level. Then construct a merged table taking the entry for each sequence with the highest bootstrap value. Something like...

{r} tax1 <- assignTaxonomy(sq, ref=ref1, outputBootstraps=TRUE, m=TRUE) tax2 <- assignTaxonomy(sq, ref=ref2, outputBootstraps=TRUE, m=TRUE) tax2.better <- tax2$boot[,6] > tax1$boot[,6] # Using bootstrap confidence at level 6 tax <- tax1 tax[tax2.better,] <- tax2[tax2.better,]

Although I'd poke that output a little bit to make sure its working as intended.

[NBBense]

Thank you very much! Just an FYI in case anyone else has this question, I had to do a little bit of data manipulation to get this to work (also ran it with tryRC=TRUE which seemed to replace a few more likely erroneous eukaryotic alignments from the silva reference that did not previously score higher in the UNITE reference):

tax1 <- assignTaxonomy(sq, "silva_nr_v132_train_set.fa.gz", outputBootstraps=TRUE, multithread=TRUE, tryRC=TRUE)
tax2 <- assignTaxonomy(sq, "sh_general_release_dynamic_01.12.2017.fasta", outputBootstraps=TRUE, multithread=TRUE, tryRC=TRUE)
tax2.better <- tax2$boot[,6] > tax1$boot[,6] # Using bootstrap confidence at level 6
tax <- as.matrix(tax1[[1]])
tax<-addSpecies(tax, "silva_species_assignment_v132.fa.gz", tryRC = TRUE)
tax2 <- as.matrix(tax2[[1]])
tax[tax2.better,] <- tax2[tax2.better,]

[benjjneb]

Thanks for cleaning up my code and letting us know it worked!