16S rRNA (bacteria and archaea)¶

16S is the default marker: config/config.yaml ships with amplicon.type: 16S.

What this page is, and what it is not

This page is documentation. It is not a settings file, it does not sit "above" or "below" config/config.yaml, and nothing written here is in force on your machine. Every YAML block below is an illustration of what to put in your own config.

Two real files decide how a 16S run behaves, and they do not overlap:

workflow/markers/16S.yaml — the marker pack. Facts about the marker: which databases to use, that the amplicon length is measured by in-silico PCR against SILVA, that Metaxa2 is the region extractor, that the rank model has seven slots (the last of which only the separate species database can fill). Data, not knobs. You normally never touch it.
config/config.yaml, plus any run config merged on top — the choices about your experiment: primers, expected length, truncation, taxon filter, paths.

The pack is loaded into its own namespace and is never merged into the config, so almost nothing can arrive from both. The one exception is a database path under references:, where a config entry wins over the pack default — that is how you point at a shared copy of SILVA instead of downloading your own. The full picture is in which file a setting actually comes from.

Setting up a run¶

Set type: 16S and supply the forward and reverse primer FASTAs for the region amplified. The rest of the block below reproduces the shipped defaults for the keys you are most likely to want to set on a 16S run — not every key the run touches; cutadapt, trunc_len and length_filter are covered further down and in the configuration reference. You do not need to copy any of it into your own config to get these values, since the template already provides them. The two exceptions are keep and discard: those ship commented out, so the values shown are the recommended 16S starting point rather than something already in effect.

amplicon:
  type: 16S

  primers:
    fwd: /path/to/fwd_primer.fasta
    rev: /path/to/rev_primer.fasta
    orientation: fixed            # fixed | mixed

  expected_length: auto           # in-silico PCR against SILVA
  probe_length_stat:
    16S: p95                      # which statistic of the probe distribution to use
  min_overlap: 12                 # bp of R1/R2 overlap the truncLen pair must preserve

  extraction:
    enabled: true                 # Metaxa2

  taxonomy:
    method: sintax                # rdp | sintax
    min_boot: 80                  # rdp only
    sintax_cutoff: 0.8            # sintax only
    try_rc: true                  # rdp only
    filter:
      enabled: true
      keep:    [k__Bacteria, k__Archaea]          # recommended, ships commented out
      discard: [o__Chloroplast, f__Mitochondria]  # recommended, ships commented out

Note

amplicon.probe_length_stat.16S must be present for a 16S run, whether or not expected_length is auto — MetaFlux checks the key at startup and exits with a message naming it if it is missing. The shipped config/config.yaml already has it.

Reference databases¶

Three SILVA files are involved. None of them needs to be listed in the config: the marker pack carries the filenames and download URLs, and MetaFlux fetches what the run needs on first use into references.refdb_root (default refdb), where it stays cached for later runs.

Pack symbol	Path under `refdb_root`	What it is used for	How it arrives
`silva_train`	`silva/silva_nr99_v138.2_toGenus_trainset.fa.gz`	RDP training set and the in-silico PCR substrate	download from Zenodo (rule `fetch_silva_train`)
`silva_species`	`silva/silva_v138.2_assignSpecies.fa.gz`	`addSpecies` reference (16S only)	download from Zenodo (rule `fetch_silva_species`)
`silva_sintax`	`silva/silva_nr99_v138.2_sintax.fa.gz`	SINTAX database for VSEARCH	built locally from `silva_train` (rule `convert_silva_sintax`)

The three headers look like this — the differences matter for what each step can return:

silva_train    >Bacteria;Pseudomonadota;Gammaproteobacteria;Enterobacterales;Vibrionaceae;Vibrio;
silva_species  >AB000390.1.1428 Vibrio halioticoli
silva_sintax   >1;tax=d:Bacteria,p:Pseudomonadota,c:Gammaproteobacteria,o:Enterobacterales,f:Vibrionaceae,g:Vibrio

The training set is a toGenus release: its lineages stop at genus, with no rank prefixes and no species field. That single fact explains both the addSpecies step and the rank depth of the two taxonomy methods, below.

An entry under references: in the run config is only needed to relocate a file that already exists elsewhere, for instance on a shared filesystem:

references:
  refdb_root: refdb
  silva:
    train:   /shared/dbs/silva_nr99_v138.2_toGenus_trainset.fa.gz
    species: /shared/dbs/silva_v138.2_assignSpecies.fa.gz
    sintax:  /shared/dbs/silva_nr99_v138.2_sintax.fa.gz

Amplicon length probe (`probe_mode: pcr`)¶

SILVA sequences are full-length 16S genes, so the length of the amplicon the primers actually produce has to be worked out rather than read off the database. The amplicon_probe rule does an in-silico PCR with two Cutadapt passes over silva_train:

keep the reference sequences in which the forward primer matches, trimming it and everything upstream of it (-g file:<fwd> --discard-untrimmed);
of those, keep the ones in which the reverse complement of the reverse primer matches, trimming it and everything downstream (-a file:<rev_rc> --discard-untrimmed --minimum-length 30).

Neither match is anchored, and that is the point: SILVA holds whole 16S genes, so the primers bind inside them. Anchoring to the ends would match nothing. Each pass strips the flanking gene sequence the primers would never have amplified, so what survives both is the amplicon body itself. Primer mismatches are tolerated at amplicon.cutadapt.max_error_rate (default 0.2). The rule then writes the length distribution — n_amplicons, min, q1, median, q3, p95, p99, max, mean, stdev — to a JSON file in the project-level cache:

refdb/cache/probe_16S_silva_v138.2_toGenus_<primer_hash>.json
refdb/cache/probe_16S_silva_v138.2_toGenus_<primer_hash>.amplicons.fa.gz

<primer_hash> is the first 12 hex characters of a SHA-256 over the two primer files. Different primers give a different filename, so the probe re-runs; the same primers hit the cache, which lives inside the repository and therefore survives across runs and across output directories.

Two steps consume the probe:

pick_trunclen takes expected_length from the statistic named in probe_length_stat.16S (default p95) and uses it in the overlap constraint.
dada_length_filter builds its ASV length window as [q1 − window_margin, p95 + window_margin], with window_margin defaulting to 50 bp.

The probe only runs when expected_length: auto. A fixed integer or a [min, max] range skips it entirely.

Warning

If no reference sequence carries both primers, the probe stops with an error rather than returning an empty distribution. The usual causes are a primer written in the wrong orientation for the database, or a degenerate primer that was collapsed to a single consensus sequence before being written to the FASTA — keep the IUPAC ambiguity codes, Cutadapt matches them. A low cutadapt.max_error_rate on a mismatch-tolerant primer has the same effect.

Truncation¶

16S uses the standard quality-based truncation picker: truncLen is cut just before the position where the aggregated per-base lower quartile first falls below trunc_len.q_threshold, subject to keeping the pair mergeable —

truncR1 + truncR2 >= expected_length + min_overlap

— with trunc_len.resolve_policy deciding what happens when quality and overlap disagree.

Three checks run in a fixed sequence. First, a coverage gate caps each sample at the length trunc_len.min_read_coverage_pct (95%) of its reads still reach, so the quality search cannot pick a position most reads never get to. Only then is the quality threshold applied, and only then the overlap constraint. See Amplicon length and truncation for the full picture, and Troubleshooting for what a run looks like when the cut still lands past the real read length.

Region extraction with Metaxa2¶

With extraction.enabled: true (the default), the target_extract rule runs Metaxa2 over the ASV sequences to keep only the part of each ASV that is genuinely small-subunit rRNA. Anything that flanks it — primer remnants, chimeric tails, non-ribosomal carry-over — is trimmed away, so the ASV lengths that reach the length filter are directly comparable to the probe distribution measured on SILVA.


Input	`5.dada2/seqs.fasta` and `5.dada2/seqtab_head_names.txt`, from `dada_seqtab`
Command	`metaxa2 -t all --fasta T --table T --plus T --silent T`
Output	`5.dada2/seqs_extracted.fasta`, `5.dada2/seqtab_extracted_head_names.txt`, `5.dada2/metaxa2_extraction.results.txt`
Raw tool output	kept under `5.dada2/metaxa2/` for inspection
Consumed by	`dada_length_filter`

-t all means all domains are scanned and every extracted sequence is kept. Metaxa2 is not used to decide what is bacterial: that decision is made later, by the taxon filter, on the classified lineage. ASVs Metaxa2 cannot extract are dropped from both the FASTA and the count table, and the drop is reported in the rule log.

For 16S the pre- and post-extraction length distributions usually look almost the same, because a clean 16S amplicon is SSU from end to end. Both are plotted in stats/dada2/asv_length_hist.png, which makes an unexpected difference easy to spot.

Metaxa2 reports origin twice, and only the second answer is the real one¶

Nothing in MetaFlux depends on this — extraction keeps every sequence whatever its origin, and the taxon filter makes the real call later from SILVA. It matters when you open the Metaxa2 files yourself to ask how much chloroplast and mitochondrial carry-over a run had, because the raw output under 5.dada2/metaxa2/ answers that question twice and the two answers disagree badly.

First pass — the HMM screen. Metaxa2 begins by sliding profile HMMs along each sequence to decide, cheaply, whether it is small-subunit rRNA at all and where the gene starts and stops. The job of this pass is to find the rRNA and mark its edges, not to identify it. Its verdict is the single letter in column 3 of seqs.extraction.results (B, A, E, C, M, O for bacteria, archaea, eukaryota, chloroplast, mitochondria, other) and the block of seqs.summary.txt headed "SSU rRNA sequences by preliminary origin". Metaxa2's own choice of word — preliminary — is the warning.

Second pass — the BLAST refinement. The extracted regions are then aligned against Metaxa2's reference database and classified on the alignment. This is the result worth quoting. It is what lands in seqs.taxonomy.txt, in the per-origin FASTAs (seqs.bacteria.fasta, seqs.mitochondria.fasta, and so on), and in the block of seqs.summary.txt headed "Number of SSU rRNA sequences assigned to each origin".

The distance between them is not a rounding difference. On the 16S test set, of the 231 ASVs Metaxa2 recognised as SSU rRNA:

Origin	First pass, HMM screen	Second pass, after BLAST
Bacteria	70	199
Chloroplast	56	0
Mitochondria	105	28
Uncertain	—	4

Taken from the first pass, this library reads as 45% mitochondrial and 24% chloroplast — a sample swamped by plant material, and a result that would send you looking for a problem with the primers. Taken from the second, it is 86% bacterial with no chloroplast at all.

The second reading is the correct one, and the disagreement is expected rather than a malfunction. Chloroplasts and mitochondria descend from free-living bacteria and their SSU genes still look bacterial, so the organellar HMM profiles and the bacterial ones overlap heavily. A fast profile scan cannot separate them; an alignment against full-length references can. The screen is deliberately permissive because its purpose is to avoid missing rRNA, not to name it.

These are ASV counts, not read counts

Every number above counts distinct sequences. A run can be 12% mitochondrial by ASV and a very different figure by abundance, since organellar ASVs are often numerous but individually rare, or the reverse in a root or leaf sample. To weight by reads, join the per-origin FASTA headers back to 6.taxonomy/asv_table.txt.

Taxonomy¶

amplicon.taxonomy.method picks one of two classifiers. Both read the length-filtered ASVs (5.dada2/seqs_lenfilt.fasta plus the matching seqtab_lenfilt_head_names.txt), both apply the same taxon filter, and both write the same three tables into 6.taxonomy/.

	`rdp`	`sintax`
Tool	DADA2 `assignTaxonomy`, then `addSpecies`	`vsearch --sintax`
Reference	`silva_train` + `silva_species`	`silva_sintax`
Confidence setting	`min_boot: 80`	`sintax_cutoff: 0.8`
Reverse complement	`try_rc: true`	`--strand both`, always
Deepest rank returned	Species, as a binomial	Genus
Conda environment	`workflow/envs/taxonomy.yaml`	`workflow/envs/vsearch.yaml`
Reproducibility	fully fixed by `amplicon.seed`	per-rank confidence can drift above 1 thread

The RDP path classifies the unique ASV sequences and then expands the result back so every ASV ID gets its own row — extraction can leave two ASVs with an identical trimmed sequence, and this keeps them from being classified twice.

VSEARCH SINTAX needs about 2–3 GB regardless of how many ASVs are in the table, which makes it the practical choice on large runs. Its bootstrap confidence is not byte-reproducible when it runs multithreaded: VSEARCH races several threads on one random-number stream, so a handful of per-rank confidence values near the cutoff can shift between otherwise identical runs. Setting resources.threads.assign_taxonomy: 1 removes that at the cost of speed.

The species step (16S only)¶

silva_train stops at genus, so the RDP path alone cannot name a species. 16S is the only marker whose pack declares a taxonomy_species_db, and for it the rdp branch of assign_taxonomy runs DADA2's addSpecies against silva_species — a separate, species-annotated SILVA file — immediately after assignTaxonomy. MetaFlux calls it with DADA2's own defaults, so a species is added only on an unambiguous exact match; ASVs with no match, or with more than one, keep an empty species slot.

addSpecies returns the bare epithet (subtilis), not the full binomial. The 16S pack sets prefix_style: bare, which tells the taxonomy-string builder two things: add the rank prefixes itself, and render the last slot as a Genus + species binomial by gluing the genus back on. A finished 16S lineage therefore reads:

k__Bacteria;p__Bacillota;c__Bacilli;o__Bacillales;f__Bacillaceae;g__Bacillus;s__Bacillus subtilis

Under sintax there is no species rank at all. The SINTAX database is converted from the same toGenus training set and carries only the six ranks d: to g:, so the Species column of taxon_seq_table.txt stays empty and no s__ segment appears in the taxonomy string. Species-level 16S calls require method: rdp.

Tip

Species assignment from a single 16S region is exact-match only and often absent — that is the method, not a failure of the run. Where species resolution is the point of the experiment, the two protein-coding markers are the alternative: gyrB has higher intra-genus resolution than 16S, and rpoB assigns far more ASVs to species than 16S V3–V4 does.

Note that the first rank is written k__ under both methods, even though SINTAX labels it d: internally: the prefixes are keyed to what the rank means, not to the letter the database happens to use. A keep: [k__Bacteria] token therefore behaves identically whichever method is selected.

Keeping and discarding taxa¶

There is no built-in filter default — the keep and discard lists in the shipped config are commented out on purpose, so a run that forgets them does no filtering rather than silently applying another marker's tokens. The recommended 16S starting point:

taxonomy:
  filter:
    enabled: true
    keep:    [k__Bacteria, k__Archaea]
    discard: [o__Chloroplast, f__Mitochondria]

keep runs first and defines what survives. [k__Bacteria, k__Archaea] keeps prokaryotes and drops both off-target eukaryotic hits and ASVs that stayed unclassified at kingdom level.
discard then prunes what is left. Chloroplast and mitochondrial 16S are amplified by most universal primers from any plant, algal or host material in the sample, and they are ribosomal enough to pass Metaxa2 — the filter is where they go.
Lab and kit contaminants a negative control revealed can be added to discard as further genus tokens, e.g. [o__Chloroplast, f__Mitochondria, g__Ralstonia].

Each token is matched against a whole rank segment of the taxonomy string, so g__Bacillus matches a genus named exactly that and never a substring somewhere else in the lineage. If the filter removes every ASV, MetaFlux stops with an error instead of writing an empty table — that combination nearly always means the tokens belong to a different marker. Full behaviour on Keeping and discarding taxa.

What a 16S run produces¶

The final tables are in 6.taxonomy/:

File	Rows	Contents
`asv_table.txt`	`ASV_#` IDs	per-sample counts plus a `taxonomy` column
`asv_table_seqs.txt`	ASV sequences	the same table, keyed by sequence
`taxon_seq_table.txt`	`ASV_#` IDs	one column per rank (Kingdom…Species) plus the sequence

Alongside them, stats/trunclen.json records the truncation decision, stats/dada2/asv_length_stats.json and asv_length_hist.png record the ASV length distributions before extraction, after extraction and after the length filter, and stats/read_tracking.txt shows how many reads each sample kept at every stage. The full layout is in Output files.