Skip to content

DADA2 taxonomy consolidation across multiple --dada_ref_taxonomy databases - #1062

Open
erikrikarddaniel wants to merge 7 commits into
nf-core:devfrom
erikrikarddaniel:dada2-consolidate-taxonomy
Open

DADA2 taxonomy consolidation across multiple --dada_ref_taxonomy databases#1062
erikrikarddaniel wants to merge 7 commits into
nf-core:devfrom
erikrikarddaniel:dada2-consolidate-taxonomy

Conversation

@erikrikarddaniel

Copy link
Copy Markdown
Member

PR checklist

  • This comment contains a description of changes (with reason).
  • If you've fixed a bug or added code that should be tested, add tests!
  • If you've added a new tool - have you followed the pipeline conventions in the contribution docs
  • If necessary, also make a PR on the nf-core/ampliseq branch on the nf-core/test-datasets repository.
  • Make sure your code lints (nf-core pipelines lint).
  • Ensure the test suite passes (nextflow run . -profile test,docker --outdir <OUTDIR>).
  • Check for unexpected warnings in debug mode (nextflow run . -profile debug,test,docker --outdir <OUTDIR>).
  • Usage Documentation in docs/usage.md is updated.
  • Output Documentation in docs/output.md is updated.
  • CHANGELOG.md is updated.
  • README.md is updated (including new tool citations and authors/contributors).

Description

Replaces the "first-listed database wins" placeholder from #1056 with two real per-ASV consolidation
algorithms, used when --dada_ref_taxonomy lists more than one database:

  • --consolidate_taxonomies most-specific: per ASV, whichever listed database resolved the deepest
    taxonomic rank wins.
  • --consolidate_taxonomies score: per ASV, whichever listed database reported the highest
    assignTaxonomy bootstrap confidence wins.

Both break ties by declaration order in --dada_ref_taxonomy. Scoped to DADA2 only for now --
consolidating across different classification tools (DADA2 vs. SINTAX vs. QIIME2 etc.) is tracked
separately in #1061, and two other consolidation strategies (driver-table routing, barrnap-based) in
#1060.

Rank vocabulary harmonization

Different databases can use different rank names for DADA2's assignTaxonomy taxLevels -- e.g. PR2
uses Domain,Supergroup,Division,Subdivision,Class,... instead of the standard
Kingdom,Phylum,Class,.... most-specific treats Kingdom/Domain and Phylum/Division as
filling the same rank slot; Supergroup/Subdivision are deliberately not counted, so a database
with more intermediate rank names doesn't win purely by having extra columns. Verified against real
downloaded reference data (see docs/usage.md for the full rationale) -- the note there also flags
the known asymmetry for databases with both Kingdom and Domain populated, and the Kingdom→Domain
rename itself is tracked separately in #1059.

Worth calling out: on the test_multidb fixture with silva,gtdb,pr2 combined, PR2 won 71 of 347
ASVs outright (GTDB 196, SILVA 80) -- it competes meaningfully even on a bacterial-focused test
dataset, not just a token presence. Would be genuinely interesting to see how this plays out on real,
larger datasets.

Testing

  • tests/multidb.nf.test extended with --consolidate_taxonomies most-specific/score blocks
    (gtdb,rdp) and a new silva,gtdb,pr2 block using the oldest database versions that still exercise
    PR2's Domain/Supergroup/Division/Subdivision scheme.
  • Regression-checked against tests/default.nf.test and the existing single-mode test_multidb test.
  • nf-core pipelines lint/nextflow config -o json currently fails on this branch (and on dev) due
    to a pre-existing, unrelated Nextflow 26.04.x bug with SAVONT_ASV's closure-in-list publishDir --
    not something this PR introduces or can fix. nextflow lint . passes clean at both the declared
    minimum (26.04.0) and latest (26.04.6) Nextflow versions.

erikrikarddaniel and others added 3 commits August 20, 2026 18:42
Replaces the "first-listed database wins" placeholder (from nf-core#1056) with
two real per-ASV consolidation algorithms across a comma-separated
--dada_ref_taxonomy list: most-specific (deepest resolved rank wins)
and score (highest assignTaxonomy bootstrap confidence wins). Scoped to
DADA2 only for now; consolidating across different classification
methods (DADA2 vs. SINTAX etc.) is left for a future PR.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
- Sort per-database tables by --dada_ref_taxonomy order before binding,
  so the consolidated TSV's column order is deterministic instead of
  depending on which parallel DADA2 task happens to finish first.
- Document ASV_tax.consolidated.<method>.tsv in the ITSx output-files
  list too, not just the non-ITSx one.
- Extend the "--consolidate_taxonomies has no effect" warning to also
  fire when --dada_ref_tax_custom is set, since it collapses the
  database list to one regardless of --dada_ref_taxonomy's value.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
most-specific now treats Kingdom/Domain and Phylum/Division as filling
the same rank slot when counting how deep a database resolved an ASV
(e.g. PR2's Domain,Supergroup,Division,Subdivision,... vs the standard
Kingdom,Phylum,...). Supergroup/Subdivision are deliberately not
counted, so a database with more intermediate rank names doesn't win
purely by having extra columns. Verified against real PR2 v5.1.0/
SILVA v138.2 reference data. Documented in docs/usage.md, including
the known asymmetry for databases with both Kingdom and Domain
populated (tracked in issue nf-core#1059), and the still-open cross-tool
consolidation question (issue nf-core#1061).

Adds a silva,gtdb,pr2 test (oldest versions that still exercise PR2's
Domain/Supergroup/Division/Subdivision scheme) alongside the existing
gtdb,rdp test.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@erikrikarddaniel
erikrikarddaniel marked this pull request as ready for review August 21, 2026 05:22

@d4straub d4straub left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had a first look.
Will let AI do the same later today.
Will also test nf-core pipelines lint again, because last time I didnt see any trouble... probably I missed something.

Comment thread docs/usage.md Outdated
[`--dada_ref_taxonomy`](https://nf-co.re/ampliseq/parameters#dada_ref_taxonomy) accepts a comma-separated list of database keys (e.g. `gtdb,silva`) to run DADA2 classification against several databases at once. Only the first-listed database feeds downstream QIIME2 analysis (filtering, diversity, barplots, ANCOM) and R objects -- consolidating results across multiple databases into one is planned as a future addition. For the exact versioned key of each database (e.g. `gtdb=R11-RS232`), see [`conf/ref_databases.config`](https://github.com/nf-core/ampliseq/blob/master/conf/ref_databases.config).
[`--dada_ref_taxonomy`](https://nf-co.re/ampliseq/parameters#dada_ref_taxonomy) accepts a comma-separated list of database keys (e.g. `gtdb,silva`) to run DADA2 classification against several databases at once. By default, only the first-listed database feeds downstream QIIME2 analysis (filtering, diversity, barplots, ANCOM) and R objects; use [`--consolidate_taxonomies`](https://nf-co.re/ampliseq/parameters#consolidate_taxonomies) to instead pick a winning database per ASV (currently DADA2-only, does not compare across different classification methods). For the exact versioned key of each database (e.g. `gtdb=R11-RS232`), see [`conf/ref_databases.config`](https://github.com/nf-core/ampliseq/blob/master/conf/ref_databases.config).

**Rank vocabulary across `--dada_ref_taxonomy` databases** (relevant to `--consolidate_taxonomies most-specific`): most databases use the standard `Kingdom,Phylum,Class,Order,Family,Genus,Species` levels, but some don't -- e.g. PR2 uses `Domain,Supergroup,Division,Subdivision,Class,Order,Family,Genus,Species`. `Kingdom` and `Domain` are really the same rank under different names (a rename to `Domain` everywhere is tracked separately as [issue #1059](https://github.com/nf-core/ampliseq/issues/1059), since e.g. GTDB's "Kingdom" values are actually domains, not the classical Animalia/Plantae/etc.) -- `most-specific` treats them as interchangeable when counting how deep each database resolved an ASV. A database whose own scheme has both a `Kingdom` and a separate `Domain` column populated (e.g. older PR2 releases before v5.0.0) would get credit for both, a minor and accepted asymmetry rather than something specifically guarded against. Of PR2's three extra levels between Domain and Class, only `Division` is counted, treated as filling the same slot as `Phylum`; `Supergroup` and `Subdivision` are deliberately not counted, so a database with more intermediate rank names doesn't win purely by having more columns. This mapping isn't perfectly accurate everywhere in PR2's own tree -- e.g. for the Metazoa branch, PR2's own `Class`-labelled column actually holds classical phylum names (`Arthropoda`, `Mollusca`, etc.), one level below where `Division` sits -- but `Division` is where protist phylum-equivalent groupings (`Alveolata`, `Stramenopiles`, `Rhizaria`, ...) are actually found, and protists are PR2's primary use case. A database that resolves an ASV shallower than `Division`/`Phylum` in this scheme still gets partial credit from any deeper ranks it does share names with (`Class` onward), it just isn't counted as reaching the "phylum-equivalent" level.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thinks thats a bit too much a wall-of-text regarding

  • maintenance, split into several lines, will be rendered as one block, easier to see diffs
  • too much details imho, too verbose, can that info be shown differently, e.g. as a table, or bullet points or such?
  • some details seem not to be placed correctly here, i.e. I would not link an issue here at all (is there a precedence?) if possible to avoid it

Comment thread CHANGELOG.md Outdated
- [#1042](https://github.com/nf-core/ampliseq/pull/1042) - Added CI test coverage for `--addsh`, which previously had none (`test_pacbio_its` now also runs DADA2 taxonomy against the UNITE database it already uses for SINTAX) (by @erikrikarddaniel)
- [#1052](https://github.com/nf-core/ampliseq/pull/1052) - Added `summary_tables/ampliseq.counts.tsv.gz`, ASV counts in long format with consistent, lower-case column names, ready for analysis in R, Python or similar without pipeline-specific parsing; also written as Parquet by default, skip with `--skip_parquet_summary` (by @erikrikarddaniel)
- [#1056](https://github.com/nf-core/ampliseq/pull/1056) - `--dada_ref_taxonomy` now accepts a comma-separated list of databases (e.g. `gtdb,silva`), running DADA2 taxonomic classification against each; one full set of output files is published per listed database, and the first-listed database feeds every downstream step that expects a single taxonomy (consolidating multiple databases into one is planned as a future addition) (by @erikrikarddaniel)
- [#NN](https://github.com/nf-core/ampliseq/pull/NN) - Added `--consolidate_taxonomies` (`most-specific` or `score`) to pick a per-ASV winning database across multiple `--dada_ref_taxonomy` databases, instead of always using the first-listed one; the winning result is published as `dada2/ASV_tax.consolidated.<method>.tsv` (by @erikrikarddaniel)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

NN = 1062 :D

@d4straub

Copy link
Copy Markdown
Collaborator

I dont really get it. I run nf-core pipelines lint with nf-core/tools version 4.0.2 and nextflow version 26.04.0.12031 on the current dev and there is no error. How are you finding there errors?

@d4straub

Copy link
Copy Markdown
Collaborator

So the AI review yielded:

Issue A — most-specific winner-selection bias (real, partially acknowledged)

In consolidate_dada2_taxonomy.nf the score is rowSums(!is.na(rank_cols)), where rank_cols includes Species (populated by addSpecies, whose non-exact species calls DADA2 itself warns are unreliable) and counts both Kingdom+Domain and Phylum+Division (so a DB carrying both, e.g. pre-5.0.0 PR2, gets double credit). The author documented the Kingdom/Domain + Phylum/Division part as an "accepted asymmetry", but the Species part is undocumented and more consequential: a database can win purely by having emitted an often-untrustworthy species label. Separately, score compares raw confidence values that aren't normalized across DBs with different rank vocabularies/confidence scales.

-> me: This might be true, but its fine as is I think.

Issue B — consolidated output schema drift (real compatibility risk)

The module writes a union of all DBs' columns — so PR2-only columns (Domain, Supergroup, Division, Subdivision) appear — plus a trailing database provenance column. Downstream consumers of ch_dada2_tax (format_taxonomy.nf, QIIME2 import, phyloseq, etc.) were written for the fixed single-DB schema ASV_ID,<ranks>,confidence,<rank>_confidence. The added nf-tests only assert header presence and that both DBs appear in the database column; they do not exercise the consolidated file through the downstream modules, so any incompatibility there is currently untested.

-> me:
The nf-test seems to cover that, the CI test docker | 26.04.0 | 9/20 is Test [ec87a6de] '-profile test_multidb --consolidate_taxonomies most-specific' and that is importing the files. However, it might be indeed that the database column leaks through but it shouldnt?
While examining the tests, I found that docker | 26.04.0 | 8/20 is ending with

       > -[nf-core/ampliseq] Pipeline completed successfully-
      PASSED (853.403s)
    Snapshots:
      Obsolete snapshots can only be checked if all tests of a file are executed successful.

i.e. the test passes, but the Snapshot seems to have failed? Not sure how thats possible, but that might warrant a look into it? Seems wrong.

Addresses review feedback from @d4straub on nf-core#1062: the CHANGELOG entry's
#NN placeholder was never replaced with the real PR number, and the
rank-vocabulary paragraph in docs/usage.md was one unbroken wall-of-text
paragraph -- rewritten as bullet points (also easier to diff going forward),
dropped the nf-core#1059 issue link per his explicit request, and trimmed some of
the hedging detail while keeping the substantive rules.

Generated by Claude
@erikrikarddaniel

Copy link
Copy Markdown
Member Author

Reproduced with nf-core/tools 4.1.0 + Nextflow 26.04.6 (build 12646) -- exact error:

ERROR    Could not parse Nextflow config for '.../nf-core-ampliseq'. Possibly because of a Nextflow version mismatch?
         Command 'nextflow config -o json ...' failed with exit code 1
         ERROR ~ Cannot invoke method endsWith() on null object

NXF_VER=25.10.0 nf-core pipelines lint avoids it, which is what we used for this PR. You're on nf-core/tools 4.0.2 + Nextflow 26.04.0.12031 -- since that works cleanly, this might be specific to a later Nextflow 26.04.x patch (or the nf-core/tools 4.1.0 side of the nextflow config -o json call), not all of 26.04.x like we'd assumed. dev itself hits this too with our versions, so it's not something this PR introduces -- but could be worth its own upstream issue once we know exactly which version(s) are affected. Could you confirm your exact Nextflow patch?

@erikrikarddaniel

Copy link
Copy Markdown
Member Author

That message ("Obsolete snapshots can only be checked if all tests of a file are executed successful") is nf-test's own benign note whenever a run has any snapshot assertion mismatch within a test file -- not a sign of a real failure. Checked the actual job logs for the current (final, green) run's 8/20 and latest-everything/8/20 shards and neither contains "obsolete" or "snapshot" at all, so whatever you saw was very likely from an earlier, now-superseded CI run during development -- already resolved by the time it went fully green (20/20 both legs).

erikrikarddaniel added a commit to erikrikarddaniel/ampliseq that referenced this pull request Aug 21, 2026
Same lessons applied from review feedback on nf-core#1062: fill in the CHANGELOG
entry's real PR number (was #NNNN), and trim the summary-tables section of
docs/output.md, which had the same wall-of-text verbosity. Also fixes a real
doc gap caught while trimming: the taxonomy table's .parquet sibling was
never actually documented, only the counts table's was.

Generated by Claude
@d4straub

d4straub commented Aug 21, 2026

Copy link
Copy Markdown
Collaborator

nf-core/tools 4.1.0

This pipeline is written with and based on nf-core/tools 4.0.2, it cannot be expected to succeed with 4.1.0.
Using NXF_VER=25.10.0 nf-core pipelines lint with nf-core/tools 4.1.0 is just not right.

@d4straub

d4straub commented Aug 21, 2026

Copy link
Copy Markdown
Collaborator

That message ("Obsolete snapshots can only be checked if all tests of a file are executed successful") is nf-test's own benign note whenever a run has any snapshot assertion mismatch within a test file -- not a sign of a real failure. Checked the actual job logs for the current (final, green) run's 8/20 and latest-everything/8/20 shards and neither contains "obsolete" or "snapshot" at all, so whatever you saw was very likely from an earlier, now-superseded CI run during development -- already resolved by the time it went fully green (20/20 both legs).

That is wrong. Check out https://github.com/nf-core/ampliseq/actions/runs/32474576241/job/96748300564?pr=1062. The issue is still there and still valid and not resolved.

@erikrikarddaniel

erikrikarddaniel commented Aug 21, 2026

Copy link
Copy Markdown
Member Author

That message ("Obsolete snapshots can only be checked if all tests of a file are executed successful") is nf-test's own benign note whenever a run has any snapshot assertion mismatch within a test file -- not a sign of a real failure. Checked the actual job logs for the current (final, green) run's 8/20 and latest-everything/8/20 shards and neither contains "obsolete" or "snapshot" at all, so whatever you saw was very likely from an earlier, now-superseded CI run during development -- already resolved by the time it went fully green (20/20 both legs).

That is wrong. Check out https://github.com/nf-core/ampliseq/actions/runs/32474576241/job/96748300564?pr=1062. The issue is still there and still valid and not resolved.

It's reported as green now and I think it was this morning. In fact, everything is green so I don't understand the linting thing either.

When Claude checked this again, it found some old issue with that test which probably what you saw. But again, it's green now.

erikrikarddaniel and others added 3 commits August 25, 2026 14:01
…port attribution

Answers the open review question on nf-core#1062 (Issue B): confirmed the consolidated
output's "database" provenance column was leaking into bin/parse_dada2_taxonomy.r's
semicolon-joined taxonomy string as a bogus extra rank, corrupting the taxonomy fed
to every downstream QIIME2 step (filtering, diversity, barplots, ANCOM) whenever
--consolidate_taxonomies is set. While fixing it, found the same script was also
missing "Species_exact" (DADA2_ADDSPECIES's own column) -- a second, pre-existing
leak independent of consolidation, live on dev/master today. Replaced the
literal-name exclusion list (already missed a new column twice) with a
suffix-pattern match ("_confidence", "_exact") so the next added metadata column
doesn't repeat this.

Also fixes modules/local/summary_report.nf, which hardcoded the report's cited
reference-database title/file/citation to the first-listed --dada_ref_taxonomy
database even when consolidation picked a different database for most ASVs --
now states plainly when a consolidated result across multiple databases is used,
listing all contributing databases instead of naming one as if it were the sole
source.

Trims the --dada_ref_taxonomy/--consolidate_taxonomies docs/output.md paragraph
into bullet points, same verbosity feedback already applied to docs/usage.md.

Added a durable regression test: qiime taxa barplot emits one level-N.csv per rank
actually present in the taxonomy string it's given, so a level-8.csv is a general
signal that some non-rank column leaked in, not just a check for these two specific
columns. All changes verified end-to-end against real data (two full nf-test runs
of the multidb --consolidate_taxonomies most-specific profile, plus direct
inspection of the rendered summary_report.html).

The bin/add_sh_to_taxonomy.py counterpart of the Species_exact leak (--addsh,
unrelated to consolidation) is filed separately as nf-core#1064 rather than folded in
here.

Generated by Claude
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
bin/parse_dada2_taxonomy.r's Species_exact fix (previous commit) removes a
taxonomic-string leak that was live on every default-addSpecies run, not just
consolidated ones, so it also drops the "level-8" barplot directory from any
snapshot that captured it. Regenerated the three affected files against real
pipeline runs (default, savont, multidb's base non-consolidated test) --
tests/sintax.nf.test.snap needed no change (skip_dada_taxonomy=true there,
its own level-8 is UNITE-fungi's real 8-rank scheme, unrelated). Each diff
verified to contain only the level-8.csv/level-8.jsonp removal plus the
auto-generated timestamp field, nothing else.

Generated by Claude
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
…ef_taxonomy help text, drop stale phytoref

- Trim bin/parse_dada2_taxonomy.r's non-rank-column comment to the durable warning,
  dropping the historical "missed twice" narrative.
- Reorganize --dada_ref_taxonomy's help text by marker (SSU rRNA incl. PR2, ITS, COI,
  nifH), put the default (SBDI-GTDB) first, and enumerate every accepted lowercase
  database=version string so GUI users can see valid values without reading the config.
  Note PR2 also carries plastid/chloroplast 16S rRNA sequences.
- Remove phytoref from docs/usage.md's database table and its CITATIONS.md entry --
  no longer a supported DADA2 database (superseded by PR2's chloroplast 16S coverage),
  a stale leftover from nf-core#674.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants