Dada2: Assigning taxonomy with new database

[carolina671]

Hi,

I would like to run the dada2 pipeline with a different database, as ai am working with a gene different to 16S, specifically i work with the nifH gene which has its own database. This database can be found in .arb format (https://wwwzehr.pmc.ucsc.edu/nifH_Database_Public/) or in a fasta format ( http://www.css.cornell.edu/faculty/buckley/nifH_database_2012.fasta).

Since the silva file for assigning the 16S taxonomy is in a fasta file .arb format is an updated version, i tried unsuccesfully to convert the .arb file to a fasta file using qiime. Maybe do you have any other ideas on how to make it readable for the pipeline?

Also I tried to assign the taxonomy directly with the fasta file but i am getting the following errors:

NIFHdatabase<- "~/DATA/NIFH/nifH_database_2012.fasta"
taxaNIFH <- assignTaxonomy(seqtab.nochimNIFH, NIFHdatabase, multithread=TRUE)
Error in assignTaxonomy(seqtab.nochimNIFH, NIFHdatabase, multithread = TRUE) :
Incorrect reference file format for assignTaxonomy.
In addition: Warning message:
In .Call2("fasta_index", filexp_list, nrec, skip, seek.first.rec, :
reading FASTA file /frontiers-shared/home/b6001443/DATA/NIFH/nifH_database_2012.fasta: ignored 4679468 invalid one-letter sequence codes

Since the fasta files are aligned and therefore contain some gaps the i tried to remove the gaps on the sequences and the I got the following errors:

NIFHdatabase<- "~/DATA/NIFH/nifH_database_2012.fasta.txt"
Warning messages:
1: replacing previous import ‘BiocGenerics::path’ by ‘Rsamtools::path’ when loading ‘GenomicAlignments’
2: multiple methods tables found for ‘rglist’
3: replacing previous import ‘BiocGenerics::path’ by ‘Rsamtools::path’ when loading ‘ShortRead’
taxaNIFH <- assignTaxonomy(seqtab.nochimNIFH, NIFHdatabase, multithread=TRUE)
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘windows’ for signature ‘"DNAString"’

Thanks!

[benjjneb]

assignTaxonomy expects a specific format of the fasta file that is outlined here: https://benjjneb.github.io/dada2/training.html#formatting-custom-databases

The fasta file you have posted doesn't conform to that format. It needs to be modified by removing the gap characters (just the sequences, not aligned sequences) and by placing the taxonomic rankings in the id lines of the fasta file (e.g. > Rank1;Rank2;Rank3).

[carolina671]

Hi,

Thanks so much for the information.

Do you know if there is any way of converting a .arb file to a fasta file using R?

Thanks,

[benjjneb]

Not that I am aware of, you may want to take a look at Pat Schloss's description of how he creates a Silva database from .arb files, in particular the beginning part about getting sequences out of .arb into a .fasta: http://blog.mothur.org/2018/01/10/SILVA-v132-reference-files/

[carolina671]

Hi Ben,

I have manage to convert the NIFH reference database into the specific format(please find attached).
seqDatabase1.fasta.txt

Nevertheless when i run the code i still get an error:

taxaNIFH <- assignTaxonomy(seqtab.nochimNIFH, "C:/Users/Gloria Betancur/Desktop/seqDatabase1.fasta.txt", multithread=TRUE)
Error in .Call2("new_XStringSet_from_CHARACTER", ans_class, ans_elementType, :
key 79 (char 'O') not in lookup table

Could you please help me?

Thanks!

[benjjneb]

Sorry for delay, was out for the holidays.

I just ran the following code, and got no error (and normal looking output):

fn <- "~/Desktop/seqDatabase1.fasta"
tax <- assignTaxonomy(st, fn, multithread=TRUE)

First, which version of dada2 are you using? packageVersion("dada2")

Second, can you recreate the error on your end? I am using a different set of query sequences, so maybe it is something query specific. What is the output of table(dada2:::C_isACGT(getSequences(seqtab.nochimNIFH)))?

Alternatively, can you also share the query sequences (i.e. your seqtab.nochimNIFH object, or just the sequences form it) so that I can reproduce the error?

[carolina671]

Hi Ben,

I tried to run the command again and didn´t get any error this time. Not sure why...

Anyway I was able to get the taxa table and the otu table!

Thanks a lot for your patience and help!

[moyn413]

Hi, this was a helpful thread to find, as I have the same question regarding the nifH arb database and dada2. @carolina671 is the file you shared (seqDatabase1.fasta.txt) the Zehr database (nifH_2017June.arb)? Thanks!

[EricRaes]

Hi all, many thanks for this thread and @carolina671 the seqDatabase1.fasta.txt file.

Just a quick comment @carolina671 the seqDatabase1.fasta.txt file doesn't have the enigmatic Candidatus Atelocyanobacterium thalassa (previously cyanobacterium UCYN-A) which seems odd?

Did you use the nifH_2017June.arb file and would you mind sharing how you converted the ARB file to a text file?

Many thanks!!

[carolina671]

Hi, this data base was obtained from the taxadiva script (https://github.com/lavanyarishishwar/taxadiva), basiclly i merged the seqDatabase.fasta and the taxonomy.tsv files. I have just checked and I am afraid the Candidatus Atelocyanobacterium thalassa is not in these files.

Sorry for the delayed answer and please let me know if you have any questions. Good luck!

[EricRaes]

Many thanks @carolina671!

[moyn413]

Thanks @carolina671 for clarifying!

@EricRaes, are you going to try to reformat the nifH_2017June.arb file to use with data2? This is something I would like to do, but haven't gotten around to it yet. I would be happy to work together on it.

[EricRaes]

Hi @moyn413
yes that's the plan :-). Currently I used @carolina671 seqDatabase1.fasta.txt file and just added a Candidatus Atelocyanobacterium thalassa sequence.

What's your email, might be easier to communicate. Mine is eric.[email protected]

[moyn413]

Hi, just to follow up on this. I managed to modify, update, and reformat the nifH ARB database for use with dada2. If you find any issues with this, please let me know as it is still new, and I'm sure there is room for improvement!
https://github.com/moyn413/nifHdada2

[benjjneb]

@moyn413 Is it OK for me to list your nifH database as a "contributed" reference file on the DADA2 taxonomy training data page? https://benjjneb.github.io/dada2/training.html

[moyn413]

Sure, @benjjneb that is fine with me.

[sankadinesh]

Hi, just to follow up on this. I managed to modify, update, and reformat the nifH ARB database for use with dada2. If you find any issues with this, please let me know as it is still new, and I'm sure there is room for improvement!
https://github.com/moyn413/nifHdada2

Hi there
I would like to know the parameters to be used while executing assignTaxonomy. Is it necessary to change minboot from 50 to any particular value for nifH. Minboot 50 is set for 16s gene. Will it be suitable for nifH as well ? For example, if one need to do classification using classify.seqs in mothur, one need to provide cut off to classify sequences
Regards,
Dinesh S L

[benjjneb]

@sankadinesh The original publication suggested (If i recall correclty) a minBoot=50 for sequences < 250 nts, and minBoot=80 for sequences > 250nts.

Subsequent publications have suggested that either is an acceptable choice, even for ~200nt sequences, with the understanding that the lower stringency of minBoot=50 results in both an increase in false positives and a decrease in false negatives. See eg: https://peerj.com/articles/4652/

It is true that almost all this work is done on 16S or ITS sequencing data. I do not know of any nifH specific evaluations, but if anyone else does please feel free to share them.

[sankadinesh]

Dear Sir,
I have used 89 as cutoff for MOTHUR (classify.seqs) and minBoot = 89 in DADA2 for same input sequences (ASVs).

The classification percentages from Kingdom to Genus levels were,

|Mothur|DADA2|
|---|---|
|41.64|43.72|
|10.94|10.14|
|6.30|4.07|
|4.03|2.46|
|3.46|1.18|
|3.03|0.33|

The types of genus and family detected was same (e.g. 4 genera) in both methods), whereas it differed only in abundance.

Note:
Both methods are based on Wang 2007, AEM paper, use kmer=8, enabled reversing sequence and cutoff=89. I am wondering regarding the variation in abundance.

Regards,
Dinesh S L

[benjjneb]

@sankadinesh

I can't tell you exactly why the numbers differ. Remember there is a random component in the nBC (the choosing of subsamples during bootstrapping) so some stochastic difference is expected. There also may be technical differences in the implementations that have a non-zero impact. Overall those look like pretty similar results though, which is also what I've seen when compaing assignTaxonomy to other nBC implementations including mothur.

[moyn413]

Gaby et al. 2018 (https://aem.asm.org/content/84/4/e01512-17), have some thresholds they use when assigning nifH taxonomy, although their method of assigning taxonomy is different.

As a next step for the database, I plan to add nifH cluster information for each reference sequence into the taxonomy file. I envision that the cluster bootstrap value from assignTaxonomy could be a useful cutoff value. I don’t think I will get around to this for another couple of months though.