diff --git a/bin/dada2_learn_errors.R b/bin/dada2_learn_errors.R index f2b18ab..c071d0d 100755 --- a/bin/dada2_learn_errors.R +++ b/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()) - diff --git a/docker/Dockerfile b/docker/Dockerfile index 5bfc0b6..ffc5c10 100644 --- a/docker/Dockerfile +++ b/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/ diff --git a/main.nf b/main.nf index 7b34a28..50f42c3 100644 --- a/main.nf +++ b/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, diff --git a/test/i100_takara/base-files.sha256 b/test/i100_takara/base-files.sha256 new file mode 100644 index 0000000..450940e --- /dev/null +++ b/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 diff --git a/test/i100_takara/dada_params.json b/test/i100_takara/dada_params.json new file mode 100644 index 0000000..6ebae68 --- /dev/null +++ b/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 + } +} diff --git a/test/i100_takara/fastq-list.txt b/test/i100_takara/fastq-list.txt new file mode 100644 index 0000000..b7dfcf1 --- /dev/null +++ b/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 diff --git a/test/i100_takara/fastq/26R150-U043_S1_L001_I1_001.fastq.gz b/test/i100_takara/fastq/26R150-U043_S1_L001_I1_001.fastq.gz new file mode 100644 index 0000000..f1dc9ab Binary files /dev/null and b/test/i100_takara/fastq/26R150-U043_S1_L001_I1_001.fastq.gz differ diff --git a/test/i100_takara/fastq/26R150-U043_S1_L001_I2_001.fastq.gz b/test/i100_takara/fastq/26R150-U043_S1_L001_I2_001.fastq.gz new file mode 100644 index 0000000..605a112 Binary files /dev/null and b/test/i100_takara/fastq/26R150-U043_S1_L001_I2_001.fastq.gz differ diff --git a/test/i100_takara/fastq/26R150-U043_S1_L001_R1_001.fastq.gz b/test/i100_takara/fastq/26R150-U043_S1_L001_R1_001.fastq.gz new file mode 100644 index 0000000..1d4b524 Binary files /dev/null and b/test/i100_takara/fastq/26R150-U043_S1_L001_R1_001.fastq.gz differ diff --git a/test/i100_takara/fastq/26R150-U043_S1_L001_R2_001.fastq.gz b/test/i100_takara/fastq/26R150-U043_S1_L001_R2_001.fastq.gz new file mode 100644 index 0000000..0b9463c Binary files /dev/null and b/test/i100_takara/fastq/26R150-U043_S1_L001_R2_001.fastq.gz differ diff --git a/test/i100_takara/params.json b/test/i100_takara/params.json new file mode 100644 index 0000000..cfb24c3 --- /dev/null +++ b/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" + } +} diff --git a/test/i100_takara/sample-information.csv b/test/i100_takara/sample-information.csv new file mode 100644 index 0000000..c3000ac --- /dev/null +++ b/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