Train a taxonomy classifier by extracting an amplicon region from a curated sequence alignment (using SILVA as an example).

Community Contributions Tutorials
, , , , ,
1

:construction: Please consider this tutorial a living document , which may change based upon community feedback and ongoing plugin development. :construction:

Overview:

There are several approaches within QIIME 2, and the RESCRIPt plugin, that enable users to extract an amplicon region of interest from larger or full-length reference sequences to make a taxonomy assignment classifier. These include the following:

  1. Primer-pair extraction approach: using primer pairs to find and extract the amplicon region from a larger reference sequence. However, this approach may not return the targeted region from all of the reference sequences, due to primer bias. That is, just as PCR primer pairs are biased towards the unintended preferential amplification of some DNA targets, and the exclusion of other DNA targets, the same is true bioinformatically. That is, not all reference sequences will contain matches to the primer pair being used, resulting in the loss of reference sequences that do contain sequence information for that amplicon region, but failed to be retained due to primer mismatches. Confusingly, PCR primer pairs might actually amplify DNA targets one would not expect, despite too many mismatches betwen the primers and the DNA target. Yet, bioinformatically, we might fail to extract that amplicon region from a potential reference sequence in our database, as the primers are not a great “string match” to the reference seuences, i.e. too many mismatches (even with adjustments to the mismatch parameters of the tool you you're using). This highlights the differences between real-word PCR chemistry and bioinformatic pattern matching. However, some consider these to not be an issue as you are, in effect, creating a reference databases with similar biases as in the real world. For more details on this approach please see: feature-classifier extract-reads.

  2. Extract sequence segments: to extract an amplicon region when the primer-pair extraction approach may not suffice, i.e. the reference sequences do not contain the primer sequence. This can occur for several reasons: a) The primers were trimmed / removed from the sequences prior to being deposited in a sequence repository (e.g. GenBank). Thus there will be no "hits" to these sequences using a primer-pair search. b) The sequence spans a region not targeted by either of the primers. c) The sequence only contains a match to one primer sequence (e.g. the sequence may have used a primer pair that is offset from the pair you've selected). d) The sequence is of low quality and finding a match is difficult. To get around this, issue RESCRIPt provides an option in which a pool of previously curated amplicon sequences (reference sequence segments) can be used to map to full-length reference sequences. Once mapped, the corresponding amplicon region will be extracted without the need of PCR primers (though PCR primer extraciton can be used to "seed" the initial amplicon-region extraction step). For more details please see: rescript extract-seq-segments.

  3. Extract from curated reference alignment: If a well curated sequence alignment is available for your gene of interest, we can use the alignment positions corresponding to the amplicon region of interest, to explicitly extract from an alignment. This approach obviates the need for any search heuristics used by the two steps outlined above. Discussed in this tutorial.

Why extract from a curated alignment?

If a trusted curated sequence alignment exists (e.g. SILVA alignment), then the amplicon region of interest can be easily and explicitly extracted from all reference sequences in the alignment that contain sequence information. Thus, mitigating potential biases of primer-pair & sequence search heuristics that are used for searching and extracting amplicon segments. If no sequence information exists for that region of the alignment, it is usually due to one of two reasons: 1) although a given reference sequence is within the curated alignment, the gene itself was not completely sequenced and is incomplete, 2) the organism to which that reference sequence belongs might not contain that region at all.

  • The caveat of using alignment positions is that it assumes you have a well-curated and trustworthy alignment, in which those alignment positions are meaningful for the group(s) under investigation! That is, your gene can be globally aligned across the majority of all known taxa, or at least the subset / group of taxa you are interested in.
  • Depending on the amplicon region under study, e.g. Fungal ITS, using primer search, or extract sequence segments, to extract your region of interest is likely your only recourse, as some marker genes are very difficult (nay impossible!) to globally align across all taxa.

Tip: extracting an amplicon region, either from PCR-primer pair search (#1) or directly from an existing alignment (#3; this tutorial) can be used to supplement the extract sequence segments approach (#2), if the goal is to find additional reference data from other sources i.e., that might not be included in an exiting reference alignment.

In this tutorial we'll use the current SILVA database as our example. But the general approach outlined here, should work for any sequence reference data for which you have a well curated, and trusted, sequence alignment.

Common hypervariable primer pairs for the 16S rRNA gene:

Mostly taken from Abellan-Schneyder et al. 2021. The table below is not an exhaustive list of available primer pairs for the SILVA alignment. Feel free to suggest other primer pairs, and their locations. Please double-check these alognment positions and let us know of any errors! :construction:

Region Primer Pair Name Reference Fwd primer (5' - 3') Rev primer (5' - 3') SILVA alignment position start and end (excluding primer region) SILVA alignment position start and end (including primer region)
V1V2 27F - 338R Salter et al. 2014 AGAGTTTGATYMTGGCTCAG GCTGCCTCCCGTAGGAGT 1043 - 6333 --
V1V3 27F - 534R Walker et al. 2015 AGAGTTTGATYMTGGCTCAG ATTACCGCGGCTGCTGG 1043 - 13126 --
V3V4 341F/357wF - 805R Herlemann et al. 2011 CCTACGGGNGGCWGCAG GACTACHVGGGTATCTAATCC 6427 - 23442 --
V3V4 341F/357wF-806R Lemons et al. 2017 CCTACGGGNGGCWGCAG GGACTACHVGGGTWTCTAAT -- --
V4 515F – 806R Caporaso et al. 2011 GTGCCAGCMGCCGCGGTAA GGACTACHVGGGTWTCTAAT 13862 - 23445 --
V4 (updated) 515F (Parada) – 806R (Apprill) Parada et al. 2016, Apprill et al. 2015 GTGYCAGCMGCCGCGGTAA GGACTACNVGGGTWTCTAAT 13862 - 23445 --
V4V5 515F (Parada) - 926R (Quince) Parada et al. 2016, Quince et al. 2011 GTGYCAGCMGCCGCGGTAA CCGYCAATTYMTTTRAGTTT 13862 - 21844 11895 - 28464
V4V5 515F - 944R Fuks et al. GTGCCAGCMGCCGCGGTAA GAATTAAACCACATGCTC 13862 - 23176 --
V6V8 939F - 1378R Lebuhn et al. 2014 GAATTGACGGGGGCCCGCACAAG CGGTGTGTACAAGGCCCGGGAACG 28601 - 41565 --
V7V9 1115F - 1492R Turner et al. 1999 CAACGAGCGCAACCCT TACGGYTACCTTGTTACGACTT 35462 - 43159 --

Note: the SILVA alignment position start and SILVA alignment position end positions are for the targeted amplicon contained within the primer pair. Thus, the extracted amplicon region will not contain the primer region / sequence! The position-start is the first base after the forward primer, and the position-end is the base just before the reverse primer.

Prepare a reference database from a SILVA alignment.

Download full SILVA alignment FASTA and import into QIIME 2

Note: We plan to enable the ability to download the full alignment files in a future release of RESCRIPt. Until then, follow the commands below for alignment based extractions of the amplicon region of interest.

Download the latest alignment from SILVA (version 138.2). Import the sequence alignment, then reverse transcribe the RNA sequence to DNA sequence.

wget https://www.arb-silva.de/fileadmin/silva_databases/current/Exports/SILVA_138.2_SSURef_NR99_tax_silva_full_align_trunc.fasta.gz

gunzip SILVA_138.2_SSURef_NR99_tax_silva_full_align_trunc.fasta.gz

qiime tools import \
  --input-path SILVA_138.2_SSURef_NR99_tax_silva_full_align_trunc.fasta  \
  --type 'FeatureData[AlignedRNASequence]' \
  --output-path silva-138-2-nr99-aln-rna.qza

qiime rescript reverse-transcribe \
    --i-rna-sequences silva-138-2-nr99-aln-rna.qza \
    --o-dna-sequences silva-138-2-nr99-aln-dna.qza

Obtain / Import SILVA Taxonomy files

If you've already generated the files from the SILVA tutorial, by running Getting SILVA data the easy way then you already have the taxonomy files you require. If not, you can run the commands under The Gritty Details, to simply download the parse the relevent files for generating QIIME 2 compatable taxonomy. That is, only the files required for qiime rescript parse-silva-taxonomy .... We'll recap these commands below to keep things easier to follow.

Import taxonomy rank file:

wget https://www.arb-silva.de/fileadmin/silva_databases/current/Exports/taxonomy/tax_slv_ssu_138.2.txt.gz

gunzip tax_slv_ssu_138.2.txt.gz

qiime tools import \
    --type 'FeatureData[SILVATaxonomy]' \
    --input-path tax_slv_ssu_138.2.txt\
    --output-path taxranks-silva-138.2-ssu-nr99.qza

Import taxonomy map file:

wget https://www.arb-silva.de/fileadmin/silva_databases/current/Exports/taxonomy/taxmap_slv_ssu_ref_nr_138.2.txt.gz

gunzip taxmap_slv_ssu_ref_nr_138.2.txt.gz

qiime tools import \
    --type 'FeatureData[SILVATaxidMap]' \
    --input-path taxmap_slv_ssu_ref_nr_138.2.txt \
    --output-path taxmap-silva-138.2-ssu-nr99.qza

Import taxonomy tree file:

wget https://www.arb-silva.de/fileadmin/silva_databases/current/Exports/taxonomy/tax_slv_ssu_138.2.tre.gz

gunzip tax_slv_ssu_138.2.tre.gz

qiime tools import \
    --type 'Phylogeny[Rooted]' \
    --input-path tax_slv_ssu_138.2.tre \
    --output-path taxtree-silva-138.2-nr99.qza

Parse SILVA taxonomy

Modify any parameters as needed. We'll use defaults for now.

qiime rescript parse-silva-taxonomy \
    --i-taxonomy-ranks taxranks-silva-138.2-ssu-nr99.qza \
    --i-taxonomy-map taxmap-silva-138.2-ssu-nr99.qza \
    --i-taxonomy-tree taxtree-silva-138.2-nr99.qza \
    --o-taxonomy silva-138.2-ssu-nr99-tax.qza

Now that we have everything we need from SILVA, let's extract our amplicon region, and then make our classifier.

Extract amplicon region from the SILVA alignment, using explicit alignment positions.

Here we'll extract the alignment region that corresponds to the V4 hyper-variable region, i.e. the "515F (Parada) – 806R (Apprill)" from the table above.

Warning: We recomend providing explicit alignment positions. Although we provide the ability to supply PCR primers for alignment region extraction, this approach is VERY inefficient for large alignments.

qiime alignment trim-alignment \
    --i-aligned-sequences silva-138-2-nr99-aln-dna.qza \
    --p-position-start 13862 \
    --p-position-end 23445 \
    --o-trimmed-sequences silva-138-2-nr99-aln-dna-V4.qza

Now that we've explicitly extracted our amplicon region of interest from the main SILVA aligment, we'll now "de-gap" (i.e. unalign) our sequences. Most downstream classifiers work on unaligned reference sequences. There are cases where sequences within the SILVA database may contain very little, or no sequence information within our targeted amplicon region (i.e. mainly consiting of alignment gaps). Also, if any reference sequences contain less than 150 nucletotide bases, after degapping, we'll remove them). Feel free to alter this value, e.g. perhaps set to 90.

qiime rescript degap-seqs \
    --i-aligned-sequences silva-138-2-nr99-aln-dna-V4.qza \
    --p-min-length 150 \
    --o-degapped-sequences silva-138-2-nr99-dna-V4.qza

Curate the database.

From here, you can optionally perform additional curation steps as described for this specific SILVA example here, or via other approaches referenced within the RESCRIPt repo.

We'll perform some simple curation below, then train our classifier.

qiime rescript dereplicate \
    --i-sequences silva-138-2-nr99-dna-V4.qza  \
    --i-taxa silva-138.2-ssu-nr99-tax.qza \
    --p-mode 'uniq' \
    --o-dereplicated-sequences silva-138-2-nr99-dna-V4-derep-uniq.qza \
    --o-dereplicated-taxa silva-138-2-nr99-tax-v4-derep-uniq.qza

qiime rescript cull-seqs \
    --i-sequences silva-138-2-nr99-dna-V4-derep-uniq.qza \
    --p-n-jobs 8 \
    --o-clean-sequences silva-138-2-nr99-dna-V4-derep-uniq-cleaned.qza

OPTIONAL: Remember, the taxonomy file can contain a superset of the IDs contined witin the sequence file. So you can skip to the training classifier step further down. However, if you'd like to keep things neat, such that your taxonomy file only contains the same IDs as in your curated sequence file, you can run:

qiime rescript filter-taxa \
    --i-taxonomy silva-138-2-nr99-tax-v4-derep-uniq.qza \
    --m-ids-to-keep-file silva-138-2-nr99-dna-V4-derep-uniq-cleaned.qza \
    --o-filtered-taxonomy silva-138-2-nr99-tax-v4-derep-uniq-cleaned.qza

Train the classifier.

We'll use the filtered taxonomy file here.

qiime feature-classifier fit-classifier-naive-bayes \
  --i-reference-reads  silva-138-2-nr99-dna-V4-derep-uniq-cleaned.qza  \
  --i-reference-taxonomy silva-138-2-nr99-tax-v4-derep-uniq-cleaned.qza \
  --o-classifier silva-138.2-ssu-nr99-classifier.qza

Note: You can also evaluate your new reference database for recent versions of QIIME 2 & RESCIPt by using qiime rescript-evaluate evaluate-fit-classifier ... and qiime rescript-evaluate evaluate-taxonomy ..., as outlined here.

Great Now you have a working clasifier! :fireworks:

If you use RESCRIPt, RESCRIPT-evaluate, or any RESCRIPt-processed data in your research, please cite the following:

  • Michael S Robeson II, Devon R O'Rourke, Benjamin D Kaehler, Michal Ziemski, Matthew R Dillon, Jeffrey T Foster, Nicholas A Bokulich. 2021. "RESCRIPt: Reproducible sequence taxonomy reference database management". PLoS Computational Biology 17 (11): e1009581.; doi: 10.1371/journal.pcbi.1009581

Be sure to also cite the reference databases which you are pulling from, see this non-exhaustive list for examples. If you use this tutorial to make your SILVA reference database, please cite the AWESOME folks over at SILVA! :

  • Maria Chuvochina, Jan Gerken, Martinique Frentrup, Yeliz Sandikci, Robin Goldmann, Heike M. Freese, Markus Göker, Johannes Sikorsji, Pablo Yarza, Christian Quast, Jörg Peplies, Frank Oliver Glöckner, Lorenz Christian Reimer. 2026. “SILVA in 2026: A Global Core Biodata Resource for rRNA within the DSMZ Digital Diversity.” Nucleic Acids Research 54 (D1): D334–41. doi: 10.1093/nar/gkaf1247.
  • Christian Quast, Elmar Pruesse, Pelin Yilmaz, Jan Gerken, Timmy Schweer, Pablo Yarza, Jörg Peplies, and Frank Oliver Glöckner. 2013. “The SILVA Ribosomal RNA Gene Database Project: Improved Data Processing and Web-Based Tools.” Nucleic Acids Research 41 (Database issue): D590-6. doi: 10.1093/nar/gks1219.
2 Likes