Adding ONT read support for ampliseq

.github/workflows/nf-test.yml

@@ -50,7 +50,7 @@ jobs:
         env:
           NFT_VER: ${{ env.NFT_VER }}
         with:
-          max_shards: 16
+          max_shards: 17
 
       - name: debug
         run: |

CHANGELOG.md

@@ -8,6 +8,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### `Added`
 
 - [#1021](https://github.com/nf-core/ampliseq/pull/1021),[#1023](https://github.com/nf-core/ampliseq/pull/1023),[#1025](https://github.com/nf-core/ampliseq/pull/1025) - Comparison of observed ASVs against expected sequences or their abundance is now available with `--expected_*` parameters (by @d4straub, reviewed by @erikrikarddaniel)
+- [#1026](https://github.com/nf-core/ampliseq/pull/1026) - Add support for Oxford Nanopore Technology (ONT) R10.4 sequencing (preferably with SUP basecalling) with Savont (by @d4straub)
 
 ### `Changed`
 

CITATIONS.md

@@ -14,23 +14,42 @@
 
 ## Pipeline tools
 
-### Core tools
-
-- [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
-
-> Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data [Online].
+### Preprocessing & QC
 
 - [Cutadapt](https://journal.embnet.org/index.php/embnetjournal/article/view/200/479)
 
   > Marcel, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet. journal 17.1 (2011): pp-10. doi: 10.14806/ej.17.1.200.
 
-- [Barrnap](https://github.com/tseemann/barrnap)
+#### Short read tools
 
-  > Seemann T. barrnap 0.9 : rapid ribosomal RNA prediction.
+- [FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
+
+  > Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data [Online].
 
 - [DADA2](https://pubmed.ncbi.nlm.nih.gov/27214047/)
+
   > Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJ, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016 Jul;13(7):581-3. doi: 10.1038/nmeth.3869. Epub 2016 May 23. PMID: 27214047; PMCID: PMC4927377.
 
+#### Long read tools
+
+- [Porechop_ABI](https://pubmed.ncbi.nlm.nih.gov/36698762/)
+
+  > Bonenfant Q, Noé L, Touzet H. Porechop_ABI: discovering unknown adapters in Oxford Nanopore Technology sequencing reads for downstream trimming. Bioinform Adv. 2022 Nov 21;3(1):vbac085. doi: 10.1093/bioadv/vbac085. PMID: 36698762; PMCID: PMC9869717.
+
+- [Chopper](https://pubmed.ncbi.nlm.nih.gov/37171891/)
+
+  > De Coster W, Rademakers R. NanoPack2: population-scale evaluation of long-read sequencing data. Bioinformatics. 2023 May 4;39(5):btad311. doi: 10.1093/bioinformatics/btad311. PMID: 37171891; PMCID: PMC10196664.
+
+- [Savont](https://doi.org/10.64898/2026.05.26.727271)
+
+  > Jim Shaw, Marie Riisgaard-Jensen, Kasper Skytte Andersen, Rasmus Kirkegaard, Morten Kam Dahl Dueholm, Heng Li. Sensitive long-read amplicon sequence variant recovery with savont. bioRxiv 2026.05.26.727271; doi: https://doi.org/10.64898/2026.05.26.727271
+
+### Default broad classification
+
+- [Barrnap](https://github.com/tseemann/barrnap)
+
+  > Seemann T. barrnap 0.9 : rapid ribosomal RNA prediction.
+
 ### Taxonomic classification and databases
 
 - Classification by [QIIME2 classifier](https://pubmed.ncbi.nlm.nih.gov/29773078/)

README.md

@@ -22,7 +22,7 @@
 
 ## Introduction
 
-**nfcore/ampliseq** is a bioinformatics analysis pipeline used for amplicon sequencing, supporting denoising of any amplicon and supports a variety of taxonomic databases for taxonomic assignment including 16S, ITS, CO1 and 18S. Phylogenetic placement is also possible. Multiple region analysis such as 5R is implemented. Supported is paired-end Illumina or single-end Illumina, PacBio and IonTorrent data. Default is the analysis of 16S rRNA gene amplicons sequenced paired-end with Illumina.
+**nfcore/ampliseq** is a bioinformatics analysis pipeline used for amplicon sequencing, supporting denoising of any amplicon and supports a variety of taxonomic databases for taxonomic assignment including 16S, ITS, CO1 and 18S. Phylogenetic placement is also possible. Multiple region analysis such as 5R is implemented. Supported is paired-end Illumina or single-end Illumina, PacBio, IonTorrent, and Oxford Nanopore data. Default is the analysis of 16S rRNA gene amplicons sequenced paired-end with Illumina.
 
 A video about relevance, usage and output of the pipeline (version 2.1.0; 26th Oct. 2021) can also be found in [YouTube](https://youtu.be/a0VOEeAvETs) and [billibilli](https://www.bilibili.com/video/BV1B44y1e7MM), the slides are deposited at [figshare](https://doi.org/10.6084/m9.figshare.16871008.v1).
 
@@ -34,7 +34,7 @@ On release, automated continuous integration tests run the pipeline on a full-si
 
 ## Pipeline summary
 
-By default, the pipeline currently performs the following:
+By default, for paired-end Illumina data, the pipeline currently performs the following:
 
 - Sequencing quality control ([FastQC](https://www.bioinformatics.babraham.ac.uk/projects/fastqc/))
 - Trimming of reads ([Cutadapt](https://journal.embnet.org/index.php/embnetjournal/article/view/200))
@@ -58,6 +58,8 @@ First, you need to know whether the sequencing files at hand are expected to con
 
 Next, the sequencing data needs to be detailed in a sample sheet or organized in a folder (see [input documentation](https://nf-co.re/ampliseq/usage#input-specifications)).
 
+Furthermore, specify what type of sequencing data you have, e.g. Illumina, PacBio, IonTorrent, and Oxford Nanopore data (see [data type documentation](https://nf-co.re/ampliseq/usage#sequencing-data-types)). In the example below paired-end Illumina data is assumed, which is default.
+
 Now, you can run the pipeline using:
 
 ```bash

assets/report_template.Rmd

@@ -74,15 +74,19 @@ params:
     input_fasta: FALSE
     input_folder: FALSE
     mqc_plot: FALSE
+    porechop_abi_log_paths: FALSE
+    chopper_stats_path: FALSE
     cutadapt_summary: FALSE
+    savont_asv_table_path : FALSE
+    savont_stats_path: FALSE
     dada_filtntrim_args: FALSE
     dada_qc_f_path: FALSE
     dada_qc_r_path: ""
     dada_pp_qc_f_path: ""
     dada_pp_qc_r_path: ""
     dada_err_path: FALSE
     dada_err_run: ""
-    asv_table_path: FALSE
+    dada_asv_table_path: FALSE
     path_asv_fa: FALSE
     path_dada2_tab: FALSE
     dada_stats_path: FALSE
@@ -272,6 +276,52 @@ _The cause is likely a bug in Kaleido, a dependency of MultiQC. There is no need
 }
 ```
 
+<!-- Subsection on Porechop_ABI -->
+
+```{r, eval = !isFALSE(params$porechop_abi_log_paths), results='asis'}
+cat(paste0("
+## Adapter removal with Porechop_ABI
+
+[Porechop_ABI](https://pubmed.ncbi.nlm.nih.gov/36698762/) removed adapters from reads.
+This is important to not have any remaining non-biological remnants in the read data.
+
+Information about changes to the data can be found in the [MultiQC report](../multiqc/multiqc_report.html) and
+Porechop_ABI log files can be found in folder [porechop_abi](../porechop_abi).
+"))
+```
+
+<!-- Subsection on Chopper -->
+
+```{r, eval = !isFALSE(params$chopper_stats_path), results='asis'}
+cat(paste0("
+## Read filtering with Chopper
+
+[Chopper](https://pubmed.ncbi.nlm.nih.gov/37171891/) filtered reads based on quality cutoffs.
+Reducing the data and number of reads can speed up subsequent steps, but is not strictely required.
+"))
+
+# import tsv
+chopper_stats <- read.table(file = params$chopper_stats_path, header = TRUE, sep = "\t")
+chopper_stats$passing_percent <- chopper_stats$chopper_output / chopper_stats$chopper_input * 100
+chopper_max_discarded <- round( 100 - min(chopper_stats$passing_percent), 1 )
+chopper_avg_passed <- round(mean(chopper_stats$passing_percent),1)
+
+# prep text
+chopper_text <- paste0("Less than ",
+    chopper_max_discarded, "% of the sequences were discarded per sample and a mean of ",
+    chopper_avg_passed, "% of the sequences (",min(chopper_stats$chopper_output)," to ",
+    max(chopper_stats$chopper_output),") per sample passed the filtering.
+")
+cat(chopper_text)
+
+# display table
+datatable(chopper_stats, options = list(
+        scrollX = TRUE,
+        scrollY = "300px",
+        paging = FALSE))
+
+```
+
 <!-- Subsection on Cutadapt -->
 
 ```{r, eval = !isFALSE(params$cutadapt_summary), results='asis'}
@@ -314,6 +364,68 @@ datatable(cutadapt_summary, options = list(
         paging = FALSE))
 ```
 
+<!-- Subsection on Savont ASVs -->
+
+```{r, eval = !isFALSE(params$savont_asv_table_path) || !isFALSE(params$savont_stats_path), results='asis'}
+cat(paste0("
+# ASV inference using Savont
+
+[Savont](https://doi.org/10.64898/2026.05.26.727271) generates Amplicon Sequence Variants (ASVs) at single-nucleotide resolution
+from long-read amplicon sequencing data such as Oxford Nanopore (ONT) R10.4 sequencing (preferably with SUP basecalling) or PacBio HiFi sequencing.
+
+Savont filters reads (98% quality), produces ASVs, and removes apparent PCR chimera.
+"))
+```
+
+<!-- Subsection on Savont ASVs -->
+
+```{r, eval = !isFALSE(params$savont_asv_table_path), results='asis'}
+cat("## Inferred ASVs\n\n")
+
+#import asv table
+asv_table <- read.table(file = params$savont_asv_table_path, header = TRUE, sep = "\t")
+n_asv <- nrow(asv_table)
+n_asv_final <- nrow(asv_table) #this is to report the original number later in the methods section
+
+# Output text
+cat("Savont generated", n_asv,
+        "amplicon sequence variants (ASVs) across all samples. ")
+if ( params$dada_sample_inference == "independent" ) {
+    cat("The ASVs can be found in [`savont/pooled_final_asvs.fasta`](../savont/). And the corresponding",
+        " quantification of the ASVs across samples is in [`savont/independent_feature-table.tsv`](../savont/).",
+        "ASVs were inferred for each sample independently",
+        "and per-sample information can be found in [`savont/per-sample](../savont/per-sample/). ")
+} else if ( params$dada_sample_inference == "pooled" ) {
+    cat("The ASVs can be found in [`savont/pooled_final_asvs.fasta`](../savont/). And the corresponding",
+        " quantification of the ASVs across samples is in [`savont/pooled_feature-table.tsv`](../savont/).",
+        "ASVs were inferred from pooled sample information. ")
+}
+```
+
+<!-- Subsection on Savont per-sample counts -->
+
+```{r, eval = !isFALSE(params$savont_stats_path), results='asis'}
+cat(paste0("
+## Read counts per sample
+
+The following table shows the read numbers per sample after ASV inference.
+"))
+
+# import tsv
+savont_stats <- read.table(file = params$savont_stats_path, header = TRUE, sep = "\t")
+
+datatable(savont_stats, options = list(
+        scrollX = TRUE,
+        scrollY = "300px",
+        paging = FALSE))
+
+cat(paste0("
+The ASV count table contained in total ",sum(savont_stats$savont_output)," counts,
+at least ",min(savont_stats$savont_output)," and at most ",max(savont_stats$savont_output)," per sample
+(average ",round(sum(savont_stats$savont_output)/length(savont_stats$savont_output),0),").
+"))
+```
+
 <!-- Subsection on DADA2 QC filtering -->
 
 ```{r, eval = !isFALSE(params$dada_filtntrim_args), results='asis'}
@@ -401,9 +513,9 @@ that position. Original plots can be found in folder [dada2/QC/](../dada2/QC/) w
 "))
 ```
 
-<!-- Section on sequences / DADA2 ASVs -->
+<!-- Section on DADA2 ASVs -->
 
-```{r, eval = !isFALSE(params$dada_err_path) || !isFALSE(params$dada_stats_path) || !isFALSE(params$asv_table_path), results='asis'}
+```{r, eval = !isFALSE(params$dada_err_path) || !isFALSE(params$dada_stats_path) || !isFALSE(params$dada_asv_table_path), results='asis'}
 cat(paste0("
 # ASV inference using DADA2
 
@@ -595,13 +707,13 @@ Samples that are very different in lost reads (per stage) to the majority of sam
 
 <!-- Subsection on DADA2 ASVs -->
 
-```{r, eval = !isFALSE(params$asv_table_path), results='asis'}
+```{r, eval = !isFALSE(params$dada_asv_table_path), results='asis'}
 cat("## Inferred ASVs\n\n")
 
 #import asv table
-asv_table <- read.table(file = params$asv_table_path, header = TRUE, sep = "\t")
+asv_table <- read.table(file = params$dada_asv_table_path, header = TRUE, sep = "\t")
 n_asv <- length(asv_table$ASV_ID)
-n_asv_dada <- length(asv_table$ASV_ID) #this is to report the original number later in the methods section
+n_asv_final <- length(asv_table$ASV_ID) #this is to report the original number later in the methods section
 
 # Output text
 cat("Finally,", n_asv,
@@ -638,7 +750,7 @@ n_asv_vsearch_cluster <- nrow(vsearch_cluster)
 cat(paste0("
 ## Clustering of ASVs
 
-[VSEARCH](https://peerj.com/articles/2584/) clustered ",n_asv_dada," ASVs into ",n_asv_vsearch_cluster,"
+[VSEARCH](https://peerj.com/articles/2584/) clustered ",n_asv_final," ASVs into ",n_asv_vsearch_cluster,"
 centroids with pairwise identity of ",params$vsearch_cluster_id,".
 Clustered ASV sequences and abundances can be found in folder [vsearch_cluster](../vsearch_cluster).
 "))
@@ -1893,6 +2005,13 @@ Data quality was evaluated with FastQC ([Andrews, 2010](https://www.bioinformati
 and summarized with MultiQC ([Ewels et al., 2016](https://pubmed.ncbi.nlm.nih.gov/27312411/)).
 "))
 ```
+```{r, eval = !isFALSE(params$porechop_abi_log_paths), results='asis'}
+cat("[Porechop_ABI](https://pubmed.ncbi.nlm.nih.gov/36698762/) removed adapters from reads. ")
+```
+```{r, eval = !isFALSE(params$chopper_stats_path), results='asis'}
+cat("[Chopper](https://pubmed.ncbi.nlm.nih.gov/37171891/) filtered reads based on quality cutoffs: ")
+cat(chopper_text)
+```
 ```{r, eval = !isFALSE(params$cutadapt_summary), results='asis'}
 cutadapt_intro = "Cutadapt ([Marcel et al., 2011](https://doi.org/10.14806/ej.17.1.200)) trimmed primers"
 if ( isFALSE(params$flag_retain_untrimmed) ) {
@@ -1901,6 +2020,15 @@ if ( isFALSE(params$flag_retain_untrimmed) ) {
     cat(paste0(cutadapt_intro, "."))
 }
 ```
+```{r, eval = !isFALSE(params$savont_asv_table_path), results='asis'}
+cat(paste0("
+[Savont](https://doi.org/10.64898/2026.05.26.727271) filtered reads (98% quality), produced ASVs, and removed PCR chimera.
+Ultimately, ",n_asv_final," amplicon sequencing variants (ASVs) were obtained across all samples.
+The ASV count table contained in total ",sum(savont_stats$savont_output)," counts,
+at least ",min(savont_stats$savont_output)," and at most ",max(savont_stats$savont_output)," per sample
+(average ",round(sum(savont_stats$savont_output)/length(savont_stats$savont_output),0),").
+"))
+```
 ```{r, eval = !isFALSE(params$dada_filtntrim_args), results='asis'}
 if ( !isFALSE(params$cutadapt_summary) && isFALSE(params$flag_retain_untrimmed) ) {
     cat("Adapter and primer-free sequences were processed ")
@@ -1952,7 +2080,7 @@ if ( isFALSE(params$flag_single_end) ) {
 }
 cat(paste0("
 and remove polymerase chain reaction (PCR) chimeras;
-ultimately, ",n_asv_dada," amplicon sequencing variants (ASVs) were obtained across all samples.
+ultimately, ",n_asv_final," amplicon sequencing variants (ASVs) were obtained across all samples.
 Between ",min(dada_stats_p$analysis),"% and ",max(dada_stats_p$analysis),"% reads per sample
 (average ",dada_stats_p_analysis_average,"%) were retained.
 The ASV count table contained in total ",sum(dada_stats_ex$analysis)," counts,
@@ -1963,7 +2091,7 @@ at least ",min(dada_stats_ex$analysis)," and at most ",max(dada_stats_ex$analysi
 
 ```{r, eval = !isFALSE(params$vsearch_cluster), results='asis'}
 cat(paste0("
-VSEARCH ([Rognes et al., 2016](https://peerj.com/articles/2584/)) clustered ",n_asv_dada," ASVs into ",n_asv_vsearch_cluster,"
+VSEARCH ([Rognes et al., 2016](https://peerj.com/articles/2584/)) clustered ",n_asv_final," ASVs into ",n_asv_vsearch_cluster,"
 centroids with pairwise identity of ",params$vsearch_cluster_id,"."))
 ```
 ```{r, eval = (params$decontam != "none"), results='asis'}

bin/savont_asv_stats.sh

@@ -0,0 +1,25 @@
+#!/bin/bash
+
+awk '
+    NR == 1 {
+        # Extract sample names from header (skip first column)
+        for (i = 2; i <= NF; i++) {
+            samples[i-1] = $i
+        }
+        num_samples = NF - 1
+        next
+    }
+    {
+        # Accumulate counts for each sample
+        for (i = 2; i <= NF; i++) {
+            totals[samples[i-1]] += $i
+        }
+    }
+    END {
+        print "sample\tsavont_output"
+        # Print in order of first appearance
+        for (i = 1; i <= num_samples; i++) {
+            print samples[i] "\t" totals[samples[i]]
+        }
+    }
+' "$1"