Binned quality option

bin/dada2_learn_errors.R

@@ -3,6 +3,7 @@
 suppressWarnings(suppressMessages(library(ggplot2, quietly = TRUE)))
 suppressWarnings(suppressMessages(library(gridExtra, quietly = TRUE)))
 suppressWarnings(suppressMessages(library(argparse, quietly = TRUE)))
+suppressWarnings(suppressMessages(library(jsonlite, quietly = TRUE)))
 suppressWarnings(suppressMessages(library(dada2, quietly = TRUE)))
 
 
@@ -13,10 +14,25 @@ errplot <- function(err, title=""){
 }
 
 
+learn_errors_args <- function(fls, multithread, errorEstimationFunction=NULL){
+  args <- list(fls=fls, multithread=multithread)
+
+  if(!is.null(errorEstimationFunction)){
+    args$errorEstimationFunction <- errorEstimationFunction
+  }
+
+  args
+}
+
+
 main <- function(arguments){
   parser <- ArgumentParser()
   parser$add_argument('--r1', help='file listing R1 fq files for this batch')
   parser$add_argument('--r2', help='file listing R2 fq files for this batch')
+  parser$add_argument('--params',
+                      help=paste(
+                          'json file containing optional parameters for ',
+                          'learnErrors (see README)'))
   parser$add_argument('--model', default='error_model.rds',
                       help='output .rds file containing the model')
   parser$add_argument('--plots',
@@ -30,9 +46,28 @@ main <- function(arguments){
   fnFs <- readLines(args$r1)
   fnRs <- readLines(args$r2)
 
+  if(is.null(args$params)){
+    params <- list()
+  }else{
+    params <- fromJSON(args$params)
+  }
+
+  binnedQs <- params$learnErrors$binnedQs
+  errorEstimationFunction <- NULL
+  if(!is.null(binnedQs) && length(binnedQs) > 0){
+    binnedQs <- as.numeric(binnedQs)
+    if(any(is.na(binnedQs))){
+      stop('learnErrors.binnedQs must contain only numeric quality values')
+    }
+    cat(gettextf('using binned quality error model with bins: %s\n',
+                 paste(binnedQs, collapse=', ')))
+    errorEstimationFunction <- dada2::makeBinnedQualErrfun(binnedQs)
+  }
+
   cat('generating error model for forward reads\n')
   errF <- tryCatch(
-      dada2::learnErrors(fnFs, multithread=multithread),
+      do.call(dada2::learnErrors,
+              learn_errors_args(fnFs, multithread, errorEstimationFunction)),
       error=function(err){
         cat('Error:', err$message, '\n')
         cat('saving NULL error model for forward reads\n')
@@ -41,7 +76,8 @@ main <- function(arguments){
 
   cat('generating error model for reverse reads\n')
   errR <- tryCatch(
-      dada2::learnErrors(fnRs, multithread=multithread),
+      do.call(dada2::learnErrors,
+              learn_errors_args(fnRs, multithread, errorEstimationFunction)),
       error=function(err){
         cat('Error:', err$message, '\n')
         cat('saving NULL error model for reverse reads\n')
@@ -65,4 +101,3 @@ main <- function(arguments){
 
 main(commandArgs(trailingOnly=TRUE))
 ## invisible(warnings())
-

docker/Dockerfile

@@ -1,6 +1,7 @@
 FROM rocker/r2u:22.04
 
-ARG DADA2_REF=v1.26
+ARG DADA2_REF=72da770 # corresponds to bioconductor 1.40.0 release and previous
+
 ENV DADA2_REF=$DADA2_REF PIP_NO_CACHE_DIR=1
 
 # https://ethicalhackx.com/speed-apt-get-update-parallel-downloads/

main.nf

@@ -267,6 +267,7 @@ process learn_errors {
 
     input:
         tuple val(sampleids), val(batch), val(orientation), path("R1_*.fastq.gz"), path("R2_*.fastq.gz")
+        path("dada_params.json")
 
     output:
         tuple val(sampleids), val(batch), val(orientation), path("error_model_${batch}_${orientation}.rds")
@@ -279,6 +280,7 @@ process learn_errors {
     non_empty_gz.sh \$(ls R1_*.fastq.gz) > R1.txt
     non_empty_gz.sh \$(ls R2_*.fastq.gz) > R2.txt
     dada2_learn_errors.R --r1 R1.txt --r2 R2.txt \
+        --params dada_params.json \
         --model error_model_${batch}_${orientation}.rds \
         --plots error_model_${batch}_${orientation}.png
     """
@@ -564,7 +566,7 @@ workflow {
 
     // squash sampleids into list and generate models by batch and orientation
     (models, ignored_model_plots) =
-        learn_errors(filtered.groupTuple(by: [1, 2]))
+        learn_errors(filtered.groupTuple(by: [1, 2]), dada_params)
     // transpose/expand out sampleids and join models back into filtered channel
     filtered = filtered.join(models.transpose(), by: [0, 1, 2]) // by: [sampleid, batch, orientation]
     (merged, r1, r2, dada_counts, overlaps, ignored_dada_rds,

test/i100_takara/base-files.sha256

@@ -0,0 +1,7 @@
+c1a182b714919c455a6836aec0d9bc746ace691bad6008c1d906fab58bd5f57c  counts.csv
+173936440bde704a29833e57fbc62d14777dde6a3d90c2ebcf14882566f6dc47  seqs.fasta
+ab62790858a9120aae700798cf43ff569a04d69f47bf8abbbaf745a5bff2dd8b  specimen_map.csv
+11caa1c0e03c7d2dec797607732064f22ae7a23f334cb484021c667694da2569  specimen_table.csv
+1755cfff7399b1eb5e7481f9c094e7d7ee0c71d81018fd5f698470d0726285ed  sv_table.csv
+87fc1af988e490de3999213f29e87474a80440dc43f7e61e2159af30a2fb1114  sv_table_long.csv
+775812ec7619a3e46d1b0804e18bfac7d90ab2cad548501dd6eb9c9164433af4  weights.csv

test/i100_takara/dada_params.json

@@ -0,0 +1,22 @@
+{
+  "fastqPairedFilter": {
+    "maxN": 0,
+    "truncQ": 2,
+    "minLen": 100
+  },
+  "learnErrors": {
+    "binnedQs": [2, 9, 23, 38]
+  },
+  "dada": {
+    "selfConsist": true,
+    "BAND_SIZE": 4,
+    "OMEGA_A": 1e-100,
+    "OMEGA_C": 1e-40
+  },
+  "mergePairs": {
+    "maxMismatch": 0
+  },
+  "removeBimeraDenovo": {
+    "minFoldParentOverAbundance": 2
+  }
+}

test/i100_takara/fastq-list.txt

@@ -0,0 +1,4 @@
+test/i100_takara/fastq/26R150-U043_S1_L001_I1_001.fastq.gz
+test/i100_takara/fastq/26R150-U043_S1_L001_I2_001.fastq.gz
+test/i100_takara/fastq/26R150-U043_S1_L001_R1_001.fastq.gz
+test/i100_takara/fastq/26R150-U043_S1_L001_R2_001.fastq.gz

test/i100_takara/params.json

@@ -0,0 +1,13 @@
+{
+  "sample_information": "test/i100_takara/sample-information.csv",
+  "fastq_list": "test/i100_takara/fastq-list.txt",
+  "output": "output-i100_takara",
+  "index_file_type": "dual",
+  "dada_params": "test/i100_takara/dada_params.json",
+  "bidirectional": false,
+  "alignment": {
+    "library": "",
+    "model": "data/SSU_rRNA_bacteria.cm",
+    "strategy": "cmsearch"
+  }
+}

test/i100_takara/sample-information.csv

@@ -0,0 +1,2 @@
+sampleid,sample_name,n_index,s_index,project,batch,controls,datadir,R1,R2,I1,I2
+26R150-U043,26R150-U043,,,i100_takara,1,,test/i100_takara/fastq,26R150-U043_S1_L001_R1_001.fastq.gz,26R150-U043_S1_L001_R2_001.fastq.gz,26R150-U043_S1_L001_I1_001.fastq.gz,26R150-U043_S1_L001_I2_001.fastq.gz