Configuration

Everything MetaFlux does is decided by one YAML file. There are no command-line switches for analysis parameters — the run is fully described by its config, which makes it easy to archive next to the results and re-run later.

The shipped template is config/config.yaml. It is a template: the /path/to/... entries are placeholders that must be replaced before anything runs.

How the file is loaded

The Snakefile hard-codes configfile: "config/config.yaml", so that file is always read as a base layer. A file passed with --configfile is merged on top of it, and --config key=value on the command line overrides both. A common pattern is to copy the template to config/config.local.yaml (git-ignored) and run with --configfile config/config.local.yaml, keeping the shipped template untouched.

snakemake --sdm conda --cores 16 --configfile config/config.local.yaml --dry-run

Why the taxon filter lists ship commented out

Because the template is always loaded underneath the run config, any concrete value left in it silently applies to runs that do not set their own. That is harmless for a thread count and destructive for amplicon.taxonomy.filter.keep: a 16S keep: [d__Bacteria, d__Archaea] leaking into an 18S run would discard every eukaryotic ASV. The keep and discard lines are therefore commented out in the template, so unset they resolve to [] and match nothing, which fails safe. (filter.organelles ships set: it needs no tokens and acts on 16S only, where it removes what Metaxa2 calls mitochondria or chloroplast.)

Paths may be absolute or relative. Relative paths resolve against the directory snakemake is invoked from, not against the location of the config file. Absolute paths are safest, especially on a cluster.

Each section below is tagged with the mode that consumes it: [shared], [amplicon], or [shotgun]. Keys belonging to the mode that is not running are ignored entirely.

Which file a setting actually comes from

Four kinds of file are mentioned on this site; three of them can change a run:

File	Where	Read by	Can it change a run?
These documentation pages	docs/**.md	people	No
The config template	config/config.yaml	Snakemake	Yes — always, as the base layer
A run config	any .yaml given to --configfile	Snakemake	Yes — merged on top
A marker pack	workflow/markers/<type>.yaml	Snakemake	Yes — but in its own namespace

Documentation pages set nothing

Everything on this site is prose rendered from Markdown, including the per-marker pages. Nothing under docs/ is ever read by Snakemake. The YAML blocks shown on these pages are illustrations of what to put in your own config — they are not in force anywhere, and they cannot override or be overridden by anything.

If two pages disagree, the documentation is wrong. What applies is config/config.yaml plus whatever you merge on top of it.

The config layers

Three layers, each beating the one before it:

1. config/config.yaml — the Snakefile hard-codes configfile: "config/config.yaml". That path is resolved against the directory snakemake runs in, not against the Snakefile, so it loads on every run launched from the repository root — the supported way to invoke MetaFlux. Started from somewhere else, Snakemake either aborts (Workflow defines configfile config/config.yaml but it is not present) or, if you passed --configfile, silently uses that file alone with no template underneath.
2. --configfile my_run.yaml — merged on top, key by key, at any depth.
3. --config key=value — command line, beats both. Only top-level keys can be named this way: --config mode=shotgun works, --config amplicon.dada2.merge.min_overlap=99 is rejected outright. Nested settings belong in a --configfile.

The merge in step 2 is nested rather than wholesale, which is what makes short run configs practical. A run config naming two settings changes exactly those two:

# my_run.yaml — an overlay that changes two settings
amplicon:
  dada2:
    merge:
      min_overlap: 10
output:
  out_dir: /data/my_run_out

Running with --configfile my_run.yaml gives min_overlap: 10 and the new output directory, while max_mismatch, just_concatenate, every trunc_len key and everything else still come from the template. You never have to restate a block to change one value inside it.

This is an overlay, not a complete run config: it inherits the template's /path/to/... placeholders for input.fastq_dir and amplicon.primers, which still have to be replaced before anything runs.

Only blocks are merged key by key

Key-by-key merging happens as long as what you write at each level is itself a block. Write a key with a list, a scalar, or no value at all, and it replaces that whole subtree in the template instead of merging into it.

The dangerous form is the last one, because it looks like an unfinished thought rather than an instruction:

amplicon:
  trunc_len:        # <- no value: YAML reads this as null

That silently deletes all eight trunc_len keys the template provides, and the run fails later with a KeyError from deep inside the workflow rather than a config error. Same for max_ee: [1], which replaces [2, 5] outright rather than merging element-wise. If you do not intend to change a block, leave it out entirely.

Marker packs live beside the config, not inside it

A marker pack holds facts about the marker: which database the classifier uses, which reference the amplicon-length probe measures, whether a region extractor exists, how many slots the marker's rank model has. A run config holds choices about this experiment: primers, expected length, truncation, taxon filter, paths.

The two are kept apart: a pack is read into MARKERS[type] and never merged into config, so a setting cannot come from both places. No key of the shipped packs appears under amplicon: in the config template; the only name they share is references, described below.

If you genuinely need a different pack — another database release, a different rank model — put a file of the same name in config/markers/. That is a whole-file replacement, not a per-key merge: your file supersedes the shipped pack entirely, so it must declare everything the pack needs, not just the part you wanted to change. That includes reference_values, the rules the reference check tests the reference file against before any classification; a pack for a marker of your own may write reference_values: {} to have only the rank shape checked (see the reference check). A copy made for an older MetaFlux has no reference_values and stops every amplicon run at start-up, naming the file, so after an upgrade copy the shipped pack again and re-apply your change.

Where a pack and the config do interact

They interact in two ways, handled differently.

A config entry overriding a pack happens in one place: database locations. Every pack declares a default filename and download URL for each database it owns, and the run config only has to name one when you want to relocate it — to point at a copy already on a shared filesystem rather than let MetaFlux download its own:

references:
  silva:
    species: /shared/dbs/silva_v144_assignSpecies.fa.gz

A non-empty references.* entry wins. An entry left null or empty counts as absent, so the pack default is used, without an error or a warning. The key that relocates each database is listed in Markers.

Two 16S settings do not override the pack but choose between variants it declares: references.silva.release picks the SILVA release (144 or 138.2), and amplicon.taxonomy.silva_other picks which build of the SILVA 144 files is used (other-verbatim or other-stripped in their names).

A pack overriding the config is the other direction, and the one case where the workflow does not do what the config says. It happens when a setting asks for something that does not exist for the chosen marker, typically in a config copied from another marker. There are two such cases, handled differently on purpose:

You asked for	On these markers	What happens
extraction.enabled: true	18S, gyrB, rpoB, COI	warning, the run continues with extraction off
taxonomy.method: sintax	gyrB, rpoB	error, the run stops
taxonomy.method: rdp	COI	error, the run stops

Region extraction

Metaxa2 cuts the ribosomal gene out of a 16S ASV; ITSx cuts ITS1 or ITS2 out of a fungal one. gyrB, rpoB and COI are protein-coding genes, so neither tool applies to them, and the 18S pack declares no extractor either. All four say extractor: none, and the workflow builds no extraction step at all for them.

If your config still says extraction.enabled: true, the pipeline would go looking for seqs_extracted.fasta — a file no rule produces — and Snakemake would stop with a missing-input error naming a filename you never asked for. Instead MetaFlux turns the setting off itself and tells you on stderr:

[MetaFlux] warning: marker 18S has no target-region extractor;
ignoring amplicon.extraction.enabled: true (no extraction step will run)

This one is a warning: the run continues without extraction, which is the intended result for a 16S config reused for 18S.

SINTAX classification

VSEARCH's --sintax needs a database in its own format. SILVA, UNITE and PR2 all have one, so 16S, ITS and 18S can be classified either way. gyrB's DD7RZ8 trainset and rpoB's FROGS build ship an RDP-style training set only, so those two packs leave taxonomy_sintax_db empty. Asking for sintax there stops the run at once, naming the alternative:

[MetaFlux] marker gyrB has no SINTAX reference database, so
amplicon.taxonomy.method: sintax is not available for it. Use method: rdp.

This case stops the run instead of falling back. Switching off a step that could not run changes no result; classifying with a method other than the one requested would change the results without notice.

COI is the mirror case: MIDORI2 is classified with SINTAX only, because DADA2's assignTaxonomy would need about 225 GB of memory for it at species level, so the COI pack leaves taxonomy_refdb empty and method: rdp stops the run in the same way:

[MetaFlux] marker COI has no rdp training set, so
amplicon.taxonomy.method: rdp is not available for it. Use method: sintax.

The five settings that define a run

Most of the file is tuning that can stay at its defaults. Five settings describe what the analysis actually is, the ones the box at the top of the template lists:

Setting	Choices
mode	amplicon or shotgun
amplicon.type	16S, ITS, 18S, gyrB, rpoB, COI (amplicon mode only)
shotgun.profiler	kraken2 (the default), sylph or both (shotgun mode only; see Profiler below)
input.fastq_dir	where the paired FASTQs are
output.out_dir	where results go

---

Mode selection [shared]

mode: amplicon          # amplicon | shotgun

Key	What it does	Default	Notes
mode	Chooses which set of rules is built into the workflow graph	amplicon	Case-insensitive. Anything other than amplicon or shotgun stops the run at start-up.

The mode is chosen here, never auto-detected from the data. Amplicon mode requires primer FASTAs; shotgun mode supports fetching reads from an NCBI BioProject and removing an arbitrary host genome. See choosing a mode.

The mode must match the data

Both modes accept {sample}_R1/_R2 filenames, so pointing amplicon rules at shotgun reads does not necessarily error — it just runs the wrong analysis. The mode can be flipped for a single run without editing the file: --config mode=shotgun.

---

Input [shared]

input:
  fastq_dir: /path/to/fastq_dir
  bioproject: null        # e.g. PRJNA603575  — shotgun only
  accession_list: null    # e.g. [SRR11605259, SRR11605260]

Key	What it does	Default	Notes
fastq_dir	Directory holding the paired FASTQs for every sample	/path/to/fastq_dir (placeholder)	Created if it does not exist, because BioProject mode downloads into it. Sample names come from globbing this directory.
bioproject	NCBI BioProject accession to fetch from SRA instead of using local files	null	Shotgun only. null means use whatever is already in fastq_dir.
accession_list	Restricts a BioProject fetch to named runs	null	A YAML list of run accessions. Useful when a BioProject mixes WGS with amplicon or RNA-Seq runs.

File naming. Amplicon mode requires {sample}_R1.fastq.gz / {sample}_R2.fastq.gz and stops with a clear message if it finds SRA-style {sample}_1/_2 naming instead. Shotgun mode accepts either convention and detects which one is in use. The reads themselves must be .fastq.gz or the uncompressed .fastq: fq and fq.gz pass the start-up check but are not resolved by the rules afterwards. The extension check only looks at the R1 files — a mixture there is rejected rather than half-processed, while a mismatched R2 extension is not caught, so the whole directory is best kept consistent.

Sample names are taken from the part of the filename before _R1/_1. Names containing any of * # @ % ^ / ! ? & : ; | < > or a space are rejected at start-up, since Snakemake treats those characters specially. Underscores are fine.

BioProject mode. When bioproject is set, MetaFlux queries NCBI e-utils, keeps only runs whose library layout is PAIRED, and caches the run table as {bioproject}_runinfo.csv inside fastq_dir so later invocations work offline. The download_sra rule then fetches each run with prefetch plus fasterq-dump and renames SRA's _1/_2 output to _R1/_R2.

---

Output [shared]

output:
  out_dir: /path/to/output_dir

Key	What it does	Default	Notes
out_dir	Root directory for every result, log and benchmark of the run	/path/to/output_dir (placeholder)	Logs land in out_dir/logs/, benchmark timings in out_dir/benchmarks/.

The directory layout inside out_dir differs by mode and is described in output files.

---

References

references:
  refdb_root: refdb
  phix:
    fasta: refdb/phix/phix.fna
    fetch_url: https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/819/615/GCF_000819615.1_ViralProj14015/GCF_000819615.1_ViralProj14015_genomic.fna.gz
  kraken_db: /path/to/kraken2_db
  sylph_databases: []
  sylph_taxonomy_dir: null

Key	What it does	Default	Notes
refdb_root	Cache directory for the marker reference databases MetaFlux downloads itself	refdb	Amplicon only. Shared across runs and across output directories, so a database is fetched once and reused.
phix.fasta	Local path of the PhiX genome	refdb/phix/phix.fna	Amplicon mode also builds a bowtie2 index next to it (same basename, .bt2 files).
phix.fetch_url	Where the PhiX genome is downloaded from if missing	NCBI FTP, GCF_000819615.1_ViralProj14015	Fetched once by the fetch_phix rule and gunzipped.
kraken_db	Directory of a prebuilt Kraken2/Bracken index	/path/to/kraken2_db (placeholder)	Shotgun only, user-supplied — MetaFlux never builds or downloads it. Needed when Kraken2 runs (shotgun.profiler: kraken2 or both); the run then stops immediately if the path does not exist.
sylph_databases	The sylph database files, one entry each with path and taxonomy	[]	Shotgun only, needed with shotgun.profiler: sylph or both, user-supplied. taxonomy is the sylph-tax name that goes with the file (GTDB_r232, FungiRefSeq-latest, …) or the path of your own lineage table. Checked at start-up; see sylph.
sylph_taxonomy_dir	Folder of sylph-tax lineage tables, filled once with sylph-tax download	null	Needed when a taxonomy above is a sylph-tax name. MetaFlux passes it to sylph-tax itself, so nothing in ~/.config/sylph-tax is read or written.

PhiX is used by both modes, but differently. In amplicon mode the FASTA above is fetched and indexed, and bowtie2 maps reads against it. In shotgun mode BBDuk uses phix2.fa.gz, BBTools' own corrected copy of the PhiX genome, bundled inside the BBTools conda package, so neither the download nor the index is needed there.

Marker databases are declared by the marker packs, not here

SILVA, UNITE, PR2, SILVA-Eukaryotic, DD7RZ8, FROGS, MIDORI2 and the set of bacterial COI homologues are not listed in the config. Each marker pack (workflow/markers/<type>.yaml) declares its own database filenames and download URLs, and MetaFlux fetches (and where needed converts) them on first use into refdb_root. An entry under references is needed only to choose the SILVA release for 16S (see SILVA release below) or to relocate a database — for instance to point at a copy already on a shared filesystem:

references:
  silva:
    species: /shared/dbs/silva_v144_assignSpecies.fa.gz

When such a key is present it wins; otherwise the path is refdb_root joined with the pack's declared filename. The full set of relocatable keys:

Config key	Marker	Default location under refdb_root
silva.train_raw	16S	silva/silva_nr99_v144_toGenus_trainset.fa.gz (SILVA 144 as published, seven ranks)
silva.train	16S	silva/silva_nr99_v144_toGenus_7rank_other-verbatim_trainset.fa.gz (built from train_raw: the same seven ranks, with the kingdom names written without --other; a built file, see below)
silva.species	16S	silva/silva_v144_assignSpecies.fa.gz
silva.sintax	16S	silva/silva_nr99_v144_7rank_other-verbatim_sintax.fa.gz (built from the trainset, not downloaded; see below)
unite.fasta	ITS	unite/sh_general_release_dynamic_s_all_19.02.2025.fasta
unite.sintax	ITS	unite/utax_reference_dataset_all_19.02.2025.fasta.gz
unite.uchime_its1	ITS	unite/uchime_its1.fasta
unite.uchime_its2	ITS	unite/uchime_its2.fasta
silva_euk.fasta	18S	silva_euk/silva_132.18s.dada2.fa.gz
pr2.dada2	18S	pr2/pr2_SSU_dada2.fasta.gz
pr2.utax	18S	pr2/pr2_SSU_UTAX.fasta.gz
pr2.sintax	18S	pr2/pr2_SSU_UTAX_domain-first.fasta.gz (built from pr2.utax with the domain moved to the SINTAX letter d:, not downloaded; method: sintax only; see below)
gyrb.dada2	gyrB	gyrb/train_set_gyrB_v6.fa.gz
rpob.frogs	rpoB	rpob/rpob_refseq_cc_20240707_domain-d_dada2.fa.gz (the DADA2 trainset converted from the FROGS archive, with the domain written d__; see below)
midori2.raw	COI	midori2/GB273/MIDORI2_UNIQ_NUC_SP_GB273_CO1_SINTAX.fasta.gz (MIDORI2's SINTAX file as published, extracted from its zip)
midori2.taxonomy	COI	midori2/GB273/taxonomy/nodes.dmp (MIDORI2's NCBI taxonomy snapshot; names.dmp is read from the same folder)
midori2.sintax	COI	midori2/GB273/midori2_gb273_uniq_nuc_sp_co1_metaflux_sintax.fasta.gz (built from the three files above and the MIDORI2 list of the bacterial set, not downloaded; see below)
bacterial_cox1.fasta	COI	bacterial_cox1/gtdb-r232.v1/bacterial_cox1_gtdb-r232.v1.fasta.gz (whole genes of the set of bacterial COI homologues)
bacterial_cox1.fragments	COI	bacterial_cox1/gtdb-r232.v1/bacterial_cox1_gtdb-r232.v1_fragments.fasta.gz
bacterial_cox1.midori2_list	COI	bacterial_cox1/gtdb-r232.v1/midori2_gb273_bacterial_records.tsv (the MIDORI2 records at least 90 % identical to a gene that passes checks 5.1–5.3 of the set's build, 186 in GB273, each marked removed (170, judged bacterial) or listed (16, kept); those marked removed are left out of midori2.sintax)

The silva.* paths are those of SILVA 144 with silva_other: verbatim; with silva_other: strip, other-verbatim becomes other-stripped.

Reference files MetaFlux builds itself

Five of these keys point at files MetaFlux builds rather than downloads:

Key	Built from	By rule
silva.train (SILVA 144 only)	silva.train_raw, SILVA's own download	convert_silva_train
silva.sintax	silva.train, the training set of the chosen release	convert_silva_sintax
pr2.sintax	pr2.utax, PR2's UTAX download	convert_pr2_sintax
rpob.frogs	the FROGS rpoB archive, downloaded by fetch_rpob_archive	convert_rpob_to_dada2
midori2.sintax	midori2.raw, midori2.taxonomy and bacterial_cox1.midori2_list	prepare_midori2_coi_sintax

A path given for one of them is where MetaFlux writes the file it builds. MetaFlux rebuilds such a file whenever its input or its conversion script is newer than it (a fresh clone, or an update that changes the script), or when the step's input or software changes. So never keep a file of your own at such a path; a write-protected file there stops the run as soon as MetaFlux wants to rebuild it. If you pin the input of a built file, pin the built file to a path of its own too (pr2.sintax next to pr2.utax, say): otherwise the default file is rebuilt from your input, or an older build from another input is reused. Whichever path the classifier reads, the reference check tests that very file before any classification (see the reference check).

A reference that has neither a config entry nor a pack default stops the run with a message naming both places it could be fixed.

SILVA release (16S)

references:
  silva:
    release: "138.2"     # unset = "144"

Key	What it does	Default	Notes
silva.release	Which SILVA release a 16S run classifies against	144	144 or 138.2. Each release has its own files under refdb_root, so switching never overwrites the other. SILVA 144 lineages have seven ranks, with a Kingdom between Domain and Phylum, and MetaFlux classifies on all of them: the 16S tables have the slots Domain, Kingdom, Phylum … Species. SILVA 138.2 has no kingdom, so its tables have Domain, Phylum … Species. See 16S rRNA.

---

Amplicon parameters [amplicon]

Ignored entirely when mode: shotgun. The one exception: if amplicon.primers.fwd or .rev is filled in during a shotgun run, MetaFlux prints a warning that they are being ignored — so a leftover value is visible rather than silently misleading.

Marker, region and seed

amplicon:
  type: 16S             # 16S | ITS | 18S | gyrB | rpoB | COI
  its_region: ITS2      # ITS1 | ITS2 (ignored unless type == ITS)
  seed: 42

Key	What it does	Default	Notes
type	Selects the marker pack: which databases are used, whether a region extractor runs, how the amplicon length is measured, and which taxonomic ranks the output carries	16S	Matched case-insensitively, so gyrb, GYRB and gyrB all resolve to the shipped gyrB pack. An unknown value lists the available markers and stops.
its_region	Which ITS subregion was amplified	ITS2	Only read when type: ITS. Selects the UNITE UCHIME probe reference and the ITSx region kept after extraction.
seed	Fixes the random number generator for the stochastic steps of the amplicon path	42	Passed to dada_seqtab, assign_taxonomy and, when phylogeny is enabled, phylo_tree (all three backends). Keep it a positive number: VSEARCH reads 0 as "pick a random seed", and sintax results would then change from run to run.

Each marker has its own page with the databases it uses and its quirks: 16S, ITS, 18S, gyrB, rpoB, COI.

How far the seed goes

DADA2's learnErrors(randomize = TRUE) and the bootstrap confidence of both taxonomy methods are reproducible with this seed at any thread count (for sintax from v2.5.0 on, which runs VSEARCH 2.32), except that assignTaxonomy breaks exact ties between equally good reference lineages with its own random number, which can move a bootstrap value by a point. Read order used to escape the seed — multithreaded bowtie2 phiX removal emitted the surviving pairs in a run-dependent order, so ASVs with equal total abundance could swap IDs — but bowtie2 now runs with --reorder, so the reads and ASVs come out the same in every run. See Troubleshooting.

PhiX removal

  decontamination:
    remove_phix: true

Key	What it does	Default	Notes
remove_phix	Maps reads against the PhiX genome with bowtie2 and keeps the pairs that do not align	true	Set false to skip. The 02.no_phix/ stage, its read count and the nophix Falco QC stage then do not exist, and Cutadapt reads the raw FASTQs directly.

PhiX is spiked into most Illumina runs as a low-diversity balancing control; those reads are not part of the sample and would be carried into ASV inference.

Primers

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

Key	What it does	Default	Notes
fwd	FASTA with the forward primer sequence	placeholder	Must exist — amplicon mode stops at start-up if either file is missing.
rev	FASTA with the reverse primer sequence	placeholder	Reverse complements of both are generated by the revcomp_primers rule; nothing needs to be supplied pre-complemented.
orientation	Whether R1 always carries the forward primer (fixed) or reads may come in either orientation (mixed)	fixed	mixed runs a second Cutadapt pass with the primers swapped, swaps the two mates of each pair it keeps so that its forward read is in R1, and concatenates the two passes. Both 5′ calls then need at least 10 bases of their primer instead of Cutadapt's default of 3, so that a pair is not kept by both passes by chance; this drops the pairs whose primer is found from only 3 to 9 of its bases, 0.006–0.15 % in the 16S, ITS, 18S, gyrB and rpoB test sets and at most 0.06 % in the COI mock libraries. fixed keeps the default and runs as before. For a library in which every pair is swapped, swap the two input files and keep fixed instead. See the mixed-orientation second pass. Anything other than fixed/mixed stops the run.

The primer pair also drives the probe cache. The cache key is a short hash of both primer files, so changing a primer invalidates the cached amplicon-length measurement automatically, and reusing the same primers across runs reuses it.

Amplicon length

  expected_length: auto     # int | [min, max] | "auto"
  probe_length_stat:
    16S: p95
    18S: p95
    rpoB: p95
    COI: p95
  min_overlap: 12

expected_length is the length of the amplicon the primers produce. It is not used to filter anything directly — it is the number that decides whether the chosen R1 and R2 truncation points still leave the pair able to merge.

Key	What it does	Default	Notes
expected_length	Amplicon size used for the truncation overlap check	auto	auto measures it from the reference database; an integer states it directly; a [min, max] pair states a window and the maximum is used for the overlap check.
probe_length_stat.<key>	Which statistic of the measured length distribution auto should take	p95 for 16S, 18S, rpoB and COI	One of min, q1, median, q3, p95, p99, max.
min_overlap	Bases of forward/reverse overlap the truncation picker must preserve	12	A design constraint only. The overlap actually enforced when pairs are merged is dada2.merge.min_overlap.

What auto does per marker. Markers whose reference is a full-length gene are measured by in-silico PCR: the primers are matched against the reference and the lengths of the surviving products form the distribution. Markers whose reference is already cut to amplicon length are measured directly.

Marker	How auto measures the amplicon	Probe reference	Statistic key read
16S	In-silico PCR; chloroplast, mitochondrial and eukaryote references are left out of the statistics	SILVA trainset	16S
ITS	Direct read of pre-extracted subregion lengths	UNITE UCHIME ITS1 or ITS2	not used
18S	In-silico PCR	SILVA-Eukaryotic 18S	18S
gyrB	Direct read — DD7RZ8 already ships trimmed amplicons	DD7RZ8	16S (deliberate reuse)
rpoB	In-silico PCR — FROGS ships full-length genes	FROGS RefSeq rpoB	rpoB
COI	In-silico PCR — MIDORI2 holds whole barcodes and whole genes	the converted MIDORI2 (midori2.sintax)	COI

Keep all four probe_length_stat keys

For every marker except ITS, MetaFlux checks at start-up that the key its pack names is present in probe_length_stat, whatever expected_length is set to, and exits with the missing key name if it is not. gyrB reads the 16S key — that is a deliberate reuse, not a typo. Deleting any of the four entries breaks the markers that point at it.

18S V9 must be sized manually

V4 primer pairs work with auto: in silico against SILVA-Euk, the probe reference, the median core amplicon ranges from 381 bp (TAReuk454FWD1 / TAReukREV3) to 596 bp (566F / 1200R, Hadziavdic et al. 2014). V9 (1391F/EukBr) does not, because EukBr sits at the 3′ end of the gene, where most references are truncated: only 41% of SILVA-Euk references (24% of PR2's) give a product, so auto would size the window from a biased minority. Set it manually, for example expected_length: [115, 150].

An integer must be between 50 and 10000 bp; a [min, max] pair must have exactly two positive numbers with min <= max. Out-of-range or malformed values are rejected with an explicit message rather than quietly used.

Details of the measurement and how it feeds truncation are on amplicon length and truncation.

Primer trimming (Cutadapt)

  cutadapt:
    max_error_rate: 0.2
    min_length: 50

Key	What it does	Default	Notes
max_error_rate	Cutadapt's -e — the fraction of the matched primer length that may differ	0.2	Cutadapt reads IUPAC codes in the primer natively, so a degenerate primer costs nothing from this budget; it covers sequencing errors in the first cycles and genuine primer/template mismatch, which is common with universal primers. Cutadapt's own default is 0.1 — 0.2 is deliberately more permissive because the 5′ pass runs with --discard-untrimmed, so an unmatched primer costs the whole pair. Check the per-sample match rates in stats/cutadapt/ rather than assuming. The same value is reused by the in-silico PCR probe.
min_length	Reads shorter than this after trimming are discarded, together with their mate	50	Cutadapt's --minimum-length, with --pair-filter=any.

No quality trimming happens at this stage. DADA2's error model is fitted from the whole run's quality scores, and pre-trimming by quality would remove the low-quality bases it needs in order to learn the quality-to-error relationship.

Truncation

  trunc_len:
    mode: auto                  # auto | manual
    q_threshold: 20
    q_floor: 15
    resolve_policy: raise_trunc # raise_trunc | relax_q | error
    subsample_files: 12
    manual_r1:                  # e.g. 260
    manual_r2:                  # e.g. 240
    min_read_coverage_pct: 95   # auto only

truncLen is one length per read direction; everything past it is discarded. A shorter cut removes more of the low-quality 3′ tail, which would otherwise degrade the error model; a cut that is too short leaves R1 and R2 without overlap, and a pair that cannot merge is lost.

Key	What it does	Default	Notes
mode	auto derives the cuts from observed quality; manual uses the two numbers below	auto	Manual mode skips the quality analysis completely.
q_threshold	Quality below which auto cuts	20	The picker reads the lower quartile (Q1) of the per-base quality at each cycle from Falco's report, takes the median across samples, and cuts just before the first position where that value falls below the threshold.
q_floor	Lowest threshold relax_q may descend to	15	Only used by the relax_q policy.
resolve_policy	What to do when the quality-based cuts leave too little overlap	raise_trunc	See below.
subsample_files	Carried in the template but not read by any rule currently	12	The picker uses every sample's Falco report, not a subsample.
manual_r1 / manual_r2	Explicit R1 and R2 truncation lengths	unset	Required when mode: manual — the run fails if either is missing. Ignored for ITS.
min_read_coverage_pct	auto only: a cut is never placed past the length this percentage of a sample's reads still reach	95	Read from Falco's "Sequence Length Distribution". See the note below — this is the guard against auto discarding an entire library.

What min_read_coverage_pct is protecting against

Falco reports quality for every position up to the longest read in a file, however few reads reach it. filterAndTrim discards any read shorter than truncLen instead of padding it, so a cut placed where only a few reads reach discards almost the whole library. On a run whose quality never drops below q_threshold, nothing else stops the cut landing there.

Before quality is looked at, each sample's quality table is therefore cut off at the length min_read_coverage_pct of its reads still reach. On the 18S test dataset this moves the R1 ceiling from 288 bp to 271 bp. At 288 / 239, the cut MetaFlux v2.1.0 and earlier chose there, only 9 of the 908,768 primer-trimmed read pairs are long enough to pass; at 271 / 239, 90.9% of the raw read pairs reach the final table.

The value is a guarantee about reads: raising it asks for a position that more reads reach, which gives a shorter cut that discards fewer reads; lowering it gives a longer cut that discards more. Between about 80 and 99 the cut barely moves, because read lengths cluster. 100 does not switch the guard off: it requires every read to reach the cut, which puts the ceiling at the shortest read in the file. No value disables the guard. logs/pick_trunclen.log prints each sample's longest read next to the length the requested percentage reaches, and emits a NOTE when the two differ by 10 bp or more.

The overlap constraint is truncR1 + truncR2 >= expected_length + min_overlap. When the quality-based cuts break it:

Policy	Behaviour
raise_trunc	Extends the cuts beyond the quality drop one base at a time, each time on the read with the higher aggregated lower quartile at the position it would gain (ties go to R1). Stops with an error if R1 and R2 together cannot gain the missing bases before reaching their ceilings.
relax_q	Lowers the quality threshold by 1 at a time, down to q_floor, re-cutting until the overlap holds. Errors if the floor is reached without success.
error	Stops immediately and reports the shortfall in bp.

The decision, including which policy resolved it and how much slack is left, is written to stats/trunclen.json.

ITS ignores truncLen

ITS amplicon length varies genuinely between taxa, so a fixed cut would throw away real short amplicons. For ITS the picker still computes and records the quality-based cuts, but dada_filter overrides truncLen to c(0, 0). The 3′ tail is not left untrimmed: dada2.trunc_q still trims each read adaptively, with dada2.max_ee and the per-read minimum length as the quality gate. This override also beats mode: manual, so manual_r1/manual_r2 have no effect for ITS.

On very clean runs the auto cut is the read-length limit

If quality never drops below q_threshold, the cut is placed at the length reached by min_read_coverage_pct (95%) of the reads in the strictest sample. Without this limit, the cut could fall above most trimmed reads, and almost every read would be discarded at the filtering step; see troubleshooting. On binned data, the cut can also be placed at this limit when quality declines towards the end of the reads; see binned quality scores.

ASV length filter

  length_filter:
    mode: auto          # auto | manual
    window_margin: 50
    range: null         # [min, max]

This filter runs on the ASVs, after target extraction, so the lengths it sees are directly comparable to the probe distribution the window was built from.

Key	What it does	Default	Notes
mode	auto builds the window from the probe; manual requires an explicit range	auto
window_margin	Bases added below and above the probe window	50	The auto window is [probe_q1 − margin, probe_p95 + margin].
range	Explicit [min, max] in bp	null	When set it wins, whatever mode says.

When mode: auto runs without a probe (because expected_length was given manually), the window falls back to the resolved expected length ±15%.

The filter writes stats/dada2/asv_length_stats.json — pre-extraction, filter-source and kept distributions side by side — plus a histogram with the window drawn on it. If the window removes every ASV the run stops rather than writing an empty table.

ITS auto needs an auto expected length

type: ITS with length_filter.mode: auto also requires expected_length: 'auto', since ITS never resolves an expected length through the truncation picker and the fallback would have nothing to size from. MetaFlux catches this combination at start-up, before the expensive DADA2 run, and the message offers two fixes: set expected_length: 'auto', or switch to length_filter.mode: manual with an explicit range.

Target-region extraction

  extraction:
    enabled: true

Key	What it does	Default	Notes
enabled	Runs the marker's region extractor over the ASVs, trimming each to the target region and dropping ASVs the extractor cannot place	true	Metaxa2 for 16S, ITSx for ITS. 18S, gyrB, rpoB and COI have no extractor. For 16S the organelle filter (taxonomy.filter.organelles) works from Metaxa2's call, so with false it does nothing and metaxa2_origin in asv_table_unfiltered.txt is NA.

Markers with no extractor (18S, gyrB, rpoB, COI) should be run with enabled: false. Left true, MetaFlux prints a warning and forces it off rather than asking the workflow graph for a file no rule produces. The ASV length filter still runs either way.

Taxonomy

  taxonomy:
    method: sintax        # rdp | sintax
    min_boot: 80          # rdp only
    sintax_cutoff: 0.8    # sintax only
    try_rc: true          # rdp only
    silva_other: verbatim # 16S on SILVA 144 only: verbatim | strip

Key	What it does	Default	Notes
method	Which classifier assigns taxonomy	sintax	rdp is DADA2's assignTaxonomy (plus addSpecies for 16S); sintax is vsearch --sintax, which uses roughly 2–3 GB regardless of ASV count (about 8 GB against MIDORI2 for COI). Both paths apply the same taxon filter and write the same tables.
min_boot	RDP bootstrap confidence a rank must reach to be reported	80	Read only by the rdp path. Range 0–100.
sintax_cutoff	SINTAX confidence a rank must reach to be reported	0.8	Read only by the sintax path. Range 0–1; 0.8 is roughly equivalent to min_boot: 80.
try_rc	Also tries the reverse complement of each ASV during classification	true	Read only by the rdp path. The sintax path always searches both strands.
silva_other	What to do with SILVA 144's X--other names (a clade that does not contain the type of taxon X but belongs to it)	verbatim	verbatim keeps the names; the filter and addSpecies then treat X--other as X. strip removes the suffix before training, so the two become one taxon. The kingdom is always written without the suffix, whichever value is set. No effect with SILVA 138.2 or other markers. See 16S rRNA.

gyrB and rpoB are rdp-only, COI is sintax-only

The DD7RZ8 and FROGS releases ship a DADA2 trainset with no SINTAX build, so those packs declare no SINTAX reference. method: sintax for either marker is refused at start-up with a message naming the marker, rather than failing later on a missing path. The COI pack declares no rdp trainset (one trained on MIDORI2 would need about 225 GB of memory at species level), so method: rdp is refused for COI in the same way.

Taxon filter

    filter:
      enabled: true
      keep:    #[d__Bacteria, d__Archaea]
      discard: #[o__Chloroplast, f__Mitochondria]
      discard_min_support: 50
      organelles: metaxa2   # metaxa2 | off; acts on 16S only

Key	What it does	Default	Notes
enabled	Master switch	true	false keeps every ASV regardless of the lists, and switches organelles off too. For COI the frame check and the off-target guard have their own switches (amplicon.coi) and are not affected by this one.
keep	Tokens defining what survives	unset → []	An ASV is kept only if one of its ranks matches a token. An empty list does nothing.
discard	Tokens pruning what keep left	unset → []	An ASV is removed if any of its ranks matches a token. Runs after keep.
discard_min_support	Best-hit check for discard, in percent	50	An ASV is also removed when the classifier's best hit carries a discard token at that token's rank with at least this support (rdp bootstrap; SINTAX confidence ×100), even if the rank is below min_boot / sintax_cutoff and missing from the reported lineage. Catches mitochondrial 16S that is reported only to phylum or class. A number from 0 to 100. null (or false / off / no) switches it off, which gives the v2.4.0 behaviour. Anything else stops the run at start-up: a value outside 0–100, text that is not a number, and also true / on / yes, which YAML reads as a yes/no value rather than a number. It does nothing while discard is empty. See the best-hit check.
organelles	Also removes the ASVs that Metaxa2 calls mitochondria or chloroplast	metaxa2	16S only, and only with extraction.enabled: true, since Metaxa2 has to have run; for every other marker it is inactive and the run prints one info line. Runs after keep, discard and the best-hit check, and catches plant mitochondria the classifier reports as bacteria (e.g. Enterobacterales), which no token can match. Metaxa2's other calls (domains, eukaryota, uncertain) remove nothing. off (or null) switches it off; any other value stops the run at start-up. Every ASV, organelles included, is still listed in 06.taxonomy/asv_table_unfiltered.txt. See the organelle filter.

Each token is matched against a whole segment of the taxonomy string (d__Bacteria, o__Chloroplast, g__Ralstonia), so it targets exactly one rank and never matches loosely somewhere else in the lineage. d__ is always the domain and k__ always the kingdom, where the marker's reference has one: on 16S with SILVA 144, and on ITS, whose top rank is the kingdom. There is no built-in default: whatever is listed is exactly what is applied.

The tokens are checked when the run starts, before any step, also with enabled: false. A token that can never match this run's taxonomy strings stops the run, and the message names the token to write instead where there is one: a prefix the strings never have (k__ on SILVA 138.2, d__ on ITS, s__ on 16S under sintax), a domain under k__ (k__Bacteria, written d__Bacteria), or a name that the marker's reference cannot hold in that slot (k__Fungi on 16S). See when MetaFlux refuses a token at start-up.

Marker	Suggested keep	Suggested discard
16S	[d__Bacteria, d__Archaea]	[o__Chloroplast, f__Mitochondria]
ITS	[k__Fungi]	[]
18S	[d__Eukaryota] (PR2 ranks) — note this also drops PR2's organelle-tagged domains (d__Eukaryota:plas, :mito, :nucl, :apic, :chro) and the few ASVs left unclassified at Domain; see 18S	[]
gyrB	no recommended default — paralogs are removed with discard: [other] on the DD7RZ8 gene-tag rank, the same keep/discard mechanism as above, not a separate filter
rpoB	no recommended default — FROGS RefSeq has no equivalent contamination to remove
COI	[k__Metazoa] for animal surveys; empty to see everything MIDORI2 can name	[] — bacterial COI homologues are moved by the off-target guard (amplicon.coi), not by a token

Contaminants a blank extraction control reveals go in discard, e.g. discard: [o__Chloroplast, f__Mitochondria, g__Ralstonia, g__Bradyrhizobium]. A bare string is accepted where a list is expected — keep: k__Fungi is treated as keep: [k__Fungi] — so a single token cannot be split into characters by accident. The full treatment is on keeping and discarding taxa.

The old regex keys are refused

Configs written before the token lists used include_pattern / exclude_pattern regexes. Those keys now stop the run with instructions, rather than being ignored and leaving the filter silently disabled.

COI screens

  coi:
    frame_check: true
    length_tolerance: 9
    offtarget_guard: true

Used only when type: COI (the value of length_tolerance is checked at start-up for every marker). Rule coi_screens measures every ASV for both screens and writes 06.taxonomy/coi_screens.tsv; the taxonomy step applies them before the taxon filter. A screen switched off still measures, so the table says what it would have done. Its verdict column applies only the screens that are on, the frame check first (an ASV that fails both is frame). For a screen that is off, read frame_ok and length_ok (frame check) or compare B with E (off-target guard). See COI.

Key	What it does	Default	Notes
frame_check	Removes ASVs with TAA or TAG in the run's reading frame (or in the frame an offset of the ASV's start implies), or whose length is too far from the run's reference length (pseudogenes, NUMTs, artefacts)	true	The run's frame and reference length are taken from the ASVs themselves (the frame open in most reads; the most frequent length among open ASVs). Removed ASVs stay in asv_table_unfiltered.txt as removed_frame.
length_tolerance	Largest length difference from the run's reference length, in bases, that the frame check accepts	9	Three codons. Must be a whole number of bases, a multiple of 3 and not negative (9.0 and "9" are read as 9); anything else, a fraction, true/false or other text included, stops the run at start-up.
offtarget_guard	Moves ASVs closer to the set of bacterial COI homologues (B) than to MIDORI2 (E) out of the main tables, into 06.taxonomy/bacterial_cox1_asvs.tsv	true	Wolbachia and other bacterial COI homologues. Moved ASVs stay in asv_table_unfiltered.txt as moved_bacterial_cox1.

DADA2 — tier 1, operational parameters

  dada2:
    pool: false             # false | "pseudo" | true
    max_ee: [2, 5]
    trunc_q: 2
    filter:
      min_len: 20
      min_len_stat: null
      max_len: null
    learn_errors:
      nbases: 1e8
      max_consist: 999

Key	What it does	Default	Notes
pool	How samples are pooled during ASV inference	false	false treats each sample independently (fastest); "pseudo" runs a second pass informed by ASVs seen across samples; true pools all samples (most sensitive to rare variants, most expensive). Strings are read case-insensitively.
max_ee	filterAndTrim's maxEE — the expected-error ceiling per read, as [R1, R2]	[2, 5]	R2 is looser because reverse reads degrade faster on Illumina chemistry.
trunc_q	filterAndTrim's truncQ — trims each read at the first base at or below this quality	2	The DADA2 default. This is the adaptive 3′ trim, distinct from the fixed truncLen cut. It is applied before truncLen; see the note below.
filter.min_len	Per-read length floor after trimming	20	DADA2's own filterAndTrim default. Used for every marker; for ITS it is the fallback when min_len_stat is unset.
filter.min_len_stat	Optionally derives min_len from the probe distribution instead	null	ITS only, opt-in. Left null, every marker uses filter.min_len. Setting a statistic (e.g. q1) puts the floor near 175 bp, above ~24% of the UNITE reference — see ITS. The other markers never take this path: their probe measures the full amplicon, which exceeds read length.
filter.max_len	Per-read length ceiling	null	null means no upper limit.
learn_errors.nbases	Bases sampled to train the error model	1e8	The DADA2 default. Reads are drawn in random order, seeded by amplicon.seed.
learn_errors.max_consist	Most rounds of error learning per read direction: DADA2 denoises the reads with its current error rates, estimates the rates again from the result, and stops when a round gives exactly the rates of an earlier one	999	Not unlimited. Most models settle within 10 rounds and the slowest seen took about 100. A model that has not settled by the last round is used as it is, and logs/dada_seqtab.log says Self-consistency loop terminated before convergence. On a large library that never settles, 999 rounds can take days (see Troubleshooting). DADA2's own default is 10.

truncQ runs before truncLen, and the leftovers are discarded rather than kept short

This is DADA2's own ordering inside filterAndTrim, not something MetaFlux changes, but it produces a failure that is easy to misdiagnose. truncQ trims each read at its first low-quality base; truncLen is then applied, and any read that truncQ has already shortened below truncLen is dropped, not passed through at its shorter length.

So a run that loses an unexpected share of reads at dada_filter may be losing them to truncQ, while the number that looks wrong is truncLen. Before lowering truncLen, check the filtering columns in stats/read_tracking.txt and consider whether a lower trunc_q is the real cause — particularly on a run whose 3′ quality degrades early.

DADA2 — tier 2, advanced tuning

    dada:
      omega_a: 1e-40
    merge:
      min_overlap: 12
      max_mismatch: 0
      just_concatenate: false
      trim_overhang: true
    chimera:
      method: consensus
      min_fold_parent: 4
      allow_one_off: false

Key	What it does	Default	Notes
dada.omega_a	OMEGA_A — how strong the evidence must be before a sequence is called a new ASV rather than an error of an existing one	1e-40	The DADA2 default. A lower value is stricter; a higher one yields more ASVs.
merge.min_overlap	mergePairs minimum overlap, in bp, enforced at runtime	12	DADA2's default. This is the one that actually decides whether a pair merges; amplicon.min_overlap is the design constraint used earlier when choosing truncation lengths.
merge.max_mismatch	Mismatches tolerated inside the overlap	0	DADA2's default, and deliberate: after denoising, a mismatch in the overlap is evidence the pair should not merge rather than a sequencing error.
merge.just_concatenate	Joins R1 and R2 with a run of Ns instead of merging on overlap	false	A safety valve for ITS amplicons too long to overlap at all. The joined sequence is R1 + R2 + spacer, not the amplicon length, so length_filter.mode: auto is no longer sizing a meaningful window — switch to manual with an explicit range sized for the concatenated length.
merge.trim_overhang	Trims bases that read past the far end of the amplicon	true	Safe for 16S; necessary for short ITS amplicons where the read runs off the end into primer sequence.
chimera.method	removeBimeraDenovo method	consensus	Chimeras are called per sample and the verdicts pooled.
chimera.min_fold_parent	How much more abundant a parent must be than its suspected chimera	4	DADA2's minFoldParentOverAbundance, raised above the package default. Higher is more permissive: parents must be more abundant to qualify, so fewer sequences are called chimeric and more ASVs survive. See the tip on the chimera settings below the table.
chimera.allow_one_off	Also flags sequences one mismatch or indel away from an exact chimera	false (DADA2's own default)	With false, only exact chimeras of two parents are flagged. true also flags sequences one mismatch or indel away from an exact chimera, so more sequences are removed.

What the min_overlap and max_mismatch defaults cost in practice

Looser values, such as 10 and 2, merge more pairs, but the extra pairs mostly turn into distinct, low-abundance ASVs, so the stricter defaults cost far more rare ASVs than reads. Filtering and denoising are not affected; the whole difference appears at the merge step. The measured cost for each marker is on the Test datasets page. Loosening them recovers those ASVs, at the cost of accepting overlap disagreements that DADA2 treats as evidence against merging.

ITS loses much more. ITS2 length varies widely between fungi, so with 2×300 reads the longer amplicons overlap by only a little — exactly the pairs a 10 → 12 bp requirement drops, and their short overlap sits in the low-quality 3′ tails where max_mismatch: 0 then removes more. On an ITS run the two settings together can cost a large share of the merged reads, so for ITS on a variable community min_overlap: 10 with max_mismatch: 2 is a defensible setting.

Where to look before changing the chimera settings

min_fold_parent: 4 is a deliberate departure from the DADA2 default, in the direction of removing fewer sequences. The reasoning, and the values DADA2's authors and users have discussed for pooled and amplicon-specific data (typically 4 to 8), are in these DADA2 issues; the right value depends on the marker and on how samples were pooled:

• DADA2 issue: Questions about chimera filtering
• DADA2 issue: High loss of reads following chimera removal
• DADA2 issue: Remove chimeras, consensus vs. pooled

Before changing these settings, check that primers were fully removed: residual primer sequence is a common cause of apparent chimeras.

Phylogeny (16S)

amplicon.phylogeny. Off by default, and 16S only — enabling it for any other marker is forced off with a warning at start-up. Builds a tree from the contaminant-filtered ASVs in 06.taxonomy/ for phylogenetic-diversity analysis you perform yourself; MetaFlux computes no diversity statistics. The full discussion, including why the other markers are excluded and how to use the tree in R, is on the Phylogeny page.

Key	What it does	Default	Notes
enabled	Master switch	false	With it off, the workflow has no phylogeny steps and builds no phylogeny conda environment.
backend	Which tree program runs	iqtree	iqtree (full ML, model selection available) | fasttree (approximate ML, much faster, no model selection) | raxml-ng (full ML alternative). All three produce the same canonical outputs.
aligner_strategy	Which MAFFT strategy	auto	auto: <200 ASVs → L-INS-i, ≥200 → FFT-NS-2, at MAFFT's own documented ceiling for its accuracy strategies. linsi / fftns2 force one — useful to keep runs of different sizes comparable. MAFFT's own --auto is never used: its thresholds are unpublished and size-dependent. The resolved choice is recorded in phylogeny.params.json.
model	Substitution model	auto	auto → GTR+F+G4 (iqtree), GTR+Γ (fasttree, raxml-ng). Pinned rather than searched because IQ-TREE's own -m MFP default makes runtime data-dependent and re-decides on every run. Model selection is still reachable: MFP (IQ-TREE ModelFinder) or DNA (RAxML-NG MOOSE), and what it picks is recorded. FastTree has no model selection — only gtr or jc are accepted, anything else is a hard error rather than a silently ignored value.
support	Branch support values	false	Off because neither UniFrac nor Faith's PD uses them and they dominate runtime. true → -B 1000 --alrt 1000 (iqtree) or FastTree's SH-like local supports (0–1, not bootstrap). Not available for raxml-ng: forced off with a warning. Raise mem_mb.phylo_tree if you enable it.
extra_args.{mafft,iqtree,fasttree,raxml_ng}	Verbatim passthrough	""	Appended to the tool's command line and recorded in phylogeny.params.json. Not validated — this is the one route by which machine-dependent behaviour MetaFlux otherwise blocks (-T AUTO, --threads auto) can re-enter a run. An unknown key here stops the run at start-up, so a typo cannot silently do nothing.

The seed is amplicon.seed, shared with the rest of the amplicon path — there is no separate phylogeny seed.

There is deliberately no rooting option: the tree is exported unrooted, only. You will prune it in R (decontam, abundance filtering) before computing any metric, and a root placed on the full ASV set stops being valid at the first pruned tip. A root: key is rejected at start-up with a pointer to the docs rather than silently ignored.

The tree is reproducible at the shipped settings

Two runs on the same data give byte-identical trees when three settings are kept: --reorder on the phiX step, an alignment built from a sequence-sorted FASTA, and phylo_tree: 1. With more phylo_tree threads, two runs on the same alignment can end on different trees of about equal likelihood. The taxonomy step that decides which ASVs reach the tree needs no such setting: with a fixed amplicon.seed it gives the same result at any thread count (sintax from v2.5.0 on; rdp apart from DADA2's rare random tie-breaks, see the seed note above). Details in the Reproducibility section of the Phylogeny page.

---

Shotgun parameters [shotgun]

Ignored entirely when mode: amplicon.

Profiler

shotgun:
  profiler: kraken2

Key	What it does	Default	Notes
profiler	Which profiler turns the cleaned reads into taxon tables	kraken2	kraken2: Kraken2 + Bracken, 03.abundance/otu_table.tsv. sylph: sylph + sylph-tax instead, 03.abundance/sylph_table.tsv, read pairs estimated from genome coverage plus an unassigned row. both: the two side by side on the same reads, never merged. Anything else stops the run at start-up.

The PhiX, host and fastp steps below are the same for all three. Under sylph, extract_taxa (it needs Kraken2's per-read output) and kmer_evidence (it filters the Kraken2 report before Bracken) are switched off with a warning, and the kraken and bracken settings are not used. taxonomy_filter acts on otu_table.tsv only; with sylph or both it warns that the sylph tables are written unfiltered. See sylph.

Decontamination

shotgun:
  decontamination:
    remove_phix: true
    host_genomes: []
    host_min_id: 0.95

Key	What it does	Default	Notes
remove_phix	BBDuk k-mer matching against the bundled phix2.fa.gz (BBTools' corrected PhiX genome), at k=31 hdist=1	true	No download or index build is involved. Keep it on unless the lab is known to have skipped the spike-in.
host_genomes	YAML list of local paths and/or URLs to host FASTA(.gz) references	[]	[] skips host removal. Multiple entries are concatenated and indexed once, so one BBMap index serves every sample.
host_min_id	BBMap minid — the alignment identity above which a read is called host-derived	0.95	Reads that do not map are kept.

Host removal is alignment-based rather than k-mer based, which is more sensitive to host reads that diverge from the reference assembly — the approach the BBTools authors recommend.

Use a masked reference for human

An unmasked human genome removes real microbial reads that happen to be homologous to ribosomal or low-complexity human sequence. The masked reference on Zenodo at https://zenodo.org/records/4116107 (Handley 2020, human_virus_masked.fasta.gz, 889 MB) is a masked hg19 — ribosomal RNA, plant, animal and fungal homology and low-entropy sequence removed — with a further layer of viral masking added. The first set of masks is what stops real microbial reads being discarded as host; the viral layer additionally keeps human regions homologous to viruses out of the filter, which is the right behaviour when the viral fraction is of interest. Because the file is roughly 900 MB, download it once to a local path and point host_genomes at it rather than re-fetching the URL on every fresh run.

BBDuk and BBMap both need correctly paired files

Both tools read R1 and R2 in lockstep and require the two files to hold the same number of records in the same order. Read length is not the issue: ragged reads, and mates of one pair with different lengths, are processed normally. With the bbmap release MetaFlux pins, a crash reading There appear to be different numbers of reads in the paired input files (BBDuk or BBMap) means the pairing is broken — usually because reads were filtered before deposition and orphaned mates were left in one file. Repair the FASTQs with the tool shipped in the same conda package rather than skipping the step.

repair.sh in1=R1.fastq.gz in2=R2.fastq.gz \
    out1=fixed_R1.fastq.gz out2=fixed_R2.fastq.gz outs=singletons.fastq.gz

A crash reading java.lang.AssertionError: List size mismatch: N vs M is not a pairing problem. It comes from newer bbmap releases (39.76 to 40.02), whose BBDuk prints it on correctly paired reads too; run with MetaFlux's pinned environment (--sdm conda) instead of repairing the files.

See troubleshooting and decontamination.

Adapter and quality trimming

  fastp:
    min_read_length: 100

Key	What it does	Default	Notes
min_read_length	fastp's --length_required — reads shorter than this after trimming are dropped	100	Comfortable for 2×150 bp reads; use around 70 for HiSeq 2×100.

fastp runs with --detect_adapter_for_pe, plus --cut_front and --cut_right sliding-window trimming from both ends. Its JSON report is read later by Bracken to pick a read-length model, and by the read-tracking table.

Kraken2

  kraken:
    confidence: 0.15
    hit_groups: 3
    memory_mapping: true
    threads: null

Key	What it does	Default	Notes
confidence	Fraction of a read's k-mers that must fall inside the clade of the taxon it is assigned to	0.15	Higher is stricter. Kraken2 does not simply reject a read that falls short — it moves the label up the taxonomy until the score clears the threshold, and calls the read unclassified only if even root fails. Raising it therefore mostly pushes assignments to coarser ranks, leaving more for Bracken to redistribute, and only secondarily increases the unclassified fraction. Wright et al. (2023) report two optima on simulated data: 0.60 by mean F1, and 0.15 by L1 distance and by recall/precision over classified reads. The default of 0.15 favours accurate relative abundances; raise it towards 0.60 to favour F1.
hit_groups	--minimum-hit-groups — distinct k-mer hit groups a read needs to keep its classification	3
memory_mapping	Pages the database from disk on demand instead of loading it into RAM	true	See below.
threads	Pins Kraken2's thread count	null	null falls back to resources.threads.kraken2.

Kraken2 is also run with --use-names, --paired, --gzip-compressed, --report-minimizer-data and --classified-out, so classified read pairs are kept alongside the report.

Memory mapping decides how much RAM a Kraken2 job asks for

Kraken2's memory request is computed at start-up, not taken from the resources.mem_mb.kraken2 entry, which is not read.

• memory_mapping: true — the kernel pages the database on demand and shares those pages between concurrent Kraken2 processes through the OS file cache, so running several samples in parallel does not multiply RAM use. Each job is budgeted a fixed 20000 MB of private workspace.
• memory_mapping: false — every process loads the whole hash privately. The budget becomes the size of hash.k2d plus 10% headroom, so it self-adjusts to the database in use (a 16 GB PlusPF, a ~100 GB PlusPF, a 200+ GB core_nt). If hash.k2d cannot be read yet — while the database is still downloading, say — a conservative 110000 MB is used instead. (With memory mapping on the fixed 20000 MB applies either way, since the index size is never consulted.)

Bracken

  bracken:
    tax_lev: S
    threshold: auto
    threshold_alpha: 5.0e-05
    threshold_min: 10

Key	What it does	Default	Notes
tax_lev	Rank at which abundance is re-estimated: P, C, O, F, G or S	S	Read only when Kraken2 runs (profiler: kraken2 or both). Drop to G for a shallow database, where species-level calls are not supported by the reference content. The OTU table then has one row per taxon at this rank, and the lineage below it is left empty (…;g__Bacteroides;s__). D and K stop the run at start-up: at those ranks Bracken counts the reads of viruses (D) or of protists (K) towards a neighbouring domain or kingdom, and on a GTDB database it finds no kingdom at all. For domain or kingdom totals, add up a lower-rank table by its d__ or k__ slot; see Troubleshooting. P, C, O and F print a warning at start-up: at those ranks Bracken can add the reads of taxa that lack the rank to a row they do not belong to (same Troubleshooting entry), so S or G is recommended. Up to v2.4.0, any value other than S gave an empty or nearly empty table; see Troubleshooting.
threshold	Minimum Kraken2 read count a taxon needs at tax_lev to enter Bracken's estimate	auto	Applied to Kraken2's own count for the taxon, before Bracken re-estimates anything. Taxa below it are removed from the model and their reads discarded rather than redistributed (Bracken reports these as "reads discarded"). auto computes it per sample from that sample's own classified read count, at run time, rather than using one fixed number for every sample regardless of depth; see Choosing confidence and threshold for the benchmark behind this. On the D6331 mock, the fixed default of 10 costs 0.427 F1 against the best threshold. Set an integer instead of auto to bypass this and use a fixed floor as Bracken does natively.
threshold_alpha	Scaling coefficient used only when threshold: auto	5.0e-05	threshold = max(threshold_min, threshold_alpha × classified_read_pairs). Set from the benchmark data. On the CAMI II marine, rhizosphere and strain-madness samples the optimum is a broad plateau around this value; on the ZymoBIOMICS D6331 standard the best value is several times higher. See Choosing confidence and threshold before changing it.
threshold_min	Floor for the auto-computed threshold	10	Prevents a near-zero threshold on a very shallow or heavily filtered sample. Also Bracken's own historic default.

Bracken redistributes the reads Kraken2 assigned to internal nodes down to the target rank. Bracken databases ship k-mer distributions for several read lengths (50, 75, 100, 150, 200, 250, 300); MetaFlux reads the post-trimming mean read length out of the fastp JSON and picks the closest one, because applying a 150 bp distribution to 100 bp reads biases the redistribution.

Species-level k-mer evidence

  kmer_evidence:
    enabled: false
    min_distinct_minimizers: 300
    min_reads: 0
    never_prune_above_fraction: 0.01

Key	What it does	Default	Notes
enabled	Whether Bracken reads the gated report (the Kraken2 report without the taxa that have too little k-mer evidence) instead of Kraken2's own	false	The evidence table, the gated report and the plot are written either way, so false changes nothing in the pipeline. Switched off with a warning under profiler: sylph.
min_distinct_minimizers	The floor: how many distinct minimizers a taxon at tax_lev needs to stay	300	A whole number is one floor for every sample of the run; 300 is the published 1,000 unique k-mers converted to Kraken2's minimizers and rounded down. auto gives each sample its own floor, from its depth: (classified read pairs − Homo sapiens read pairs) / 1,000, at most 6,500, worked out after Kraken2 and written to logs/kmer_evidence/{sample}.log. auto is opt-in: its gain on MetaFlux's benchmark rests mostly on one real mock community; it changes the table much more than 300 does, and read pairs leave the table (2.7–15.1 % on the test dataset, against at most 0.09 % at 300); and it assumes paired 2×150 bp reads. Any other value stops the run at start-up. See Scaling the floor with depth.
min_reads	A clade-read floor, applied together with the one above: a taxon must clear both	0	0 judges on k-mer evidence alone; 200 reproduces aMeta's pair of filters.
never_prune_above_fraction	Safety rail: a taxon holding more than this share of the classified reads is never removed	0.01	The share is Kraken2's clade reads over all classified reads, not the taxon's share of the final Bracken table, so the rail rarely fires.

The gate sits between Kraken2 and Bracken: it filters Bracken's input and does not replace Bracken, so with enabled: true the corrected counts reach otu_table.tsv. A removed taxon's reads go to its parent, for Bracken to share among the relatives that passed. See Species-level k-mer evidence.

OTU-table taxon filter

  taxonomy_filter:
    enabled: false
    keep:    []
    discard: []

Key	What it does	Default	Notes
enabled	Master switch	false	Off by default: metagenomic sample types vary too much for a universal default to be safe.
keep	Tokens defining what survives	[]	e.g. [d__Bacteria, d__Archaea] for prokaryotes only; this drops viruses too, because a virus carries its realm in d__ (d__Duplodnaviria), or nothing there.
discard	Tokens pruning what keep left	[]	e.g. [g__Homo] as an explicit host safety net.

This is the shotgun counterpart of amplicon.taxonomy.filter and uses the same rank-aware token matching: keep first, then discard, empty list means that direction does nothing. It is applied to the final OTU table after Bracken, so it changes which OTUs appear in the table and never how reads were classified. A token must equal one whole slot of the table's eight, d__ k__ p__ c__ o__ f__ g__ s__.

When Kraken2 runs, the tokens are checked at start-up, also with enabled: false, against those slots and against the Kraken2 database: a wrong prefix, a domain under k__, Viruses under any prefix, a virus realm under k__, and (with the filter on) a keep list that no row of this database can match stop the run; a token this database merely lacks, such as any k__ token on a GTDB database, gives a warning. See keeping and discarding taxa.

If the filter removes every taxon from a table that had some, rule finalize_otu_table stops with an error that names the tokens and lists the d__ values the table holds, and no table is written; see Troubleshooting. A table that is empty before the filter, because no sample has a taxon at tax_lev, is still written, with a warning in logs/finalize_otu_table.log.

Per-taxon read extraction

  extract_taxa: []

Key	What it does	Default	Notes
extract_taxa	NCBI scientific names whose reads should be pulled out into their own FASTQ files	[]	[] skips extraction entirely. Example: [Bacteria, Archaea].

Each name must be exactly the name of one row in the sample's Kraken2 report (02.classification/kraken2/{sample}_report.txt), and it is looked up at whatever rank it appears, so a broad group (Bacteria) and a single species (Escherichia coli) are named the same way. KrakenTools then extracts every read under that taxid, including its children, strains among them. A name found in several rows of a report stops the run, unless the rows are a genus and a subgenus of the same name (then the genus is extracted), and a single species gives only an approximate read set; the read extraction page explains both. Output files are named after a slug of the taxon (lower-case, spaces to underscores). A taxon absent from a sample yields empty files rather than a failure. See per-taxon read extraction.

---

Resources [shared]

resources:
  threads_default: 4
  mem_mb_default: 2000
  threads:
    kraken2: 16
    # … one entry per rule
  mem_mb:
    kraken2: 20000
    # … one entry per rule

Per-rule CPU and RAM live in one block with global fallbacks. Rules look their values up by name; a name that is absent falls back to threads_default or mem_mb_default. These values feed straight into Snakemake's resource accounting, so they are what a cluster executor turns into job requests.

Key	What it does	Default
threads_default	Threads for any rule with no threads entry	4
mem_mb_default	Memory in MB for any rule with no mem_mb entry	2000

Threads

Rule name	Mode	Default	Notes
multiqc	shared	2
falco	amplicon	2
bowtie2	amplicon	8	Covers both build_phix_index and the per-sample rm_phix mapping.
cutadapt	amplicon	4	Covers both trim_primers and the in-silico PCR probe.
pick_trunclen	amplicon	4
dada_quality_plots	amplicon	4
dada_filter	amplicon	4
dada_seqtab	amplicon	8
dada_length_filter	amplicon	2
target_extract	amplicon	8	Metaxa2 or ITSx.
assign_taxonomy	amplicon	8	With a fixed amplicon.seed the result does not depend on this value (for sintax from v2.5.0 on).
coi_screens	amplicon	8	COI only: VSEARCH against MIDORI2 (E) and the set of bacterial COI homologues (B).
aggregate_read_counts	shared	2	Present in the config template but not read — the rule declares no threads; see the note below.
phylo_input	amplicon	1	Reads two text tables, writes a FASTA.
phylo_align	amplicon	4	MAFFT.
phylo_tree	amplicon	1	Applies to the iqtree and raxml-ng backends only. One thread is a reproducibility setting, not a performance one: IQ-TREE's search depends on when each thread finishes, so two runs of the same alignment on several threads can end on different, about equally good trees, while two runs on one thread give the same tree. At amplicon scale one thread costs little: a minute or two for a few hundred ASVs. Raise it if you have thousands of ASVs and can accept a tree that is not exactly reproducible; the gain is limited anyway, because IQ-TREE parallelises across alignment columns and a 16S alignment is only a few hundred wide. This is the one key whose fallback is not threads_default — a rule that does not name it here gets 1, not 4. RAxML-NG clamps it further to what its own --parse step recommends (usually 1 for an amplicon alignment; it scales with the number of distinct alignment patterns) because it terminates when given far too many threads for a short alignment. The fasttree backend ignores this key — see the note below.
phylo_export	amplicon	1	Validates the tree, writes the Newick.
phylo_qc	amplicon	1	Reads the tree once for the QC report.
decontam_phix	shotgun	6	BBDuk scales poorly past 4–6 worker threads on a 5 kb reference; run more samples in parallel instead.
build_host_index	shotgun	16
decontam_host	shotgun	16
fastp	shotgun	16	Read by the trim_adapters rule. fastp scales poorly past about 16 threads.
kraken2	shotgun	16	Overridden by shotgun.kraken.threads when that is not null.
kmer_evidence	shotgun	1	Only when Kraken2 runs, whether the k-mer gate is on or off. One job per sample; the script reads that sample's Kraken2 report and uses one thread.
bracken	shotgun	4
finalize_otu_table	shotgun	1	One pass over the Bracken reports, building otu_table.tsv and otu_table.biom.
extract_taxon_reads	shotgun	4
compress_extracted	shotgun	6	pigz parallelism for the extracted read files.
download_sra	shotgun	6	fasterq-dump parallelism; past about 8 the bottleneck is disk, not CPU.
sylph_sketch	shotgun	1	Only when sylph runs. One thread is a reproducibility setting: with two or more threads sylph removes duplicate k-mers in parallel and the sketch differs slightly from run to run, while one thread gives the same sketch every time and still takes only seconds per sample. Snakemake sketches several samples at once. Like phylo_tree, it falls back to 1, not threads_default.
sylph_profile	shotgun	16	Only when sylph runs. One job for all samples and databases; the thread count changes no number.
sylph_taxprof	shotgun	1	Only when sylph runs. sylph-tax uses one thread.
sylph_tables	shotgun	1	Only when sylph runs.
sylph_read_tracking	shotgun	1	Only with profiler: sylph; under kraken2 and both, aggregate_read_counts writes read_tracking.txt.

Memory

Rule name	Mode	Default (MB)	Notes
multiqc	shared	2000
bowtie2	amplicon	4000
pick_trunclen	amplicon	4000
dada_quality_plots	amplicon	4000
dada_filter	amplicon	2000
dada_seqtab	amplicon	16000	The heaviest amplicon step: error learning, denoising, merging and chimera removal all happen here across every sample at once.
dada_length_filter	amplicon	4000
target_extract	amplicon	8000
assign_taxonomy	amplicon	16000	The rdp path is the demanding one; sintax stays around 2–3 GB whatever the ASV count, except against MIDORI2 for COI, where it holds 8.43 GB at the peak. That figure is set by the reference, not by the ASVs: it was the same, to within 1 MB, for 767, 804 and 888 ASVs on 8 threads. The step took 71 to 86 s for 683 to 909 ASVs in the validation runs.
coi_screens	amplicon	12000	COI only: VSEARCH holds MIDORI2 GB273 (2.1 billion bases) and its k-mer index, 8.51 to 8.57 GB at the peak for 767 to 888 ASVs on 8 threads; the step takes about a minute (55 to 87 s for 683 to 909 ASVs in the validation runs).
aggregate_read_counts	shared	2000	Present in the config template but not read — the rule declares no mem_mb; see the note below.
phylo_input	amplicon	2000	Two small text tables in, one FASTA out.
phylo_align	amplicon	4000
phylo_tree	amplicon	8000	Sized for the most demanding backend. FastTree was measured at ~181 MB on a 10,000-sequence, 400-column alignment, one thread; RAxML-NG estimates ~3.3 GB at 30,000 ASVs by its own formula and prints its estimate in the phylo_tree log; IQ-TREE is not separately measured but its likelihood vectors scale comparably. Raise this if you set support: true, which was not part of that sizing.
phylo_export	amplicon	2000
phylo_qc	amplicon	2000
decontam_phix	shotgun	4000
build_host_index	shotgun	30000	BBMap index build for the host reference; a masked human genome needs roughly 24–28 GB.
decontam_host	shotgun	30000	Mapping against that index.
fastp	shotgun	8000	Read by the trim_adapters rule.
kraken2	shotgun	20000	Not read — the real figure is computed at start-up from memory_mapping and the actual hash.k2d size.
kmer_evidence	shotgun	2000	Only when Kraken2 runs. Holds one sample's Kraken2 report and draws its evidence plot.
bracken	shotgun	4000
finalize_otu_table	shotgun	2000
extract_taxon_reads	shotgun	2000
sylph_sketch	shotgun	4000	Only when sylph runs. Making a sketch uses a few tens of MB of memory for a few million read pairs (the sketch file itself is a few MB); deep libraries need more, hence the headroom.
sylph_profile	shotgun	null	Only when sylph runs. null means worked out at start-up: 1.2 × the size of the database files plus 4096 MB, because sylph holds every database in memory (about 34,600 MB for GTDB r232 c200 with the fungal database). The start-up line prints the figure. A number here overrides it, for example on a cluster.
sylph_taxprof	shotgun	2000	Only when sylph runs.
sylph_tables	shotgun	2000	Only when sylph runs.
sylph_read_tracking	shotgun	2000	Only with profiler: sylph.

Names here are lookup labels, not always rule names

Most entries match a rule name exactly. A few are lookup labels rather than rule names: bowtie2 covers build_phix_index and rm_phix, cutadapt covers trim_primers and amplicon_probe, and fastp is the label used by the trim_adapters rule — there is no rule called fastp, and a trim_adapters key here is ignored. Rules that declare no threads or mem_mb at all — link_reads, count_reads_*, revcomp_primers, the fetch_* and convert_* reference rules, prepare_midori2_coi_sintax, check_taxonomy_reference, samples_manifest and aggregate_read_counts — take Snakemake's own defaults and ignore any entry of the same name. Adding a key for a rule that does not look it up has no effect. Two keys break the other way: phylo_tree and sylph_sketch fall back to 1, not threads_default, because there the thread count changes the result itself (the tree, the sketch) and not just how long it takes.

One more special case: threads.phylo_tree is read by the iqtree and raxml-ng backends but ignored by fasttree, which is pinned to a single thread in the rule itself. This is deliberate: FastTree's parallel build is documented as non-deterministic, both between runs and against the serial build, and it does not speed up the maximum-likelihood phase beyond about three cores. The pin also keeps Snakemake from reserving cores the job cannot use. mem_mb.phylo_tree applies to all three backends normally.