Add taxonomy summary table with per-ASV barrnap/decontam/filter annotations - #1063
Conversation
…ations Adds summary_tables/ampliseq.taxonomy.<classifier>.<database>.tsv.gz (+ Parquet), a slim, consistent-schema reformat of every classifier/database in ch_tax_tsv (Kraken2 excluded), joined with per-ASV annotations: barrnap's winning rRNA domain, decontam's contaminant call, individual filter pass/fail booleans, and ampliseq_accept (whether an ASV survived the whole standard filtering chain end to end). Individual filter columns are kept alongside ampliseq_accept so the pipeline's own filtering choices can be recombined differently if wanted. Fixes an existing gap where DUCKDB_TABLE2PARQUET only converted the counts table to Parquet, since it ran before any other summary table existed. Generated by Claude
sintax.nf.test/vsearch_lca.nf.test snapshots regenerated via real nf-test runs (verified real content: correct fungal ITS taxonomy, sh column preserved, barrnap_domain=euk, and this run exercised passed_length_filter_itsx for the first time). The other 11 affected profiles (every profile that runs a non-Kraken2 classifier, except default/multidb/pplace already covered in the previous commit) are hand-edited rather than run locally: for each, traced the exact classifier/database name(s) through the actual channel-building code (including the "user" database key used for every *_ref_tax_custom path, and PPLACE's two distinct database tags for --pplace_tree vs --pplace_sheet) and added the same two new snapshot entries (BUILD_ASV_ANNOTATIONS/ SUMMARY_TABLE_TAXONOMY versions, new summary_tables/ file paths) already confirmed correct by the real runs. Safe because nf-test's snapshot().match() is structural JSON equality, not byte-for-byte text comparison. Not locally re-verified -- CI will confirm or correct any misprediction. Also fixes a real YAML bug in summary_table_taxonomy/meta.yml: an unquoted "[classifier: ...]" in a description string was parsed by the YAML/prettier tooling as an attempted flow-sequence, caught by prek only after this change (meta.yml isn't read by Nextflow at runtime, so no test run had exercised it). Generated by Claude
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
…y CI The prior push's CI run (7 failed shards) surfaced a genuine bug, not just snapshot mispredictions: BUILD_ASV_ANNOTATIONS crashed whenever barrnap's summary.tsv had zero data rows -- a real, common case for ITS-only amplicon runs (barrnap targets rRNA, so zero hits across every ASV is expected, not a freak edge case), caused by a dplyr rowwise()/mutate()/c_across() corner case on empty input. Fixed by guarding nrow(evals) > 0 before that block. Also corrected 6 snapshots the earlier hand-edit pass missed, all verified by real local nf-test runs rather than hand-edited again: - doubleprimers/multiregion/novaseq: missing BUILD_ASV_ANNOTATIONS' versions entry -- it runs unconditionally regardless of classifiers, so these weren't actually "unaffected" as assumed. - fasta/savont/savont_independent: missing the bare "summary_tables" directory entry (never existed before in these -- PR A scoped SUMMARY_TABLE_COUNTS to asv_calling == "dada2" only, and fasta uses --input_fasta with no count table), plus a DUCKDB_TABLE2PARQUET versions entry that only running for real (not hand-editing again) caught. Generated by Claude
BUILD_ASV_ANNOTATIONS and SUMMARY_TABLE_TAXONOMY version entries were missing from default.nf.test and glosed.nf.test's snapshots -- stale from before this branch was rebased onto current dev. No functional change; verified all 21 test files pass locally. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
8445c6b to
8541192
Compare
ampliseq lists Added entries by ascending PR number, so nf-core#1063 belongs above nf-core#1065. The rebase onto dev resolved the conflict the other way round. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
d4straub
left a comment
There was a problem hiding this comment.
Didnt check the R script details, just in-person-code-browsing. Looks good so far I think. Just a few comments
# Conflicts: # CHANGELOG.md # tests/default.nf.test.snap # tests/multidb.nf.test.snap # tests/savont.nf.test.snap # tests/sintax.nf.test.snap # tests/vsearch_lca.nf.test.snap
The consolidated DADA2 table (--consolidate_taxonomies, added in nf-core#1062) reaches SUMMARY_TABLE_TAXONOMY through ch_tax_tsv like any other classifier/database entry, so it needs no wiring of its own. Its "database" column names the database that won each ASV, which reads as the table's own database once the file name already carries that, so it is renamed on the way into the summary table. tests/multidb_mostspecific.nf.test asserts the consolidated summary table exists and carries the column. The three multidb snapshots pick up the summary tables the test profiles now produce. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
- conf/modules.config: null-safe saveAs filter, as suggested in review. - docs/output.md: the Parquet note moves out of the collapsed "Output files" block, where it was easy to miss. - docs/output.md: state that the barrnap/decontam/filter annotations are joined onto the taxonomy tables only, and how to combine them with the counts table. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
Both R scripts carried multi-paragraph comments that traced the mechanism the code already shows, then justified the justification. Each is cut to the one fact a reader needs at that line: why booleans stack in their own table, why filter membership is diffed rather than read, why the whole chain is diffed again for ampliseq_accept, why Species and Species_exact are folded into one column. These sit in script: blocks, so they are copied verbatim into the .command.sh a user inspects at runtime. The meta.yml descriptions lose the channel-wiring rationale, which belongs in workflows/ampliseq.nf rather than in a module's interface documentation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
The entry named the schema, the Parquet copies and the parameter that turns them off, all of which docs/output.md covers. It now says what the files are and links there. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
|
❌ nf-test failed with latest Nextflow versionNote Tests with Nextflow's latest version failed but it will not cause a CI workflow failure.
See the full run for details. |
There was a problem hiding this comment.
Looks almost good to me. Two nitpicks:
- the fruitless execution of BUILD_ASV_ANNOTATIONS (details below)
- now that we have those summery_tables folder, maybe we should also store "overall_summary.tsv" there (and rename it to maybe something like read_counts_passed.tsv)?
| // MODULE: Per-ASV annotations (barrnap domain call, decontam contaminant call, per-filter | ||
| // pass/fail, and the whole-chain "ampliseq_accept") for the taxonomy summary tables below | ||
| // | ||
| BUILD_ASV_ANNOTATIONS ( |
There was a problem hiding this comment.
I saw in the .snap files that three tests execute BUILD_ASV_ANNOTATIONS but not SUMMARY_TABLE_TAXONOMY and produce no additional output. Additional module execution without benefit seems pointless. Maybe that could be avoided?
The three tests in question (with the settings that I suspect causing that):
- doubleprimers (skip_dada_taxonomy = true & kraken2_ref_taxonomy = "greengenes")
- multiregion (skip_dada_taxonomy = true & sidle_ref_taxonomy = "greengenes88")
- novaseq (skip_taxonomy = true)
There was a problem hiding this comment.
Good catch, and you identified the right three.
BUILD_ASV_ANNOTATIONS is called unconditionally, but the only thing that reads its output is SUMMARY_TABLE_TAXONOMY, which consumes ch_tax_tsv with Kraken2 filtered out.
In all three of those profiles that filtered channel is empty -- Kraken2 only in doubleprimers, SIDLE never enters ch_tax_tsv in multiregion, and novaseq skips taxonomy altogether -- so the annotations get built and then dropped.
I would rather not gate it on the parameters, since ch_tax_tsv is filled from six separate places under six different conditions and a parameter expression would have to restate all of them, which is the maintenance trap you pointed at in the other thread.
Gating on the channel itself keeps the two in step:
ch_summary_tax = ch_tax_tsv.filter { meta, _tsv -> meta.classifier != "KRAKEN2" }
ch_annot_barrnap = ch_summary_tax
.map { _meta, _tsv -> true }
.first()
.combine( ch_barrnapsummary.ifEmpty( [] ) )
.map { _gate, summary -> summary }The gated channel replaces the module's first input and the other twelve stay as they are, so the module itself does not change.
I checked the operator behaviour on its own before wiring it in: an empty taxonomy channel leaves the process unexecuted, and several taxonomy tables still run it exactly once.
I will push this and confirm the three snapshots lose the process without anything else moving.
There was a problem hiding this comment.
Done in ecd0115.
BUILD_ASV_ANNOTATIONS now takes its first input through a gate derived from the channel SUMMARY_TABLE_TAXONOMY consumes, so it runs only when there is a taxonomy table to join the annotations onto.
The module itself is unchanged, and SUMMARY_TABLE_TAXONOMY reuses the same filtered channel rather than filtering a second time.
All three profiles you named now lose the module's entry from the versions file and nothing else moves in their snapshots:
doubleprimerspassedmultiregionpassednovaseqpassed
fasta was run as a control, since it does write a summary table, and its snapshot is byte-identical.
SUMMARY_TABLE_TAXONOMY is the only consumer of BUILD_ASV_ANNOTATIONS, and it takes ch_tax_tsv with Kraken2 filtered out. Where that filtered channel is empty the combine yields nothing, so no summary table is written and the annotations table is built for nobody. Three profiles are in that position: doubleprimers classifies with Kraken2 alone, multiregion goes through SIDLE, which never enters ch_tax_tsv, and novaseq skips taxonomy altogether. The gate is the filtered channel itself rather than a parameter expression. ch_tax_tsv is filled from six places under six different conditions, and a parameter expression would have to restate all of them and then stay in step with them. Taking the first emission of the channel the module's output is joined onto cannot drift the same way. The three snapshots lose the module's entry from the versions file and nothing else. Profiles that do write summary tables are unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01PYuhWC4HE6BP63fnsVQZHN
PR checklist
nf-core pipelines lint).nextflow run . -profile test,docker --outdir <OUTDIR>).nextflow run . -profile debug,test,docker --outdir <OUTDIR>).docs/usage.mdis updated.docs/output.mdis updated.CHANGELOG.mdis updated.README.mdis updated (including new tool citations and authors/contributors).Description
Closes the last piece of issue #1033: a per-classifier taxonomy summary table, joined with per-ASV
annotations from barrnap, decontam, and the optional post-processing filters.
summary_tables/ampliseq.taxonomy.<classifier>.<database>.tsv.gz(+ Parquet)One file per classifier/database already carried in
ch_tax_tsv(Kraken2 excluded -- its rankvocabulary and output shape are too different to normalize). A slim, consistent-schema reformat of
that classifier's own native taxonomy table:
asv_id,kingdom..species,confidence.QIIME2's
k__/p__-prefixedTaxonstring and phylogenetic placement's unranked, semicolon-joinedstring are both parsed into the same seven rank columns (best-effort for the latter).
sequenceisdropped (recoverable via
asv_id) and, for DADA2, the per-rank*_confidencecolumns are droppedtoo (DADA2-specific detail, not comparable across classifiers) -- both remain available in full in
the native per-classifier files elsewhere in the output directory.
Joined onto every one of these files, whenever the corresponding step ran:
barrnap_domain: the rRNA domain barrnap called with the lowest e-value for this ASV, blank ifnone was significant.
decontam_contaminant/decontam_not_contaminant: decontam's own contaminant call.passed_ssu_filter,passed_length_filter_asv,passed_codon_filter,passed_length_filter_itsx: whether this ASV survived each individual optional filter.ampliseq_accept: whether this ASV survived the whole standard filtering chain (decontam throughthe ITSx-region length filter) end to end -- a convenience alongside the individual columns above,
not instead of them, so the pipeline's own filtering choices can be recombined differently.
New module
BUILD_ASV_ANNOTATIONSbuilds these by diffing each filter's ASV table immediatelybefore/after it runs (none of
FILTER_SSU/FILTER_LEN/FILTER_CODONSemit a per-ASV pass/failcolumn natively) -- purely additive plumbing, no change to the filter chain's own logic, order, or
gating. Tidyverse throughout (per-type long-format stacking +
pivot_wider, split by value typebefore stacking so Parquet gets real
BOOLEAN/UTF8columns rather than everything coerced throughone generically-typed column).
Also fixes a gap in the existing
ampliseq.counts.tsv.gz/.parquetParquet conversion (#1052):DUCKDB_TABLE2PARQUETused to run immediately afterSUMMARY_TABLE_COUNTS, before this PR'staxonomy tables existed yet, so it silently never converted them. Moved to run once after every
summary table has been mixed into the shared channel.
Two real bugs found by review + testing, both fixed
DADA2_ADDSPECIESalways renames its own species call toSpecies_exact, only re-adding a nativeSpeciescolumn when assignTaxonomy's own taxlevels included one -- databases that rely onaddSpecies alone for species-level calls (e.g. RDP) would have silently had no
speciescolumn atall. Fixed by coalescing
species/species_exact(prefer the native assignTaxonomy call, fallback to addSpecies), verified against real RDP output.
SUMMARY_TABLE_TAXONOMYneeded.combine()against the annotations table -- without it, Nextflowsilently dropped every
ch_tax_tsventry past the first, since a queue channel that only everemits once isn't auto-broadcast against a multi-item channel the way a true value channel is.
Testing
tests/default.nf.test,tests/multidb.nf.test,tests/pplace.nf.test,tests/sintax.nf.test,tests/vsearch_lca.nf.testall run for real (docker), with manual content spot-checks beyond theautomated assertions -- e.g. confirming QIIME2's
k__/p__string and PPLACE's unranked stringboth parse into the right rank columns, and that the RDP
speciesfix actually populates real data.fasta,glosed,multi,savont,savont_independent,single,pacbio_its,pplace_sheet,qiimecustom,reftaxcustom,failed,iontorrent) have their.nf.test.snapsnapshots updated by hand rather than a fulllocal run, based on tracing the exact classifier/database name each profile produces through the
channel-building code (including the
userdatabase key used by every*_ref_tax_custompath).Flagging this explicitly in case CI surfaces a misprediction there.
nf-core pipelines lintandprek run -aboth clean, at the same baseline asdev(thepre-existing
nextflow config -o json/SAVONT_ASV issue and a couple of unrelatedfiles_unchangedchecks aren't introduced by this PR).