MetaFlux: a unified short-read multi-marker amplicon and shotgun taxonomic profiling workflow

Abstract

MetaFlux v2.4.0 Two things in this release. An optional phylogeny stage for 16S data, which is off by default, and a set of changes that make an amplicon run reproducible from the raw reads to the final tree — that second part applies to every amplicon run whether or not the phylogeny stage is ever switched on. Reproducibility Until now, two runs of the same data with the same configuration could disagree in ways that had nothing to do with the data. Bowtie2 strips phiX on several threads and wrote the surviving pairs in whatever order the threads finished, so each run handed DADA2 the same reads in a different sequence. DADA2 numbers ASVs by decreasing abundance and keeps the input order for ties, so two ASVs with exactly the same total count could trade IDs between runs. Nothing biological changed, but the names did, and downstream that turned into a different tree. Three changes close it: --reorder on the phiX step. Bowtie2 now emits the surviving pairs in input order, which is what a single-threaded run produces. Over six alternating repeats the time cost was below measurement noise; memory rose by a few MB. The phylogeny input is sorted by sequence, so the alignment no longer depends on how the ASVs happen to be numbered. The tree search runs on one thread by default. A multithreaded maximum-likelihood search visits candidate trees in an order that depends on thread timing and can settle in a different optimum: two runs of the identical alignment at four threads ended 52 of 416 splits apart. resources.threads.phylo_tree is therefore the one resource key whose fallback is 1 rather than threads_default; raise it if you have thousands of ASVs and can accept a tree that is not exactly reproducible. Verified by deleting the 16S test output and running the pipeline twice from raw reads, two independent runs launched together: every FASTQ file at every stage, the ASV sequences, the per-sample counts, the alignment and the Newick tree came back byte-identical. One gap remains, and it is documented rather than fixed: with taxonomy.method: sintax on more than one thread, VSEARCH races its threads on a single random-number stream, so confidence values near the cutoff drift and a lineage can gain or lose trailing ranks — 16 of 211 ASVs between the two runs above, with every read count identical. Because the contaminant filter tests tokens at order and family rank, that can in principle change which ASVs reach the tree. Set resources.threads.assign_taxonomy: 1 if the ASV set itself has to be reproducible. The phylogeny stage This stage exists because @kwakikwak asked for it in discussion #2: phylogenetic beta diversity is a routine downstream step, and building the tree afterwards in a single-threaded R session is slow enough on a diverse community to be a real obstacle. The suggestion there was MAFFT plus FastTree. Both are in; the shipped default ended up being IQ-TREE instead, for reasons set out in that thread and on the Phylogeny page. Stage 85, 7.phylogeny/. After taxonomy, the retained 16S ASVs (the post-filter set in 6.taxonomy/, never the raw DADA2 output) are aligned with MAFFT and a maximum-likelihood tree is built. The tree is exported unrooted, as asv_16s.unrooted.nwk, with every ASV of the abundance table present exactly once as a tip — the workflow refuses to write the file otherwise. Diversity statistics are deliberately not computed: users prune, root and analyse in R, and the documentation carries a complete, tested R example (ape, phangorn, rbiom, vegan; no phyloseq) covering pruning after contaminant removal, midpoint rooting, rarefaction with distance matrices averaged over iterations, exact Faith's PD, UniFrac and Bray–Curtis ordination and their Procrustes comparison. Three backends, one config key. IQ-TREE 3 is the default (GTR+F+G4, fixed seed, --safe); FastTree and RAxML-NG are alternatives. Model selection (IQ-TREE ModelFinder, RAxML-NG MOOSE) is available on request and its choice is recorded, but never the default, so run time stays predictable. Alignment strategy follows the ASV count (L-INS-i below 200 sequences, FFT-NS-2 above) and can be forced. All thread counts are explicit; no tool is left to pick its own. Provenance and QC. phylogeny.params.json records the resolved aligner strategy, model, seed, thread counts, tool versions, input checksums and any pass-through arguments. stats/phylogeny/phylogeny_qc.json reports tip accounting, alignment statistics and a long-branch screen (pendant edge above Q3 + 3×IQR; root-to-tip above 1.5× the median on an in-memory midpoint-rooted copy). Flagged tips are listed, never removed. With fewer than four eligible ASVs the stage records a documented skip instead of failing. Validated on real data. Defaults were checked on the 16S test dataset (PRJNA305879, 211 ASVs): ModelFinder prefers a FreeRate model (TPM3u+R5), but the patristic distances that UniFrac and Faith's PD use are practically unchanged (r = 0.998), so the pinned default stands; absolute tree length, however, shifts 20–35 % with the rate model and must not be compared across models. Alignment masking was tested and rejected (28 % of tree length lost, no gain in monophyly). Backends differ more than settings within a backend: swapping two rows of the alignment moved IQ-TREE by 194 of 416 splits and RAxML-NG by 102, while FastTree moved by 10 and kept its branch lengths exactly — worth knowing when reading any single tree's topology. Per-sample diversity was far steadier than topology throughout (Faith's PD r ≥ 0.99 in every case). Scope. 16S only. 18S is deferred pending validation; ITS, gyrB and rpoB are excluded for marker-specific, cited reasons — the documentation explains each rather than promising them. Documentation New pages Phylogeny (16S) (with a 33-term glossary) and Phylogeny validation; updated configuration, output, troubleshooting and overview pages. Citations for MAFFT, IQ-TREE 3, ModelFinder, FastTree 2, RAxML-NG, DendroPy and the R packages used in the example. The Choosing a mode page now states plainly what shotgun mode is for and when a marker-based profiler is the better tool. Acknowledgements added: the MICROBE project (Horizon Europe, GA 101094353) and the Anthropic Open Source Program. New dependencies (all bioconda / conda-forge) mafft ≥ 7.505, iqtree ≥ 3.1, fasttree ≥ 2.1.11, raxml-ng ≥ 2.0, dendropy ≥ 5.0 — one conda environment per tool, created only when the phylogeny stage is enabled. Upgrading from v2.3.1 Existing configuration files parse unchanged, and a leftover root: key under phylogeny is rejected at parse time with an explanation, because the tree is exported unrooted only. Re-running an existing dataset will not reproduce a v2.3.1 run byte for byte, and this is intended. --reorder changes the order of the reads written by the phiX step, and therefore of everything derived from it, so ASVs that happened to be tied in abundance may be numbered differently than they were before. The sequences, the per-sample counts and the taxonomy are unaffected — compare an old run with a new one by sequence, not by ASV ID. No numbers change; their labels and their order on disk may.