16S-nf

A Nextflow (DSL2) pipeline for classifying Oxford Nanopore 16S rRNA amplicon data from microbial communities: vsearch identity clustering + spoa/racon/medaka consensus + reference/BLAST taxonomy assignment, plus Bray-Curtis/PCoA beta diversity. Sibling pipeline to edna-ont-nf -- same architecture, retargeted at 16S microbiome profiling instead of eDNA species detection, and emitting the same abundance_table.tsv schema so the iolani-bioinformatics frontend's abundance/rarefaction charts work against either pipeline's output unmodified.

Loosely follows wf-metagenomics' stages (ingest -> filter -> classify -> abundance + diversity report), but classifies via alignment against a targeted reference (wf-metagenomics' "minimap2 mode", implemented here via BLAST) rather than Kraken2 -- a better fit once noisy per-read ONT data has been consensus-polished, and it keeps the reference a plain FASTA a non-specialist can swap out, same as edna-ont-nf's --taxdb. wf-metagenomics computes alpha diversity/rarefaction per sample but not beta diversity -- the PCoA step here fills that gap.

Requirements

• Nextflow >= 23.10.0
• One of Docker, Singularity, or Conda (every process is containerized; pick a profile below)

Quickstart

nextflow run main.nf \
  --input samplesheet.csv \
  --taxdb silva_16s.fasta \
  --fwd_primer AGAGTTTGATCMTGGCTCAG --rev_primer TACGGYTACCTTGTTACGACTT \
  -profile docker

samplesheet.csv is a two-column CSV:

sample,fastq
sample_a,data/barcode17/*.fastq.gz
sample_b,data/barcode18/*.fastq.gz

fastq may be a single file or a glob matching multiple part-files (they're merged per sample before filtering). Glob paths resolve relative to the launch directory, not the samplesheet's location.

silva_16s.fasta is a plain FASTA of 16S reference sequences (e.g. exported from SILVA or NCBI's 16S targeted-loci database); a BLAST db is built from it at the start of every run.

Instead of supplying your own, --tax_db_source can auto-download a reference for you:

nextflow run main.nf \
  --input samplesheet.csv \
  --tax_db_source silva138 \
  -profile docker

--tax_db_source	Reference
custom (default)	Whatever FASTA --taxdb points to
silva138	SILVA 138.1 SSURef, restricted to Bacteria/Archaea, converted from RNA to DNA. Hits are called at genus level (see below)
refseq16s	NCBI RefSeq targeted-loci 16S export (Bacteria + Archaea)

silva138/refseq16s need network access from wherever the pipeline runs, and cache the downloaded fasta under --tax_db_cache (default assets/tax_db_cache/) so it's fetched once and reused on later runs. silva138 also fetches SILVA's taxonomy rank map (silva_138_taxmap.txt, alongside the fasta) to name each level of a hit's lineage. A cache from before that was fetched gets refreshed once. When --tax_db_source is anything other than custom, --taxdb is ignored (and not required).

Optionally, --metadata metadata.csv (sample,<any columns>, e.g. sample,group,site) gets joined onto the PCoA output for coloring points by group/site/etc. in the frontend -- it has no effect on classification.

Dev/debug entry point to run just the merge step in isolation, e.g. while iterating on a module:

nextflow run main.nf -entry MERGE_ONLY --input samplesheet.csv -profile docker

Parameters

Parameter	Default	Description
--input	(required)	Samplesheet CSV (sample,fastq)
--tax_db_source	custom	Reference database: custom (use --taxdb), silva138, or refseq16s (the latter two auto-download and cache)
--taxdb	(required if --tax_db_source custom)	16S reference sequences FASTA; a BLAST db is built from this each run
--tax_db_cache	assets/tax_db_cache	Where auto-downloaded reference databases are cached between runs
--outdir	results	Output directory
--fwd_primer / --rev_primer	null	Primer sequences (both 5'->3', as ordered) for cutadapt trimming; skipped if unset. Reads in either orientation are trimmed to what lies between the primers, and reads without both are dropped
--min_len / --max_len / --min_qual	1200 / 1800 / 10	chopper length/quality filtering thresholds -- sized for full-length 16S; narrow for a single V-region amplicon
--enable_read_stats	true	Read length and Q-score summary before vs. after filtering (final_report/read_qc_summary.html, read_stats.tsv, read_length_qscore_hist.tsv); purely informational
--cluster_id	0.86	vsearch identity threshold for read clustering -- tune per amplicon/primer set
--min_cluster	5	Minimum reads in a read-level cluster to attempt consensus -- a compute/quality floor, kept low so fragments can still be merged
--spoa_max_reads	500	Most reads per cluster used for the spoa draft, taken evenly through the cluster (0 = all). RACON still polishes with every read, and cluster_size counts every read
--merge_id	0.97	Identity at which per-cluster consensus sequences are merged (second vsearch pass, MERGE_CONSENSUS); 0 skips the merge
--enable_medaka	false	Use medaka-polished consensus instead of the racon consensus downstream
--container_registry / --dockerhub_registry	quay.io / docker.io	Registries the tool images are pulled from. The google_batch profile points both at the project's Artifact Registry mirror (see Container mirror); pass the defaults to pull from upstream
--blast_chunk_size	300	Consensus sequences per BLAST_TAX task, so a sample with thousands of clusters is searched on several VMs side by side (0 = one task per sample). The hits are the same either way
--min_pident	90	BLAST hits below this %identity are flagged low_identity, not dropped; also the cutoff for what counts as an identified taxon in the rarefaction curve and PCoA
--min_abundance	20	Merged clusters with fewer reads than this get low_abundance = true in the abundance table, not dropped (0 = off)
--min_rel_abundance	0	Merged clusters holding a smaller fraction of their sample's clustered reads than this get low_abundance = true, e.g. 0.001 = under 0.1% (0 = off)
--metadata	null	Optional CSV (sample,<columns>) joined onto the PCoA output for coloring

All defaults live in nextflow.config, not main.nf (see the comments there if you're adding a new one).

Pipeline

flowchart TD
    tax_db_source{"--tax_db_source"}
    tax_db_source -->|"custom (default)"| taxdb[/"--taxdb silva_16s.fasta"/]
    tax_db_source -->|"silva138"| FETCH_TAXDB_SILVA
    tax_db_source -->|"refseq16s"| FETCH_TAXDB_REFSEQ
    taxdb --> MAKEBLASTDB
    FETCH_TAXDB_SILVA --> MAKEBLASTDB
    FETCH_TAXDB_REFSEQ --> MAKEBLASTDB
    MAKEBLASTDB --> blastdb[("BLAST db")]

    reads[/"--input samplesheet.csv"/] -->|"sample,fastq rows"| fastqs[/"fastq(.gz) files\n(per-sample, referenced by each row)"/]
    fastqs --> MERGE_FASTQ --> CHOPPER --> CUTADAPT --> VSEARCH_CLUSTER
    MERGE_FASTQ --> READ_STATS
    CHOPPER --> READ_STATS
    READ_STATS --> READ_STATS_REPORT
    VSEARCH_CLUSTER -->|"per sample"| SPOA_CONSENSUS --> MINIMAP2_ALIGN --> RACON

    RACON --> medaka_check{"--enable_medaka?"}
    medaka_check -->|"false (default)"| consensus["consensus fasta"]
    medaka_check -->|"true"| MEDAKA --> consensus

    consensus --> merge_check{"--merge_id > 0?"}
    merge_check -->|"yes (default 0.97)"| MERGE_CONSENSUS --> merged["merged consensus fasta"]
    merge_check -->|"0"| merged

    merged --> BLAST_TAX
    blastdb --> BLAST_TAX

    BLAST_TAX --> BUILD_REPORT

    metadata[/"--metadata (optional)"/] --> PCOA
    BUILD_REPORT --> PCOA

    BUILD_REPORT --> report[["final_report/"]]
    PCOA --> report
    BUILD_REPORT --> confident[["confident/"]]
    BUILD_REPORT --> lowconf[["low_confidence/"]]
    BUILD_REPORT --> nohit[["no_hit/"]]

Loading

1. MAKEBLASTDB -- build a BLAST db (once per run) from --taxdb, or from FETCH_TAXDB_SILVA/FETCH_TAXDB_REFSEQ's output when --tax_db_source is silva138/refseq16s
2. MERGE_FASTQ -- merge multi-part fastq(.gz) files per sample
3. CHOPPER -- length/quality filter
   • READ_STATS (side branch, --enable_read_stats, on by default) -- read length and mean Q-score before vs. after filtering
4. CUTADAPT -- primer trimming (skipped if no primers supplied): one linked adapter with --revcomp, so reads from either strand are trimmed and come out forward-oriented
5. VSEARCH_CLUSTER -- de novo identity clustering
6. SPOA_CONSENSUS -- draft consensus per cluster from at most --spoa_max_reads of its reads (one task per sample, clusters in parallel across its cpus)
7. MINIMAP2_ALIGN + RACON -- alignment-based polish
8. MEDAKA -- ONT-specific polish, opt-in via --enable_medaka
   • MERGE_CONSENSUS (--merge_id, default 0.97, 0 skips) -- second vsearch pass that folds near-identical cluster consensus sequences back together, summing their read counts
9. BLAST_TAX -- taxonomy assignment against the run's BLAST db (one blastn per chunk of --blast_chunk_size consensus sequences; BUILD_REPORT joins each sample's chunks back into 07_taxonomy/{sample}.hits.tsv)
10. BUILD_REPORT -- per-run abundance table + QC summary html, and every consensus fasta gathered by BLAST-hit confidence
11. PCOA -- Bray-Curtis distance + classical PCoA across samples, optionally joined with --metadata

Consensus fasta headers are stamped as >{sample}_{cluster_id} sample={sample} cluster_size={n_reads} by whichever of RACON/MEDAKA produces the final consensus for a cluster.

Output layout

results/
  {sample}/
    00_merged/                merged fastq
    01_filtered/               chopper output
    02_trimmed/                 cutadapt output
    03_clusters/                vsearch read clusters (+ clusters.uc)
    04_draft/                   spoa draft consensus, one fasta per cluster
    05_racon/                   racon consensus, one fasta per cluster
    06_merged/                  merged consensus fastas + merge_map.tsv (unless --merge_id 0)
    06_consensus/               medaka consensus (only if --enable_medaka)
    07_taxonomy/                BLAST hits ({sample}.hits.tsv, all of the sample's clusters)
  blastdb/                     BLAST db built from --taxdb
  consensus_by_confidence/
    confident/                 best BLAST hit >= --min_pident
    low_confidence/            best BLAST hit < --min_pident
    no_hit/                    no BLAST hit at all
  final_report/
    abundance_table.tsv        one row per cluster: sample, cluster_size, best hit, tied_taxa (species tied for the best bitscore, if more than one), flag_reason, low_abundance
    run_qc_summary.html        cluster counts, per-sample flagged and low_abundance counts
    read_qc_summary.html       read length / Q-score, before vs. after filtering (--enable_read_stats)
    read_stats.tsv             the same numbers per sample and stage
    read_length_qscore_hist.tsv  length and Q-score histograms behind the platform's Read QC charts
    pcoa_coordinates.tsv       one row per sample: PC1-PC3 (+ % variance explained), plus any --metadata columns
  pipeline_info/                Nextflow timeline/report/trace

abundance_table.tsv is byte-for-byte the same schema edna-ont-nf produces (seq_id, sample, cluster_size, subject_id, pident, length, evalue, bitscore, stitle, tied_taxa, flag_reason, low_abundance) -- this is what lets the platform's existing abundance and rarefaction charts work against 16S-nf's output with no frontend changes. pcoa_coordinates.tsv is new to this pipeline (edna-ont-nf doesn't compute beta diversity).

SILVA hits. A SILVA hit's stitle is a lineage (Bacteria;Proteobacteria;...;Pseudomonas;<organism>), not a binomial, and SILVA's curated taxonomy stops at genus: the organism name is whatever the source record said, uncultured bacterium for most records. So when the reference is SILVA-style (any stitle containing ;), the table gains superkingdom, phylum, class, order, family, genus columns between tied_taxa and flag_reason, where edna-ont-nf puts its --taxdump lineage. Each hit is labelled Genus sp., or, when SILVA has no genus name for it, by its lowest named rank (Anaerolineaceae (family)). Placeholder names (uncultured, Incertae Sedis, Unknown ...) are left blank rather than treated as taxa. tied_taxa and the PCoA use the same labels, and so does the frontend.

Repo structure

main.nf                     entry point, samplesheet parsing, --help
workflows/16s.nf             subworkflow chaining all steps
modules/*.nf                one process per tool, one container each
bin/build_report.py         abundance table + QC html
bin/pcoa.py                 Bray-Curtis + classical PCoA, optional metadata join
bin/read_stats.py           read length/Q-score summary (stdlib only; also runs stand-alone)
scripts/mirror_images.sh    copies the pipeline's container images into our Artifact Registry (see Container mirror)
nextflow.config              param defaults, profiles, resource labels
nextflow_schema.json         JSON Schema describing every --param (for UIs/validation tooling)
conf/test.config             -profile test overrides (small synthetic dataset)
assets/                      your own local samplesheet/taxdb/metadata go here (gitignored, except NO_METADATA)
assets/tax_db_cache/         auto-downloaded SILVA/RefSeq reference fasta, cached across runs (gitignored)
tests/data/                  small synthetic dataset used by -profile test
.github/workflows/ci.yml     runs -profile test on push/PR

Container mirror

On Google Batch every task starts on a fresh VM and pulls its image, so one timed-out pull from quay.io or Docker Hub can end a long run. The google_batch profile therefore pulls from our own Artifact Registry repo (us-central1-docker.pkg.dev/iscc-400300/containers), which keeps each image's upstream path (biocontainers/blast:..., ontresearch/medaka:..., library/ubuntu:...) so only the registry host differs. The repo is shared with edna-ont-nf; images both pipelines use (same tags) are already there.

scripts/mirror_images.sh fills it. It reads the image list from the pipeline's own container directives, so it can't drift:

gcloud auth configure-docker us-central1-docker.pkg.dev   # once, for pushing
scripts/mirror_images.sh --check     # read-only: is every image there, same digest as upstream?
scripts/mirror_images.sh             # copy whatever is missing

Run it again whenever an image tag is added or changed, and before releasing (--check exits 1 if anything is missing). Every container must go through ${params.container_registry} or ${params.dockerhub_registry}; the script refuses to run if one doesn't, since it would silently bypass the mirror.

Testing

tests/data/ holds a small synthetic dataset (two samples, a handful of reads each -- the same fixtures edna-ont-nf uses for its own test profile, since this only needs to exercise pipeline mechanics, not real taxonomy) plus a metadata.csv so -profile test exercises the PCoA metadata join too:

nextflow run main.nf -profile test,docker

This populates results/ exactly like a normal run (see Output layout above), just against the small checked-in fixtures instead of your own data.

Note: the medaka container is ONT's own multi-arch image (ontresearch/medaka), not biocontainers. The biocontainers build is amd64-only and SIGILLs under Docker's emulation on Apple Silicon.