AmplicONT

ONT wf-16S Amplicon Post-Processing Pipeline

Modular, config-driven downstream post-processing, statistical analysis, diversity quantification, and Kraken/Pavian/Krona export for Oxford Nanopore Technologies (ONT) full-length 16S rRNA and ITS amplicon sequencing generated by EPI2ME wf-16s using either Minimap2 or Kraken2/Bracken classifiers.

Current pipeline version: 0.5.0.

The kreport pathway provides the canonical Kraken report, an opt-in Krona export, and an opt-in original offline Kraken-report explorer. The .kreport file remains the official Pavian interoperability artifact; the built-in HTML is not the upstream Pavian Shiny application.

Architecture & Overview

This repository transforms primary outputs from ONT's Nextflow-based wf-16s pipeline into reproducible statistical figures, tables, and Pavian-compatible Kraken reports (.kreport). Supports both minimap2 alignment-based and kraken2/bracken k-mer/Bayesian estimation upstream contracts across 16S and ITS databases. Treat species-level calls and richness estimates as conditional on the upstream classifier, reference database, and abundance threshold—not as independent validation of organism presence.

analysis/
├── 00_run_pipeline.R        # Orchestrator & CLI entry point
├── 01_qc_diagnostics.R      # Classification donut, read length distributions, filter diagnostics
├── 02_alpha_diversity.R     # Registry-driven alpha metrics, rarefaction, optional phylogeny
├── 03_beta_diversity.R      # Bray-Curtis, Jaccard, PCoA, PERMANOVA (adonis2) & betadisper
├── 04_taxa_composition.R    # Standardized long-format taxa tables, bars, stacks, heatmap
├── 05_ordination.R          # Hellinger PCA & Bray-Curtis NMDS (cohort mode)
├── 06_shared_taxa.R         # Core/unique taxa by group prevalence, UpSet plot
├── 07_kreport_pavian.R      # Kraken report, Krona/Pavian exports, taxonomy provenance
├── 08_faprotax.R            # Opt-in FAPROTAX 1.2.12 taxon-based inference
├── 09_faprotax2.R           # Opt-in official FAPROTAX-mapper 2.0 taxonomy adapter
├── install_packages.R       # Dependency checker & installer (--install)
├── utils/merge_wf16s_cohort.R # Standalone multi-bundle cohort merger utility
├── utils/classifier_contract.R # Classifier contract validation & population accounting
├── utils/krona_builder.py   # Stdlib-only offline Krona-compatible HTML builder
├── utils/verify_krona_correspondence.py # Exact Krona HTML/TSV semantic verifier
├── utils/kraken_report_viewer.py        # Stdlib-only offline Kraken-report explorer
├── utils/verify_kraken_viewer_correspondence.py # Independent viewer verifier
├── utils/taxonomy_sankey_renderer.py    # Stdlib-only canonical Sankey renderer
├── utils/verify_taxonomy_sankey_correspondence.py # Independent Sankey verifier
├── utils/taxonomy_sankey_client.js       # Offline browser model and SVG renderer
└── vendor/krona-2.8.1/     # Byte-pinned, unmodified Krona browser assets

Repository-level supporting material is grouped by purpose:

wf16s-inputs/ contains local producer bundles; only the compact Ambar Ayunda regression contract is tracked.
patch/ contains implementation plans, review evidence, and release-audit notes, indexed by version.
archived/ contains historical scripts, analyses, and legacy outputs.

Analytical Modes & Gating

The pipeline supports two execution modes:

Single-Sample Mode (single):
- Focuses on within-sample characteristics: read-length distribution, classification yield, alpha diversity indices, analytical rarefaction curve, rarefaction resamples, taxonomic composition, and Kraken report generation.
- Multi-sample comparative modules (03_beta_diversity.R, 05_ordination.R, 06_shared_taxa.R) are gracefully skipped with structured records in the manifest.
Cohort Mode (cohort):
- Requires $\ge 2$ samples and a valid metadata mapping (SampleID and Group).
- Executes PCoA ordination, group-level taxonomic stacked bars, composition heatmaps, UpSet plots, and non-parametric group comparisons.
- Statistical Gates:
  - PERMANOVA (adonis2) requires at least 2 groups with $\ge 2$ samples per group and residual degrees of freedom. It is paired with betadisper homogeneity of multivariate dispersion testing. A warning is issued if any group has $<3$ samples.
  - Optional beta rarefaction stability is marked completed only when the configured minimum success fraction (default 0.50) of iterations produces valid Procrustes-aligned coordinates; otherwise the diagnostic is skipped with an explicit reason.
  - Non-parametric alpha group tests (Wilcoxon rank-sum or Kruskal-Wallis with Benjamini-Hochberg FDR) require at least 2 groups with $\ge 3$ biological samples each.

Input Data Contracts

Supported upstream contract

Version 0.5.0 is validated against the tracked wf-16s minimap2 and Kraken2/Bracken abundance, assignment, taxonomy-cache, and parameter schemas represented by the Ambar Ayunda regression fixture and synthetic Kraken2/ITS fixtures. Inputs from historical or future wf-16s releases must satisfy the contracts below; compatibility with arbitrary schemas is not implied.

1. Species Abundance Table (`abundance_table_species.tsv`)

Tab-delimited table where the first column is tax.
tax contains exactly eight semicolon-delimited ranks: superkingdom;kingdom;phylum;class;order;family;genus;species
Every rank must be non-empty after trimming and must not have leading or trailing whitespace; use an explicit label such as Unknown when applicable.
One or more sample columns containing non-negative, integer counts (e.g. 80556 or 80556.0).
Synthetic unclassified row: Unclassified;Unknown;Unknown;Unknown;Unknown;Unknown;Unknown;Unknown.
Optional aggregate column total must reconcile exactly to the row sum across sample columns.

2. Read Assignments (`*.minimap2.assignments.tsv`)

Optional per-sample mapping. If omitted, abundance-only modules proceed normally and read length diagnostics are cleanly skipped.
Five tab-delimited columns without header:
1. Status (C = aligned/classified, U = unaligned/unclassified)
2. Read ID (unique string)
3. TaxID (canonical decimal string; TaxID 0 denotes read aligned but failed downstream identity/coverage QC)
4. Length field (exactly a positive decimal length such as 1481, or a numeric prefix and positive length such as 0|1481)
5. Lineage path (pipe-delimited)
Effective Classification: Defined as taxid > 0. Raw status C alone overstates classification; reads with status C and TaxID 0 are filtered out by upstream reference coverage/identity thresholds.
Plain-text and .gz assignment inputs are read in bounded two-pass chunks. Schema validation, physical error line numbers, and the final typed QC data frame are preserved while avoiding a whole-file text buffer.

3. Producer contract and optional bamstats

input.params_json is required. Version 0.5.0 accepts minimap2 and kraken2, species-rank output from the bundled NCBI database sets without custom reference or taxonomy overrides.
input.wf16s_output_root is optional. When supplied, the pipeline discovers each sample's bamstats.readstats.tsv.gz by its internal sample_name and partitions C + TaxID 0 reads into identity-only, coverage-only, and both-threshold failures.
Missing optional bamstats produces explicit NA fields; a configured invalid or ambiguous bamstats contract fails rather than selecting an arbitrary file.

4. Optional phylogenetic alpha inputs

input.phylogenetic_tree may name a rooted Newick tree with unique tips and finite, non-negative branch lengths. Branch lengths are never imputed.
input.phylogenetic_tip_map is an optional two-column TSV with the exact header TaxonPath<TAB>TipLabel. Without it, positive classified TaxonPath values must exactly equal tree-tip labels.
Every positive classified taxon must map. Unclassified cannot be mapped; unused tree tips are excluded from metric traversal while the original root and stem paths are retained; the excluded-tip count, method, and input hashes are recorded in alpha_phylogeny_status.tsv and the run manifest.
Omitting the tree is an intentional module-level skip for Faith PD, PSR, PSE, and 02e; it emits alpha_phylogeny_skipped.tsv rather than a partial PNG.

Analytical Denominators & Metrics

Classified-Only Denominator:
- Relative abundance for community composition and alpha diversity is computed relative to classified reads ($\sum \text{counts}[\text{taxid} > 0]$).
- Unclassified reads are kept separate and analyzed in the QC diagnostics module.
Alpha Diversity Indices:
- Observed species richness $S$
- Chao1 estimated richness
- Shannon index $H$ and Effective Number of Species ($e^H$)
- Simpson's index $D$ ($1 - \sum p_i^2$) and Inverse Simpson ($1 / \sum p_i^2$)
- Pielou's evenness $J = H / \log(S)$
- Fisher's log-series alpha
- Berger-Parker dominance ($\max(N_i) / N$)
- ACE, Heip's evenness, Smith-Wilson Evar, and the McIntosh diversity index
- Explicitly ordered Rényi entropy (default $q=1$) and Hill numbers (default $q=0,1,2$); $N_0=S$, $N_1=e^H$, and $N_2$ is inverse Simpson
- Optional Faith phylogenetic diversity, phylogenetic species richness (PSR), and phylogenetic species evenness (PSE) when the tree contract is satisfied
Rarefaction Resamples:
- Analytical rarefaction curves are computed with vegan::rarefy().
- Subsampling iterations with vegan::rrarefy() represent rarefaction resamples without replacement to assess community sampling sensitivity. They are never presented as independent biological replicates.
- rarefaction_resamples.tsv retains its legacy columns first and appends registry metrics. rarefaction_resamples_long.tsv is the canonical plotting table with explicit metric keys, labels, q parameters, validity, and reasons.
- alpha_metric_definitions.tsv, alpha_metric_status.tsv, and alpha_plot_registry.tsv make formulas, edge rules, and facet membership inspectable. The fixed figures are 02b (six core metrics), 02c (four estimators/evenness metrics), 02d (four evenness/entropy metrics), and optional 02e (Hill q=0/1/2 plus Faith PD, PSR, and PSE).
Richness overview:
- 02_richness_overview.tsv reports positive taxa, singletons, taxa with at most 10 reads, and the reads represented by that low-count tail.
- These summaries remain conditional on upstream filtering, database choice, taxonomy, and the abundance threshold; low count alone does not prove artifact.

Unclassified-read handling by module

Module	Handling
QC	Included in abundance accounting; all assignment reads are shown when assignments exist
Alpha diversity	Excluded
Beta diversity	Excluded
Taxonomic composition	Excluded from rank tables; included in the classification-fraction output
Ordination	Excluded
Shared taxa	Excluded
Kraken/Pavian/Krona	Included as the first `U` line or top-level `Unclassified` contribution and in total-read arithmetic

Composition figures and sidecars

This release adds deterministic classified-read-relative composition figures at phylum, family, and genus ranks. The new vertical stacked bars are written as matching PNG/TSV pairs:

04_phylum_stacked, 04_family_stacked, and 04_genus_stacked for single samples or cohorts;
04_<rank>_group_mean_stacked for cohort group means; and
04_heatmap_<rank> for each configured heatmap rank.

The existing single-sample horizontal artifacts (04a through 04d) remain unchanged. New sidecars retain TaxonPath identity, contextual display labels, sample/group order, stack/row order, denominator validity, and explicit IsOther flags. The all-read Unclassified row is not added to these figures; it remains in classification_fraction.tsv and its corresponding plot.

The defaults are:

composition:
  stacked_bar_ranks: ["phylum", "family", "genus"]
  stacked_bar_max_taxa: 15
  stacked_bar_min_taxa: 5
  stacked_bar_min_mean_relative: 0.005
  heatmap_rank: null
  heatmap_ranks: ["phylum", "family", "genus"]
  heatmap_top_n_taxa: 10
  heatmap_include_other: true
  heatmap_transform: "log10_relative"

For compatibility, a non-null legacy scalar heatmap_rank is migrated to a one-element heatmap_ranks vector with one warning. Supplying both keys is an error; the pipeline never guesses precedence. Heatmap columns are never clustered. Cohort rows use complete-linkage clustering only when there are at least two display rows; single-sample rows are abundance-ordered and are not clustered. A rank with no positive classified abundance produces a structured *_skipped.tsv sidecar instead of a fabricated image.

Zero-classified samples are represented as undefined composition: stacked and heatmap sidecars keep finite numeric zero abundances with ValidDenominator = FALSE, and those samples are excluded from taxon-selection means, cohort means, heatmap pseudocounts, and clustering distances. Heatmap columns are plotted as neutral NA cells with an explicit denominator-status annotation. The residual Other row is never positive for an invalid sample. When a cohort group has no valid samples, its group-mean rows remain present with SamplesUsed = 0 and MeanRelativeAbundance = NA; the plot marks that group as undefined. Other group means are arithmetic means of valid per-sample classified- read proportions, not read-depth-pooled proportions.

A positive-total sample with zero classified reads can be processed when the requested modules are limited to qc, composition, and kreport. Modules whose statistics require a positive classified denominator (alpha, beta, ordination, shared, and faprotax) fail preflight with the affected sample and module names instead of receiving undefined input. For a single zero-classified sample, the historical 04a-04d composition figures are also withheld and 04_legacy_single_sample_skipped.tsv records why.

FAPROTAX Functional Inference (opt-in)

Request the module explicitly with --modules faprotax, or include faprotax in a comma-separated module list. The default module list is unchanged. The module checks the optional microeco dependency only when it is requested and requires microeco >= 2.3.0 with its embedded FAPROTAX version exactly 1.2.12. This is FAPROTAX 1.2.12 via microeco, not FAPROTAX2-db, PICRUSt2, or a gene/pathway abundance predictor.

Only positive-count Bacteria and Archaea rows enter the mapping. The canonical unclassified row and classified non-prokaryotic rows are excluded from mapping but retained in the explicit read accounting. Pipeline superkingdom is passed to microeco as Kingdom; the NCBI kingdom rank is intentionally omitted. Function read counts are the sums of classified taxon counts carrying each function, with no rarefaction, copy-number correction, or cross-function renormalization. Functions overlap, so their percentages are independent fractions and must not be stacked or expected to sum to 100%.

The module writes 08_FAPROTAX/ with the following artifacts:

faprotax_function_abundance.tsv: long SampleID, Function, FunctionReadCount, PctEligibleProkaryoticReads, PctClassifiedReads, and PctTotalReads table. Each percentage uses the denominator named in its column; a zero eligible denominator yields zero percent.
faprotax_mapping_coverage.tsv: per-sample TotalReads, UnclassifiedReads, ClassifiedReads, EligibleProkaryoticReads, ExcludedNonProkaryoticClassifiedReads, FunctionMappedReads, FunctionUnmappedEligibleReads, EligibleProkaryoticTaxa, and FunctionMappedTaxa. The equations are ClassifiedReads + UnclassifiedReads = TotalReads, EligibleProkaryoticReads + ExcludedNonProkaryoticClassifiedReads = ClassifiedReads, and FunctionMappedReads + FunctionUnmappedEligibleReads = EligibleProkaryoticReads. FunctionMappedReads counts each eligible read once when its taxon has at least one function; it is not the sum across overlapping functions. Every function count is bounded by 0 <= FunctionReadCount <= EligibleProkaryoticReads for its sample.
faprotax_taxon_function_assignments.tsv: positive mappings with stable FeatureID values and the original TaxonPath.
faprotax_top_functions.png: independent bars faceted by sample for the configured top functions; it is not a compositional stack.
faprotax_provenance.json: package/database versions, method, source, denominator definitions, input filtering, and the overlap warning.

If no positive-count Bacteria or Archaea rows exist, the module writes only faprotax_skipped.tsv and records a structured skipped module result. A missing/incompatible microeco installation or embedded database version is a failure, not a skip.

These are taxon-based functional assignments/inference. They do not establish gene presence, pathway completeness, expression, or activity. Cite the original database/software paper: Louca et al. (2016). The tracked-fixture comparison with the official collapse_table.py workflow is recorded in faprotax1.2.12-concordance-review.md; the two engines must not be treated as byte-for-byte equivalent.

FAPROTAX2 taxonomy inference (opt-in, unreleased)

--modules faprotax2 runs the official FAPROTAX-mapper 2.0 beta in taxonomy mode. The existing faprotax backend remains available, and both functional backends remain outside the default module selection. The mapper distribution includes its own reference tree and model data; FAPROTAX-db 2.0 alone is not that runtime. The upstream audit records the official download, license evidence and SHA-256 identities.

The wrapper consumes the existing pipeline abundance contract. It selects positive-count Bacteria/Archaea rows, translates the eight NCBI ranks to the mapper's seven ranks by omitting the intervening kingdom, and retains the original TaxonPath and count basis. Internally generated Taxon_###### IDs identify taxonomy rows; they are not representative sequences.

Native Windows setup and example

Use native Windows R 4.5.3, Python 3.12, and the existing synchronized renv library. The optional bootstrap checks that baseline and provisions mapper assets and optional binary packages separately. It preserves the baseline library. Run from the repository root in PowerShell:

$env:RENV_CONFIG_AUTOLOADER_ENABLED = 'FALSE'
$env:R_LIBS = (Resolve-Path 'renv/library/windows/R-4.5/x86_64-w64-mingw32').Path
$env:LC_ALL = 'English_United States.utf8'
Rscript --version
python --version
Rscript --vanilla tests/setup_faprotax2_windows.R
if ($LASTEXITCODE -ne 0) { throw 'Optional mapper setup failed' }

Rscript --vanilla analysis/00_run_pipeline.R `
  --config tests/fixtures/faprotax2/config.taxonomy.yml `
  --modules faprotax2 `
  --output-dir 'output/faprotax2-taxonomy-repro' --validate-only
if ($LASTEXITCODE -ne 0) { throw 'Validation failed' }

Rscript --vanilla analysis/00_run_pipeline.R `
  --config tests/fixtures/faprotax2/config.taxonomy.yml `
  --modules faprotax2 `
  --output-dir 'output/faprotax2-taxonomy-repro'
if ($LASTEXITCODE -ne 0) { throw 'Taxonomy mapper failed' }
Rscript --vanilla tests/verify_faprotax2_run.R 'output/faprotax2-taxonomy-repro'

The baseline path above is the standard existing Windows project library; resolve and inspect the actual installed library before using an alternative path. Setup rejects a missing or mismatched baseline and reuses only hash-matched cached downloads. It does not restore the baseline automatically. The setup library is an overlay: shared project/R dependencies remain available, and each run records their resolved versions, locations and DESCRIPTION hashes. Keep result directories separate from mapper, reference and package caches. Choose a fresh result directory, or explicitly request --overwrite for owned results you intend to replace.

For all nine modules, use --modules qc,alpha,beta,composition,ordination,shared,kreport,faprotax,faprotax2. Single-sample comparative analyses retain their structured skips. To compare both functional backends on the same input, request --modules faprotax,faprotax2 and run tests/compare_faprotax_backends.R on that completed output. Detailed comparison evidence must go under a separate ignored output/ directory.

Configuration, outputs and units

Add the following block to an existing pipeline config; paths resolve relative to that YAML file. The public fixture config supplies a runnable example after the optional setup.

faprotax2:
  mode: taxonomy
  mapper_path: path/to/FAPROTAX-mapper_v2.0/faprotax.R
  reference_dir: path/to/FAPROTAX-mapper_v2.0/FAPROTAX_ref_data
  mapper_library: path/to/optional-mapper-library
  nthreads: 1
  timeout_seconds: 600
  top_n_functions: 20

The adapter pins the audited script and all six reference files, checks optional dependencies before output mutation, and invokes Rscript with argument vectors. It uses the reference's per-function defaults, normalize=none and no probability weighting. One thread is the verified default. The pipeline seed is not forwarded to the external mapper; two local unseeded fixture repeats were byte-identical, which does not guarantee determinism across inputs, platforms or binary builds.

09_FAPROTAX2/ contains the long function-abundance table, coverage table, positive feature/function assignments, top-function PNG and provenance JSON, plus the input/identity crosswalk, original mapper matrices/report and process stdout/stderr. The raw states are 1 / 0 / NA for predicted presence / predicted absence / unknown; annotation probabilities are unitless. Estimated function abundance is t(states with NA replaced by 0) %*% eligible taxon counts. Unknown states remain visible in the raw tables. Functions overlap and their independent abundances can sum above the input total.

Coverage separates eligible/excluded input and taxa with positive function states. KnownAbsenceOnlyAbundance means rows whose model states are all zero; absence is still a prediction. Coverage is not a placement-accuracy measure: unrecognized taxonomies can receive broad fallback predictions. With no eligible prokaryotes the module writes one structured skip record. Missing software or altered pinned references fail preflight.

These phenotype predictions do not establish genes, pathway completeness, expression, measured activity or pathogenicity. See the native production validation and independent backend comparison for exact tested commands, results and limits. CI exercises offline adapter contracts and real child-process failure handling; the licensed official-mapper integration was verified locally rather than added as an automatic online CI download.

Sequence route remains blocked

The official mapper accepts representative 16S FASTA plus matching feature abundances. The supplied local ASV datasets had exact FASTA/table correspondence and were Illumina/DADA2-derived. The native one-feature trial reached alignment and APPLES placement, then failed the mapper's reference-tree consistency check; no function table was produced. faprotax2.mode: sequence is therefore explicitly rejected. It has no verified pipeline usage example or representative-feature input schema exposed by this implementation. The sequence trial report records the portability fixes, failed gate and local reproducibility artifacts.

Kraken Report & Pavian Interoperability

Kraken2 abundance tables contain Bracken species estimates (E) plus raw unclassified reads (U), so their total is T = E + U. Raw Kraken assignment counts use N = C + U; the classified population absent from retained species estimates is recorded separately as D = C - E. Alpha diversity and taxonomic composition use the species estimates. Builtin HTML exports state this basis.

To combine compatible single-sample results, run Rscript analysis/utils/merge_wf16s_cohort.R --input RUN1 --input RUN2 --output-dir COHORT. Each repeated --input accepts a literal path, including commas and semicolons. The merger validates report and assignment accounting before publication or --validate-only success. Its generated config references input.cohort_merge_manifest; the pipeline checks the merged abundance hash, original source hashes, and per-sample file mappings, and inventories those inputs throughout the run. Metadata and analysis settings remain editable. Merger --overwrite requires an unchanged merger bundle without unrelated files; choose a new destination once analysis outputs have been added.

For local review of uncommitted code, pipeline runs and tests/verify_release_run.R accept explicit --allow-dirty. Verification still checks exact source hashes and artifacts; clean release verification remains the default.

Generates standard 6-column .kreport files directly from the abundance tree:
1. Percentage of total reads rooted at this clade
2. Number of reads in clade (direct + descendants)
3. Number of reads assigned directly to taxon
4. Rank code: U, R, D, K, P, C, O, F, G, S (Standard Kraken rank codes; kingdom is K)
5. NCBI Taxonomy ID
6. Indented scientific name (2 spaces per depth level)
NCBI Taxonomy Resolution:
- Offline default (network_mode: cache_only): resolves TaxIDs using local taxonomy_cache.json without internet requests.
- Assignment-derived TaxIDs are written to a run-local resolved cache; an offline run does not mutate the configured source cache.
- Opt-in refresh (--refresh-taxonomy): performs a local preflight and then one NCBI Entrez E-utilities resolution pass with bounded retry/backoff, per-request rate limiting, and transactional cache replacement. Multiple processes or clients sharing an IP must coordinate their aggregate NCBI request rate externally.
- Local refresh validation (--validate-only --refresh-taxonomy) makes no network requests or filesystem changes. It reports taxonomy_resolution=pending_online when nodes still require resolution; unresolved_policy: error is deferred until an online resolution pass has occurred.
- Explicit full online validation (--validate-only --refresh-taxonomy --online-preflight) performs the same read-only resolution and failure-policy checks as execution. --online-preflight is rejected during normal execution to prevent an accidental duplicate full-network pass.
- --online-preflight requires kreport among the requested modules; a validation request for modules such as qc alone is rejected instead of silently skipping the taxonomy check.
- Fatal request failures and malformed/inconsistent NCBI responses fail under both unresolved policies. Valid not-found, ambiguous, or context-mismatched records remain unresolved under warn and fail under error.
- Current NCBI domain records are accepted only as an alias for legacy superkingdom at the canonical Bacteria, Archaea, and Eukaryota roots; scientific-name and ordered-ancestry validation remains strict, and Kraken rank codes are unchanged.
- Refresh execution writes flushed, credential-safe operational evidence under the output parent's .wf16s-diagnostics/<transaction-id>/ directory. Non-empty evidence survives ordinary staging cleanup; normal transaction cleanup removes only an empty directory owned by that transaction, never the shared diagnostics root. These logs are not scientific outputs, and an empty transaction directory may remain after hard operating-system termination outside the managed process tree.
- Unresolved and conflicting mappings are exported explicitly. Conflicting assignment mappings retain the modal TaxID (minimum numeric TaxID on a tie) for backward-compatible Kreport output, but are labeled Conflicted/assignment_conflict, emit a run warning, and must not be treated as authoritative without review. Named nodes with TaxID 0 should be resolved before treating the report as taxonomy-complete in Pavian.
- The tracked cache-only reference run intentionally reports 46 unresolved taxonomy nodes and 26 conflicting lineage-to-TaxID mappings. These diagnostics are surfaced rather than silently discarded.
Interactive Sankey Visualization:
- Generated .kreport files can be uploaded to Pavian for interactive Sankey and sunburst diagrams.

Krona export (opt-in)

Krona export is an extension of the kreport module and does not add a separate module name. Enable it with either the CLI switch or YAML:

Rscript analysis/00_run_pipeline.R \
  --config config.yml \
  --modules kreport \
  --krona \
  --overwrite

krona:
  enabled: true
  render_html: true
  html_renderer: "builtin"  # builtin, kronatools, or auto
  executable: "ktImportText"

For each selected sample, R always writes a portable 07_Kreport/krona/<sample>.krona.tsv in the tab-delimited lineage format accepted by KronaTools ktImportText. When render_html: true, the default html_renderer: builtin uses analysis/utils/krona_builder.py; no Perl or KronaTools installation is required. It writes a deterministic, standalone <sample>.krona.html containing the Krona-compatible XML, JavaScript, and image resources. html_renderer: kronatools is strict and fails if the configured executable is unavailable. html_renderer: auto uses an available external executable and otherwise falls back to the builtin renderer. Set render_html: false to retain TSV/provenance only.

The vendored Krona 2.8.1 browser files are unmodified and SHA-256 pinned in analysis/vendor/krona-2.8.1/SOURCE.json; the wrapper is described as Krona-compatible, not as official Krona software. The module writes 07_Kreport/krona/krona_provenance.json with renderer policy/version, vendor manifest hashes, HTML status, and exact read-accounting reconciliation.

Builtin Krona XML records both magnitude (the clade total) and magnitudeUnassigned (the direct assignment at that node), including internal nodes and the dataset root. Krona artifact paths in provenance use the path_basis = "run_dir" contract: they are POSIX-style paths relative to the completed run root and carry SHA-256 hashes, so they remain valid after private module staging is removed or the completed run is moved as a unit. The pipeline renderer identity comes from VERSION; the upstream Krona identity comes from the pinned vendor manifest.

Krona magnitudes use direct positive-count taxonomy rows plus one explicit top-level Unclassified contribution when applicable. Ancestor clade counts are not emitted separately, and Krona generates its own root. The output is a visualization of the classifier-conditioned abundance table, not independent taxonomic validation.

Optional Pavian integration and offline Sankey renderer

Set pavian.enabled: true or pass --pavian with kreport requested. The pipeline keeps <sample>.kreport as the six-column Kraken report intended for official Pavian upload/use. It also writes taxonomy_resolution.tsv and, when enabled, creates 07_Kreport/pavian/<sample>.pavian.json, an optional <sample>.pavian.html, and pavian_provenance.json.

The primary built-in visualization output is the deterministic Sankey artifact under 07_Kreport/pavian/sankey/<sample>.sankey.json and, when enabled, <sample>.sankey.html. The Sankey model is a classified clade-read flow with persistent Other and Assigned above residual lanes. Every selected rank column and transition conserves classified reads, including hidden taxa, intermediate residual carries, direct counts, and zero-classified runs. Use pavian.sankey.ranks and pavian.sankey.max_taxa_per_rank to choose the default view; the standalone HTML also supports transient rank/top-N views.

pavian:
  enabled: true
  render_html: true
  sankey:
    enabled: true
    render_html: true
    ranks: ["D", "K", "P", "C", "O", "F", "G", "S"]
    max_taxa_per_rank: 10

The HTML is an original dependency-free static explorer, not the upstream Pavian Shiny application. It embeds the normalized JSON as Base64, uses a restrictive offline CSP, performs no taxonomy lookup or network request, and preserves resolved, unresolved, and assignment-conflict status from the resolution sidecar. The .kreport remains the real Pavian-compatible upload artifact; the built-in explorer and Sankey are reproducible local renderers, not a bundled Pavian server or a claim of upstream Pavian implementation equivalence. Set pavian.render_html: false to retain deterministic JSON and provenance without HTML. Sankey HTML is emitted only when both pavian.render_html and pavian.sankey.render_html are true; the effective states are recorded in provenance and manifest revision 4.

Quickstart

config.yml is the tracked Ambar Ayunda regression demonstration. config.example.yml is a copy-and-edit template pre-populated for that single-sample demonstration. Give cohort projects their own abundance table, metadata, assignment mapping, and a fresh output directory.

1. Audit Dependencies

On Linux, the locked igraph binary also requires the GLPK runtime library (libglpk40 on the Ubuntu release runner).

# Check core installed packages (exits non-zero if missing)
Rscript analysis/install_packages.R

# Install missing core packages from CRAN
Rscript analysis/install_packages.R --install

# Check optional module dependencies (e.g. FAPROTAX / microeco)
Rscript analysis/install_packages.R --modules faprotax

# Install optional module dependencies from CRAN
Rscript analysis/install_packages.R --install --modules faprotax

# Audit all dependencies including optional modules
Rscript analysis/install_packages.R --all

2. Zero-Mutation Validation Run

Rscript analysis/00_run_pipeline.R --config config.yml --validate-only

3. Execute Full Pipeline

# Run pipeline with default config (writes to output/)
Rscript analysis/00_run_pipeline.R --config config.yml

# Override output directory and allow overwrite
Rscript analysis/00_run_pipeline.R --config config.yml --output-dir output --overwrite

Use a new output root for publication runs. A non-empty output root may only be reused with --overwrite when it contains a valid pipeline-owned manifest; the rerun is staged transactionally and replaces the complete owned product set, while unrelated files are preserved. Normal publication runs also require a synchronized renv.lock environment and a clean maintained source tree. The development-only --allow-unlocked and --allow-dirty flags record the exception in provenance.

Cohort sample identity

For a cohort, the same sample identifiers must appear in all applicable inputs:

abundance columns:  Sample_A, Sample_B
metadata SampleID:  Sample_A, Sample_B
assignment keys:    Sample_A, Sample_B

Set mode: cohort, provide metadata with unique SampleID and non-empty Group columns, and ensure those IDs match the selected abundance columns exactly. metadata.example.tsv demonstrates the metadata format only; it is not paired with the tracked one-sample abundance fixture.

Validation

Rscript analysis/install_packages.R --restore
Rscript -e "files <- list.files('analysis', pattern='[.]R$', recursive=TRUE, full.names=TRUE); invisible(lapply(files, parse))"
python -m compileall -q analysis tests
Rscript tests/testthat.R
python -m unittest -v tests/test_ncbi_taxonomy.py tests/test_check_committed_whitespace.py
Rscript analysis/00_run_pipeline.R --config config.yml --validate-only

The test suite creates a temporary synthetic cohort at runtime to verify software behavior only; no biologically validated cohort dataset is included. The committed renv.lock pins the R 4.5.3 package graph, including test tooling and microeco; use analysis/install_packages.R --restore before release or publication runs. The lock covers this post-processing environment, not the upstream wf-16s workflow, ONT basecaller, reference database archive, or optional KronaTools executable.

Run-manifest schema v2, revision 5

Version 0.5.0 uses manifest schema v2 revision 5. It retains strict artifact accounting, SHA-256 and byte-size digests for every owned product, explicit preserved unowned output tracking, and physical output census verification. Each run also records an ordered source_files path/SHA-256 inventory for the maintained producer scope and an aggregate digest derived from exact path<TAB>sha256 records. exports.krona and exports.pavian state the resolved effective export contract, including the Sankey ranks, top-N limit, and parent/child HTML gates; Pavian is still an extension of kreport, not a ninth pipeline module. Cardinality-dependent fields such as samples, requested modules, package versions, input collections, module outputs/warnings, artifacts, and Krona provenance samples are always JSON arrays, including zero and one item. Failed fail-fast runs retain later requested modules as not_run records with null execution timestamps. Every v2 manifest records whether the environment is locked, the SHA-256 of the repository renv.lock used for the run, the supplied configuration file's SHA-256 and byte size, and available package identity metadata from the lockfile and installed DESCRIPTION files.

Publication safety includes exclusive file locking (.wf16s_output.lock) preventing concurrent pipeline execution against the same output root, isolated per-module private staging directories, and hashed sibling transaction journals (.*.wf16s_journal.json) that guarantee deterministic rollback or forward-recovery across aborted writes. Online taxonomy refresh is computed into the run-local resolved cache first, then commits the shared source cache under an inter-process lock; a sibling taxonomy transaction journal and hash comparisons restore or forward-recover an interrupted commit without overwriting an unrelated cache revision. True zero-mutation --validate-only checks perform in-memory AST syntax validation without writing bytecode or touch files.

The cohort merger uses the same destination lock and its own schema-validated sibling journal. It fingerprints the complete prior bundle before replacement, verifies the renamed backup and staged bundle, and retains recovery material whenever either changes. Resource preflight uses the requested modules in both validate-only and execution modes. resource_budgets names the units it limits: sample-by-sample matrix cells, sample-pairs multiplied by taxa, samples multiplied by admissible permutation rows, sample-by-taxon resampling cells multiplied by iterations, alpha long-table rows, and original-tree covariance cells.

The tracked assignment fixture retains its typed-column digest under the streamed reader. On this Windows R 4.5.3 host, the manual 250,000-row benchmark completed in 13.33 seconds, returned a 30.53 MiB data frame, and observed a 343.27 MiB process peak working set. Run Rscript tests/benchmark_assignment_parser.R --rows=250000 on each target host for local performance evidence.

Command-Line Options

Option	Description
`-c, --config PATH`	Path to YAML config (default: `config.yml`)
`-o, --output-dir PATH`	Override output directory
`-m, --modules LIST`	Comma-separated modules: `qc,alpha,beta,composition,ordination,shared,kreport,faprotax,faprotax2`
`--validate-only`	Run the same requested-module preflight checks without writing to disk
`--keep-going`	Continue executing remaining modules if one fails
`--refresh-taxonomy`	Enable online NCBI queries to resolve missing TaxIDs
`--online-preflight`	With `--validate-only --refresh-taxonomy`, perform full read-only online resolution; rejected during execution
`--allow-unlocked`	Development-only escape hatch for a missing/mismatched lock (recorded as unlocked)
`--allow-dirty`	Development-only escape hatch for maintained source changes (recorded as dirty)
`--allow-large-workload`	Override a configured resource budget after reviewing validate-only preflight estimates
`--krona`	Enable Krona TSV export and builtin or policy-selected HTML rendering (requires `kreport`)
`--pavian`	Enable official-Pavian-compatible `.kreport` integration plus the builtin offline explorer and Sankey (requires `kreport`)
`--overwrite`	Replace the complete set of pipeline-owned products in a validated output root

Legacy Scripts

Original standalone scripts (analyze_16s_improved.R, analyze_16s.R, convert_to_kreport.R, fetch_ncbi_taxonomy.py) remain preserved under archived/ for retrospective comparison. New analyses should use analysis/00_run_pipeline.R.

Release Scope

Version 0.5.0 covers validated NCBI/minimap2/Kraken2/Bracken species parsing, exact separate read and estimate accounting, single-sample summaries, deterministic multi-rank classified-read and estimate composition figures and heatmaps, synthetically and real-bundle tested cohort modules, offline-by-default taxonomy resolution, Pavian-ready Kraken reports, the opt-in Krona export, original offline Kraken-report explorer, the conservation-checked offline taxonomy Sankey renderer, an opt-in FAPROTAX 1.2.12 functional inference module (via locked microeco 2.3.0), and the standalone merge_wf16s_cohort.R multi-bundle merger utility. It also includes manifest schema v2 revision 5, producer-source inventory, bounded assignment parsing, physical output census verification, transaction journals, and multi-process cache locking. The project does not include differential abundance, phylogenetic UniFrac, FAPROTAX2-db, PICRUSt2, the upstream Pavian Shiny application, or biological validation of cohort statistics or species-level calls. Version 0.5.0 follows 0.4.9.

Unreleased changes add the opt-in FAPROTAX-mapper 2.0 taxonomy backend and manifest revision 6. Sequence mode remains gated by the failed native proof; no version bump or release tag is included in this migration.

See CHANGELOG.md for release history and CITATION.cff for citation metadata. The code is available under the MIT License.

Navigation Menu

Repository navigation

About

Repository files navigation

AmplicONT

Architecture & Overview

Analytical Modes & Gating

Input Data Contracts

Supported upstream contract

1. Species Abundance Table (`abundance_table_species.tsv`)

2. Read Assignments (`*.minimap2.assignments.tsv`)

3. Producer contract and optional bamstats

4. Optional phylogenetic alpha inputs

Analytical Denominators & Metrics

Unclassified-read handling by module

Composition figures and sidecars

FAPROTAX Functional Inference (opt-in)

FAPROTAX2 taxonomy inference (opt-in, unreleased)

Native Windows setup and example

Configuration, outputs and units

Sequence route remains blocked

Kraken Report & Pavian Interoperability

Krona export (opt-in)

Optional Pavian integration and offline Sankey renderer

Quickstart

1. Audit Dependencies

2. Zero-Mutation Validation Run

3. Execute Full Pipeline

Cohort sample identity

Validation

Run-manifest schema v2, revision 5

Command-Line Options

Legacy Scripts

Release Scope

About

Releases

Packages

Contributors

Languages

Name		Name	Last commit message	Last commit date
.github/workflows		.github/workflows
analysis		analysis
archived		archived
patch		patch
renv		renv
tests		tests
wf16s-inputs/output_AAy		wf16s-inputs/output_AAy
.Rprofile		.Rprofile
.gitattributes		.gitattributes
.gitignore		.gitignore
.renvignore		.renvignore
CHANGELOG.md		CHANGELOG.md
CITATION.cff		CITATION.cff
LICENSE		LICENSE
README.md		README.md
VERSION		VERSION
config.example.yml		config.example.yml
config.yml		config.yml
config_helga.yml		config_helga.yml
metadata.example.tsv		metadata.example.tsv
renv.lock		renv.lock

Latest commit

History

Folders and files

Repository files navigation

AmplicONT

Architecture & Overview

Analytical Modes & Gating

Input Data Contracts

Supported upstream contract

1. Species Abundance Table (abundance_table_species.tsv)

2. Read Assignments (*.minimap2.assignments.tsv)

3. Producer contract and optional bamstats

4. Optional phylogenetic alpha inputs

Analytical Denominators & Metrics

Unclassified-read handling by module

Composition figures and sidecars

FAPROTAX Functional Inference (opt-in)

FAPROTAX2 taxonomy inference (opt-in, unreleased)

Native Windows setup and example

Configuration, outputs and units

Sequence route remains blocked

Kraken Report & Pavian Interoperability

Krona export (opt-in)

Optional Pavian integration and offline Sankey renderer

Quickstart

1. Audit Dependencies

2. Zero-Mutation Validation Run

3. Execute Full Pipeline

Cohort sample identity

Validation

Run-manifest schema v2, revision 5

Command-Line Options

Legacy Scripts

Release Scope

About

Resources

Stars

Watchers

Forks

Releases

Packages

Contributors

Languages

1. Species Abundance Table (`abundance_table_species.tsv`)

2. Read Assignments (`*.minimap2.assignments.tsv`)