adding the capacity for decontamination of arbitrary reference datasets with deacon

Author: nrminor

.gitignore

@@ -35,8 +35,14 @@
 !/conf/*.yml
 
 # bin of executable scripts
-!/bin
+!/bin/
 !/bin/*.py
+!/bin/*.R
+!/bin/*.r
+!/bin/*.pl
+!/bin/*.lua
+!/bin/*.sh
+!/bin/*.awk
 
 # groovy libraries
 !/lib

bin/make_primer_patterns.lua

@@ -0,0 +1,115 @@
+#!/usr/bin/env lua
+
+-- Usage:
+-- lua make_primer_patterns.lua -i input.fasta [-o output_prefix] [-f forward_pattern] [-r reverse_pattern]
+
+local function parse_args()
+	local opts = {
+		input_fasta = nil,
+		output_prefix = "primer_patterns",
+		forward_pattern = "^(.*?)",
+		reverse_pattern = "^(.*?)",
+	}
+
+	local i = 1
+	while i <= #arg do
+		local a = arg[i]
+		if a == "-i" or a == "--input_fasta" then
+			i = i + 1
+			opts.input_fasta = arg[i]
+		elseif a == "-o" or a == "--output_prefix" then
+			i = i + 1
+			opts.output_prefix = arg[i]
+		elseif a == "-f" or a == "--forward_pattern" then
+			i = i + 1
+			opts.forward_pattern = arg[i]
+		elseif a == "-r" or a == "--reverse_pattern" then
+			i = i + 1
+			opts.reverse_pattern = arg[i]
+		else
+			error("Unknown argument: " .. a)
+		end
+		i = i + 1
+	end
+
+	assert(opts.input_fasta, "You must provide --input_fasta (-i)")
+	return opts
+end
+
+local function read_fasta_lines(path)
+	local file, err = io.open(path, "r")
+	assert(file, "Could not open FASTA file: " .. (err or "unknown error"))
+	local lines = {}
+	for line in file:lines() do
+		lines[#lines + 1] = line:gsub("%s+", "") -- strip whitespace
+	end
+	file:close()
+	return lines
+end
+
+local function generate_patterns(fasta_path, label, fwd_prefix, rev_suffix)
+	local lines = read_fasta_lines(fasta_path)
+
+	-- Extract sequences and headers
+	local seqs, headers, entries = {}, {}, {}
+	local accumulator = nil
+	for _, line in ipairs(lines) do
+		if line:sub(1, 1) == ">" then
+			if accumulator then
+				entries[#entries + 1] = accumulator
+			end
+			accumulator = { header = line, seq = "" }
+		else
+			assert(accumulator)
+			accumulator.seq = accumulator.seq .. line
+		end
+	end
+	if accumulator then
+		entries[#entries + 1] = accumulator
+	end
+
+	-- Crash if the number of parsed sequences isn't exactly 2
+	assert(#seqs == 2)
+
+	-- Parse start coordinates from header lines
+	local starts = {}
+	for _, header in ipairs(headers) do
+		local start = header:match(":(%d+)%-%d+")
+		assert(start, "Invalid header format (expected bedtools-style): " .. header)
+		starts[#starts + 1] = tonumber(start)
+	end
+
+	-- Heuristic check for orientation assumption
+	if starts[1] > starts[2] then
+		io.stderr:write(
+			"⚠️  Warning: Please double check that the provided FASTA is formatted like an output from `bedtools getfasta`, e.g.\n\n'>PP599462.1:0-16'"
+		)
+	end
+
+	-- Build patterns
+	local fwd_pattern = fwd_prefix .. seqs[1]
+	local rev_pattern = seqs[2] .. rev_suffix
+
+	-- Write output
+	local out, err = io.open(label .. ".txt", "w")
+	assert(out, "Failed to write output: " .. (err or "unknown error"))
+	out:write(fwd_pattern .. "\n")
+	out:write(rev_pattern .. "\n")
+end
+
+local function main()
+	local opts = parse_args()
+
+	-- make sure the provided file exists
+	local file = io.open(opts.input_fasta, "r")
+	assert(file, "Input FASTA file does not exist: " .. opts.input_fasta)
+	file:close()
+
+	-- generate the patterns and write them to a text file
+	generate_patterns(opts.input_fasta, opts.output_prefix, opts.forward_pattern, opts.reverse_pattern)
+end
+
+-- Run main
+if debug.getinfo(1, "S").short_src == arg[0] then
+	main()
+end

bin/split_primer_combos.lua

@@ -0,0 +1,42 @@
+#!/usr/bin/env lua
+
+-- Usage: lua split_primers.lua input.bed _LEFT _RIGHT
+-- Default suffixes if not provided
+local input_bed = arg[1] -- remember lua is 1-based!
+local forward_suffix = arg[2] or "_LEFT"
+local reverse_suffix = arg[3] or "_RIGHT"
+
+-- Make sure the input bed file is provided
+assert(input_bed, "Usage: lua split_primers.lua <input.bed> [_LEFT] [_RIGHT]")
+
+-- Initialize a table to store each primer
+local primers = {}
+
+-- Read BED line by line
+local file = io.open(input_bed, "r")
+assert(file, "Failed to open input BED file.")
+
+for line in file:lines() do
+	local ref, start_pos, stop_pos, name, index, sense =
+		line:match("([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)")
+	assert(name, "Malformed BED line: " .. line)
+
+	local base_name = name
+	base_name = base_name:gsub(forward_suffix, "")
+	base_name = base_name:gsub(reverse_suffix, "")
+
+	primers[base_name] = primers[base_name] or {}
+	table.insert(primers[base_name], line)
+end
+
+file:close()
+
+-- Write each group to its own BED file
+for base_name, records in pairs(primers) do
+	assert(#records == 2, "Expected 2 records for " .. base_name .. ", found " .. #records)
+	local out = assert(io.open(base_name .. ".bed", "w"))
+	for _, rec in ipairs(records) do
+		out:write(rec .. "\n")
+	end
+	out:close()
+end

main.nf

@@ -63,6 +63,10 @@ workflow {
         Channel.fromPath( params.ref_gbk ) :
         Channel.empty()
 
+    ch_contam_fasta = params.contam_fasta && file(params.contam_fasta).isFile()
+        ? Channel.fromPath( params.contam_fasta )
+        : Channel.empty()
+
     ch_snpeff_config = params.snpEff_config ?
         Channel.fromPath( params.snpEff_config ) :
         Channel.empty()
@@ -74,6 +78,7 @@ workflow {
             ch_primer_bed,
             ch_refseq,
             ch_ref_gbk,
+            ch_contam_fasta,
             ch_snpeff_config,
         )
 

modules/deacon.nf

@@ -0,0 +1,49 @@
+process INDEX_CONTAMINANTS {
+
+    /* */
+
+	errorStrategy { task.attempt < 3 ? 'retry' : 'ignore' }
+	maxRetries 2
+
+	cpus 1
+
+    input:
+    path contam_fasta
+
+    output:
+    path "*.idx"
+
+    script:
+    def dbName = file(contam_fasta).getSimpleName()
+    """
+    deacon index build ${contam_fasta} > ${dbName}.idx
+    """
+
+}
+
+process DECONTAMINATE {
+
+    /* */
+
+	tag "${barcode}"
+
+	errorStrategy { task.attempt < 3 ? 'retry' : 'ignore' }
+	maxRetries 2
+
+	cpus 1
+
+    input:
+	tuple val(barcode), path(reads), path(contam_index)
+
+    output:
+	tuple val(barcode), path("${barcode}.decontam.fastq.gz")
+
+    script:
+    def dbName = file(contam_index).getSimpleName()
+    """
+    zcat ${reads} \
+    | deacon filter ${contam_index} --log ${dbName}.log.json \
+    | pigz > ${barcode}.decontam.fastq.gz
+    """
+
+}

modules/fastqc.nf

@@ -15,6 +15,7 @@ process FASTQC {
 	output:
     path "*.html", emit: html
     path "*.zip", emit: zip
+    tuple val(barcode), path(reads), emit: fastq
 
 	script:
 	"""
@@ -49,4 +50,4 @@ process FASTQC_RS {
 	mv fastqc_data.txt ${label}/fastqc_data.txt
 	"""
 
-}
+}

nextflow.config

@@ -139,6 +139,9 @@ params {
     // nextclade dataset
     nextclade_dataset = null
 
+	// Contamination FASTA dataset to scrub from reads
+	contam_fasta = null
+
 	// devider haplotype phasing preset
 	devider_preset = "nanopore-r10" // old-long-reads, nanopore-r9, nanopore-r10, hi-fi