Add truncq parameter to support DADA2 quality score-based read truncation

Author: zmudge3

CHANGELOG.md

@@ -6,6 +6,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## nf-core/ampliseq version 2.15.0dev
 
 ### `Added`
+- [#909](https://github.com/nf-core/ampliseq/pull/909) - Parameter `--truncq` allows read truncation by quality score, where each read is truncated at the first instance of a quality score less than or equal to `truncq` (default value is 2). This quality score-based truncation occurs before read length truncation. If `--trunc_qmin` and `--trunc_rmin` are used to automatically calculate the values for `trunclenf` and `trunclenr` used for read length truncation, those calculations are based on the read metrics before quality score-based truncation (using `truncq`) is performed. `truncq` is passed directly as `truncQ` into DADA2's [filterAndTrim](https://rdrr.io/bioc/dada2/man/filterAndTrim.html) method.
 
 ### `Changed`
 

assets/report_template.Rmd

@@ -30,6 +30,7 @@ params:
     barplot: FALSE
     abundance_tables: FALSE
     alpha_rarefaction: FALSE
+    truncq: ""
     trunclenf: ""
     trunclenr: ""
     max_ee: ""
@@ -299,6 +300,9 @@ Additional quality filtering can improve sequence recovery.
 Often it is advised trimming the last few nucleotides to avoid less well-controlled errors that can arise there.
 "))
 
+if (params$truncq) {
+    cat("Each read was truncated at the first instance of a quality score less than or equal to ", params$truncq, ". ", sep = "")
+}
 if (params$trunc_qmin) {
     f_and_tr_args <- readLines(params$dada_filtntrim_args)
     trunc_len <- strsplit(gsub(".*truncLen = c\\((.+)\\),maxN.*", "\\1",
@@ -313,7 +317,7 @@ if (params$trunc_qmin) {
         "forward reads at ", tr_len_f, " bp and reverse ",
         "reads at ", tr_len_r, " bp, reads shorter than this were discarded. ", sep = "")
 } else if (params$trunclenf == "null" && params$trunclenr == "null") {
-    cat("Reads were not trimmed. ")
+    cat("Reads were not trimmed to a specific length. ")
 } else if (params$trunclenf != 0 && params$trunclenr != 0) {
     cat("Forward reads were trimmed at ", params$trunclenf,
             " bp and reverse reads were trimmed at ", params$trunclenr,
@@ -1700,6 +1704,9 @@ if ( params$dada_sample_inference == "independent" ) {
 
 cat("with DADA2 ([Callahan et al., 2016](https://pubmed.ncbi.nlm.nih.gov/27214047/)) to eliminate PhiX contamination, ")
 
+if (params$truncq) {
+    cat("truncate reads at the first instance of a quality score less than or equal to ", params$truncq, ", ", sep = "")
+}
 if (params$trunc_qmin) {
     cat("trim reads (before median quality drops below ", params$trunc_qmin,
         " and at least ",params$trunc_rmin*100, "% of reads are retained; ",

conf/modules.config

@@ -139,8 +139,9 @@ process {
     max_len = params.max_len ?: "Inf"
     withName: DADA2_FILTNTRIM {
         ext.args = [
-            'maxN = 0, truncQ = 2, trimRight = 0, minQ = 0, rm.lowcomplex = 0, orient.fwd = NULL, matchIDs = FALSE, id.sep = "\\\\s", id.field = NULL, n = 1e+05, OMP = TRUE',
+            'maxN = 0, trimRight = 0, minQ = 0, rm.lowcomplex = 0, orient.fwd = NULL, matchIDs = FALSE, id.sep = "\\\\s", id.field = NULL, n = 1e+05, OMP = TRUE',
             "qualityType = \"${params.quality_type}\"",
+            "truncQ = ${params.truncq}",
             params.pacbio || params.iontorrent || params.single_end ? "maxEE = ${params.max_ee}" : "maxEE = c(${params.max_ee}, ${params.max_ee})",
             params.pacbio ? "trimLeft = 0, minLen = ${params.min_len}, maxLen = $max_len, rm.phix = FALSE" :
                 params.iontorrent ? "trimLeft = 15, minLen = ${params.min_len}, maxLen = $max_len, rm.phix = TRUE" :

modules/local/summary_report.nf

@@ -86,6 +86,7 @@ process SUMMARY_REPORT  {
         cutadapt_summary ?
             params.retain_untrimmed ? "flag_retain_untrimmed=TRUE,cutadapt_summary='$cutadapt_summary'" :
             "cutadapt_summary='$cutadapt_summary'" : "",
+        "truncq=$params.truncq",
         find_truncation_values ? "trunc_qmin=$params.trunc_qmin,trunc_rmin=$params.trunc_rmin" : "",
         "trunclenf='$params.trunclenf'",
         "trunclenr='$params.trunclenr'",

nextflow.config

@@ -29,6 +29,7 @@ params {
     trunc_rmin                                = 0.75
     trunclenf                                 = null
     trunclenr                                 = null
+    truncq                                    = 2
     max_ee                                    = 2
     max_len                                   = null
     ignore_failed_filtering                   = false

nextflow_schema.json

@@ -203,6 +203,13 @@
             "description": "Read trimming and quality filtering is supposed to reduce spurious results and aid error correction",
             "default": "",
             "properties": {
+                "truncq": {
+                    "type": "integer",
+                    "default": 2,
+                    "fa_icon": "fas fa-greater-than-equal",
+                    "minimum": 0,
+                    "description": "Truncate each read at the first instance of a quality score less than or equal to truncq. Applied before read length truncation."
+                },
                 "trunclenf": {
                     "type": "integer",
                     "description": "DADA2 read truncation value for forward strand, set this to 0 for no truncation",
@@ -219,7 +226,7 @@
                     "type": "integer",
                     "default": 25,
                     "description": "If --trunclenf and --trunclenr are not set, these values will be automatically determined using this median quality score",
-                    "help_text": "Automatically determine `--trunclenf` and `--trunclenr` before the median quality score drops below `--trunc_qmin`. The fraction of reads retained is defined by `--trunc_rmin`, which might override the quality cutoff.\n\nFor example:\n\n```bash\n--trunc_qmin 35\n```\n\nPlease note:\n\n1. The code choosing `--trunclenf` and `--trunclenr` using `--trunc_qmin` automatically cannot take amplicon length or overlap requirements for merging into account, therefore use with caution.\n2. A minimum value of 25 is recommended. However, high quality data with a large paired sequence overlap might justify a higher value (e.g. 35). Also, very low quality data might require a lower value.\n3. If the quality cutoff is too low to include a certain fraction of reads that is specified by `--trunc_rmin` (e.g. 0.75 means at least 75% percent of reads are retained), a lower cutoff according to `--trunc_rmin` superseeds the quality cutoff.",
+                    "help_text": "Automatically determine `--trunclenf` and `--trunclenr` before the median quality score drops below `--trunc_qmin`. The fraction of reads retained is defined by `--trunc_rmin`, which might override the quality cutoff.\n\nFor example:\n\n```bash\n--trunc_qmin 35\n```\n\nPlease note:\n\n1. The code choosing `--trunclenf` and `--trunclenr` using `--trunc_qmin` automatically cannot take amplicon length or overlap requirements for merging into account, therefore use with caution.\n2. A minimum value of 25 is recommended. However, high quality data with a large paired sequence overlap might justify a higher value (e.g. 35). Also, very low quality data might require a lower value.\n3. If the quality cutoff is too low to include a certain fraction of reads that is specified by `--trunc_rmin` (e.g. 0.75 means at least 75% percent of reads are retained), a lower cutoff according to `--trunc_rmin` superseeds the quality cutoff.\n4. The calculations for the values of `--trunclenf` and `--trunclenr` are made before any quality score-based read truncation (using `--truncq`) is performed.",
                     "fa_icon": "fas fa-greater-than-equal"
                 },
                 "trunc_rmin": {

workflows/ampliseq.nf

@@ -202,7 +202,7 @@ workflow AMPLISEQ {
     trunclenr = params.trunclenr ?: 0
     if ( !single_end && !params.illumina_pe_its && (params.trunclenf == null || params.trunclenr == null) && !params.input_fasta ) {
         find_truncation_values = true
-        log.warn "No DADA2 cutoffs were specified (`--trunclenf` & `--trunclenr`), therefore reads will be truncated where median quality drops below ${params.trunc_qmin} (defined by `--trunc_qmin`) but at least a fraction of ${params.trunc_rmin} (defined by `--trunc_rmin`) of the reads will be retained.\nThe chosen cutoffs do not account for required overlap for merging, therefore DADA2 might have poor merging efficiency or even fail.\n"
+        log.warn "No DADA2 read truncation cutoffs were specified (`--trunclenf` & `--trunclenr`), therefore reads will be truncated where median quality drops below ${params.trunc_qmin} (defined by `--trunc_qmin`) but at least a fraction of ${params.trunc_rmin} (defined by `--trunc_rmin`) of the reads will be retained.\nThe chosen cutoffs do not account for required overlap for merging, therefore DADA2 might have poor merging efficiency or even fail.\nThe cutoffs are chosen before any quality score-based read truncation (using `--truncq`) is performed.\n"
     } else { find_truncation_values = false }
 
     // save params to values to be able to overwrite it