From bbe2ac2f22930d6503c248c9567fa76938ac6f83 Mon Sep 17 00:00:00 2001 From: Nick Minor Date: Tue, 28 Jul 2026 16:48:54 -0500 Subject: [PATCH 01/10] feat: retain BLAST hits across taxa Remove the post-BLAST query-group gates that required a non-phage viral hit after MEGABLAST and BLASTN. Target enrichment and depletion continue to control which sequences reach BLAST, while every annotated top-hit group now proceeds directly to merging, LCA annotation, and reporting. MEGABLAST partitioning and existing hit-retention limits remain unchanged. Delete the obsolete Nextflow processes, Python filter and summary scripts, unit tests, filter decision artifacts, and sequence-flow rows. Rename affected channels around annotated-hit semantics, publish annotated MEGABLAST and BLASTN TSVs at their result-family roots, and remove the virus-only result directories. Describe the existing enrichment parameters as target-enrichment controls in executable help without renaming the 3.x API. Normalize the sequence-flow process input to the codebase's scalar-or-list convention because removing the filter summaries allows valid runs with a single evidence file. --- bin/build_sequence_flow.py | 20 ----- bin/filter_non_virus_blast_nodes.py | 67 --------------- bin/finalize_blast_results.py | 2 +- bin/summarize_blast_filtering.py | 67 --------------- bin/test_filter_non_virus_blast_nodes.py | 55 ------------ bin/test_summarize_blast_filtering.py | 102 ----------------------- conf/results.config | 38 +-------- lib/py_nvd/cli/commands/preset.py | 14 ++-- lib/py_nvd/cli/commands/run.py | 14 ++-- lib/py_nvd/models.py | 18 ++-- modules/blast.nf | 84 ++----------------- modules/reporting.nf | 3 +- nextflow.config | 4 +- subworkflows/classify_with_blastn.nf | 14 +--- subworkflows/classify_with_megablast.nf | 9 +- workflows/nvd_main.nf | 4 +- 16 files changed, 42 insertions(+), 473 deletions(-) delete mode 100755 bin/filter_non_virus_blast_nodes.py delete mode 100755 bin/summarize_blast_filtering.py delete mode 100644 bin/test_filter_non_virus_blast_nodes.py delete mode 100644 bin/test_summarize_blast_filtering.py diff --git a/bin/build_sequence_flow.py b/bin/build_sequence_flow.py index 53921278..9ec62f3f 100755 --- a/bin/build_sequence_flow.py +++ b/bin/build_sequence_flow.py @@ -181,22 +181,6 @@ def megablast_partition_rows(path: Path) -> list[LedgerRow]: return rows -def blast_filter_rows(path: Path) -> list[LedgerRow]: - return [ - LedgerRow( - sample_id=source["sample_id"], - stage=source["stage"], - input_class=source["query_class"], - output_class="virus_only_hits", - decision="retain", - reason="viral_taxonomy_match", - sequences_in=source["queries_in"], - sequences_out=source["queries_retained"], - ) - for source in read_tsv(path) - ] - - def rows_for_input(path: Path) -> list[LedgerRow]: name = path.name adapters = ( @@ -206,10 +190,6 @@ def rows_for_input(path: Path) -> list[LedgerRow]: (("_mapped_counts.txt",), mapback_rows), ((".blast_query_batches.tsv",), query_batch_rows), ((".megablast_query_partition.tsv",), megablast_partition_rows), - ( - (".megablast_query_filtering.tsv", ".blastn_query_filtering.tsv"), - blast_filter_rows, - ), ) for suffixes, adapter in adapters: if name.endswith(suffixes): diff --git a/bin/filter_non_virus_blast_nodes.py b/bin/filter_non_virus_blast_nodes.py deleted file mode 100755 index 77ad4919..00000000 --- a/bin/filter_non_virus_blast_nodes.py +++ /dev/null @@ -1,67 +0,0 @@ -#!/usr/bin/env python3 -# /// script -# requires-python = ">=3.11" -# dependencies = [ -# "pandas", -# ] -# /// - -import argparse -import os - -import pandas as pd - -LINEAGE_TOKEN_PARTS = 2 - - -def parse_args() -> argparse.Namespace: - parser = argparse.ArgumentParser( - description="Keep BLAST query groups with non-phage viral hits", - ) - parser.add_argument("input_file", help="Annotated BLAST TSV") - parser.add_argument("output_file", help="Filtered output TSV") - return parser.parse_args() - - -def is_virus_lineage(lineage: object) -> bool: - """Return true when a lineage token names Viruses exactly.""" - for token in str(lineage).split(";"): - parts = token.strip().split(":", maxsplit=1) - if len(parts) == LINEAGE_TOKEN_PARTS and parts[1].strip() == "Viruses": - return True - return False - - -def contains_non_phage_viruses(group: pd.DataFrame) -> bool: - virus_hits = group[group["rank"].apply(is_virus_lineage)] - non_phage_viruses = virus_hits[ - ~virus_hits["stitle"].str.contains("phage", case=False, na=False) - ] - return len(non_phage_viruses) > 0 - - -def filter_non_virus_blast_nodes(input_file: str, output_file: str) -> None: - - # BLAST processes may emit zero-byte files when no hits are found. - if os.path.getsize(input_file) == 0: - # Create an empty output file if the input is empty - open(output_file, "w").close() - return - - # Read the input file - all_hits_df = pd.read_csv(input_file, sep="\t") - - # Filter the dataframe - filtered_df = all_hits_df.groupby("qseqid").filter(contains_non_phage_viruses) - - # Write the filtered data to the output file - filtered_df.to_csv(output_file, sep="\t", index=False) - - -def main() -> None: - args = parse_args() - filter_non_virus_blast_nodes(args.input_file, args.output_file) - - -if __name__ == "__main__": - main() diff --git a/bin/finalize_blast_results.py b/bin/finalize_blast_results.py index aa2261b6..fdee9271 100755 --- a/bin/finalize_blast_results.py +++ b/bin/finalize_blast_results.py @@ -98,7 +98,7 @@ def main(argv: list[str] | None = None) -> None: parser.add_argument( "--virus-index-version", required=True, - help="STAT k-mer database version used for virus enrichment index", + help="Target-enrichment index version", ) parser.add_argument( "--run-id", diff --git a/bin/summarize_blast_filtering.py b/bin/summarize_blast_filtering.py deleted file mode 100755 index 930e2d90..00000000 --- a/bin/summarize_blast_filtering.py +++ /dev/null @@ -1,67 +0,0 @@ -#!/usr/bin/env python3 -"""Summarize query retention across a BLAST taxonomic filter.""" - -from __future__ import annotations - -import argparse -import csv -from pathlib import Path - - -def parse_args() -> argparse.Namespace: - parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument("--sample-id", required=True) - parser.add_argument("--query-class", required=True) - parser.add_argument("--stage", required=True) - parser.add_argument("--input", type=Path, required=True) - parser.add_argument("--output", type=Path, required=True) - parser.add_argument("--summary", type=Path, required=True) - return parser.parse_args() - - -def query_ids(path: Path) -> set[str]: - if path.stat().st_size == 0: - return set() - - with path.open(newline="", encoding="utf-8") as handle: - reader = csv.DictReader(handle, delimiter="\t") - if reader.fieldnames is None or "qseqid" not in reader.fieldnames: - message = ( - f"annotated BLAST TSV {path.name!r} is missing required column 'qseqid'" - ) - raise ValueError(message) - return {row["qseqid"] for row in reader} - - -def main() -> None: - args = parse_args() - input_ids = query_ids(args.input) - output_ids = query_ids(args.output) - with args.summary.open("w", newline="", encoding="utf-8") as handle: - writer = csv.DictWriter( - handle, - fieldnames=( - "sample_id", - "query_class", - "stage", - "queries_in", - "queries_retained", - "queries_removed", - ), - delimiter="\t", - ) - writer.writeheader() - writer.writerow( - { - "sample_id": args.sample_id, - "query_class": args.query_class, - "stage": args.stage, - "queries_in": len(input_ids), - "queries_retained": len(output_ids), - "queries_removed": len(input_ids - output_ids), - }, - ) - - -if __name__ == "__main__": - main() diff --git a/bin/test_filter_non_virus_blast_nodes.py b/bin/test_filter_non_virus_blast_nodes.py deleted file mode 100644 index 4f098491..00000000 --- a/bin/test_filter_non_virus_blast_nodes.py +++ /dev/null @@ -1,55 +0,0 @@ -"""Tests for filtering annotated BLAST rows to viral query groups.""" - -from __future__ import annotations - -import pandas as pd -from filter_non_virus_blast_nodes import contains_non_phage_viruses - - -def group_with(*, rank: str, stitle: str = "viral reference") -> pd.DataFrame: - """Build a minimal qseqid group for filter predicate tests.""" - return pd.DataFrame( - { - "qseqid": ["contig1"], - "rank": [rank], - "stitle": [stitle], - }, - ) - - -def test_superkingdom_viruses_lineage_is_kept() -> None: - """Legacy superkingdom lineage strings are still viral.""" - group = group_with(rank="root:cellular organisms; superkingdom:Viruses") - - assert contains_non_phage_viruses(group) - - -def test_acellular_root_viruses_lineage_is_kept() -> None: - """Current taxonomy strings can identify Viruses at acellular root.""" - group = group_with(rank="acellular root:Viruses; realm:Riboviria") - - assert contains_non_phage_viruses(group) - - -def test_non_viral_lineage_is_rejected() -> None: - """Non-viral taxonomy strings should not keep the query group.""" - group = group_with(rank="domain:Bacteria; phylum:Pseudomonadota") - - assert not contains_non_phage_viruses(group) - - -def test_viral_phage_only_group_is_rejected() -> None: - """Phage-only viral query groups keep the existing exclusion behavior.""" - group = group_with( - rank="acellular root:Viruses; realm:Duplodnaviria", - stitle="Escherichia phage lambda", - ) - - assert not contains_non_phage_viruses(group) - - -def test_substring_lookalike_is_not_viral() -> None: - """The viral predicate should match exact token names, not substrings.""" - group = group_with(rank="note:NotVirusesMaybe; domain:Bacteria") - - assert not contains_non_phage_viruses(group) diff --git a/bin/test_summarize_blast_filtering.py b/bin/test_summarize_blast_filtering.py deleted file mode 100644 index ae696c89..00000000 --- a/bin/test_summarize_blast_filtering.py +++ /dev/null @@ -1,102 +0,0 @@ -"""Tests for BLAST query-filtering summaries.""" - -from __future__ import annotations - -import csv -import sys -from typing import TYPE_CHECKING -from unittest.mock import patch - -import pytest -from summarize_blast_filtering import main, query_ids - -if TYPE_CHECKING: - from pathlib import Path - - -def test_query_ids_accepts_zero_byte_and_header_only_inputs(tmp_path: Path) -> None: - empty = tmp_path / "empty.tsv" - empty.touch() - header_only = tmp_path / "header_only.tsv" - header_only.write_text("sample\tqseqid\ttask\n", encoding="utf-8") - - assert query_ids(empty) == set() - assert query_ids(header_only) == set() - - -def test_query_ids_resolves_reordered_qseqid_and_deduplicates_hits( - tmp_path: Path, -) -> None: - annotated = tmp_path / "annotated.tsv" - annotated.write_text( - "sample\ttask\tqseqid\n" - "sample_A\tmegablast\tquery_001\n" - "sample_A\tmegablast\tquery_001\n" - "sample_A\tmegablast\tquery_002\n", - encoding="utf-8", - ) - - assert query_ids(annotated) == {"query_001", "query_002"} - - -def test_query_ids_rejects_header_without_qseqid(tmp_path: Path) -> None: - annotated = tmp_path / "annotated.tsv" - annotated.write_text( - "sample\ttask\nsample_A\tmegablast\n", - encoding="utf-8", - ) - - with pytest.raises(ValueError, match="qseqid"): - query_ids(annotated) - - -def test_summary_counts_queries_by_named_qseqid(tmp_path: Path) -> None: - annotated = tmp_path / "annotated.tsv" - annotated.write_text( - "sample\tqseqid\ttask\n" - "sample_A\tquery_001\tmegablast\n" - "sample_A\tquery_001\tmegablast\n" - "sample_A\tquery_002\tmegablast\n", - encoding="utf-8", - ) - retained = tmp_path / "retained.tsv" - retained.write_text( - "task\tsample\tqseqid\nmegablast\tsample_A\tquery_001\n", - encoding="utf-8", - ) - summary = tmp_path / "summary.tsv" - - with patch.object( - sys, - "argv", - [ - "summarize_blast_filtering.py", - "--sample-id", - "sample_A", - "--query-class", - "single_read", - "--stage", - "megablast_virus_filter", - "--input", - str(annotated), - "--output", - str(retained), - "--summary", - str(summary), - ], - ): - main() - - with summary.open(newline="", encoding="utf-8") as handle: - rows = list(csv.DictReader(handle, delimiter="\t")) - - assert rows == [ - { - "sample_id": "sample_A", - "query_class": "single_read", - "stage": "megablast_virus_filter", - "queries_in": "2", - "queries_retained": "1", - "queries_removed": "1", - }, - ] diff --git a/conf/results.config b/conf/results.config index e56a55cf..7f8bc390 100644 --- a/conf/results.config +++ b/conf/results.config @@ -449,30 +449,13 @@ process { // ========================================================================= withName: 'ANNOTATE_MEGABLAST_RESULTS' { publishDir = [ - path: { params.megablast_hits + "/annotated" }, + path: { params.megablast_hits }, mode: 'copy', saveAs: { filename -> filename.equals('versions.yml') ? null : filename }, enabled: true ] } - withName: 'FILTER_NON_VIRUS_MEGABLAST_NODES' { - publishDir = [ - [ - path: { params.megablast_hits + "/virus_only" }, - mode: 'copy', - pattern: '*.txt', - enabled: true - ], - [ - path: { params.megablast_hits + "/decisions" }, - mode: 'copy', - pattern: '*.megablast_query_filtering.tsv', - enabled: true - ] - ] - } - withName: 'PARTITION_MEGABLAST_QUERIES' { publishDir = [ [ @@ -495,30 +478,13 @@ process { // ========================================================================= withName: 'ANNOTATE_BLASTN_RESULTS' { publishDir = [ - path: { params.blastn_hits + "/annotated" }, + path: { params.blastn_hits }, mode: 'copy', saveAs: { filename -> filename.equals('versions.yml') ? null : filename }, enabled: true ] } - withName: 'FILTER_NON_VIRUS_BLASTN_NODES' { - publishDir = [ - [ - path: { params.blastn_hits + "/virus_only" }, - mode: 'copy', - pattern: '*.txt', - enabled: true - ], - [ - path: { params.blastn_hits + "/decisions" }, - mode: 'copy', - pattern: '*.blastn_query_filtering.tsv', - enabled: true - ] - ] - } - // ========================================================================= // 07 — Merged BLAST results // ========================================================================= diff --git a/lib/py_nvd/cli/commands/preset.py b/lib/py_nvd/cli/commands/preset.py index 1b8e83b6..5dba9f45 100644 --- a/lib/py_nvd/cli/commands/preset.py +++ b/lib/py_nvd/cli/commands/preset.py @@ -199,17 +199,17 @@ def preset_register( virus_index: Path | None = typer.Option( None, "--virus-index", - help="Path to prebuilt vertebrate-infecting virus deacon index (.idx file)", + help="Path to a prebuilt Deacon target-enrichment index (.idx file)", ), virus_index_url: str | None = typer.Option( None, "--virus-index-url", - help="URL to download a prebuilt vertebrate-infecting virus deacon index", + help="URL to download a prebuilt Deacon target-enrichment index", ), virus_reference_fasta: Path | None = typer.Option( None, "--virus-reference-fasta", - help="Custom vertebrate-infecting virus FASTA for enrichment", + help="Custom target FASTA for building a Deacon enrichment index", ), no_enrichment: bool | None = typer.Option( None, @@ -219,22 +219,22 @@ def preset_register( virus_kmer_size: int | None = typer.Option( None, "--virus-kmer-size", - help="K-mer size for building a custom virus enrichment index", + help="K-mer size for building a custom target-enrichment index", ), virus_window_size: int | None = typer.Option( None, "--virus-window-size", - help="Minimizer window size for building a custom virus enrichment index", + help="Minimizer window size for building a custom target-enrichment index", ), virus_abs_threshold: int | None = typer.Option( None, "--virus-abs-threshold", - help="Minimum absolute minimizer hits for virus read enrichment", + help="Minimum absolute minimizer hits for target enrichment", ), virus_rel_threshold: float | None = typer.Option( None, "--virus-rel-threshold", - help="Minimum relative minimizer proportion for virus read enrichment", + help="Minimum relative minimizer proportion for target enrichment", ), filter_reads: bool | None = typer.Option( None, diff --git a/lib/py_nvd/cli/commands/run.py b/lib/py_nvd/cli/commands/run.py index aa75496c..1d421f04 100644 --- a/lib/py_nvd/cli/commands/run.py +++ b/lib/py_nvd/cli/commands/run.py @@ -195,19 +195,19 @@ def run( virus_index: Path | None = typer.Option( None, "--virus-index", - help="Path to prebuilt vertebrate-infecting virus deacon index (.idx file)", + help="Path to a prebuilt Deacon target-enrichment index (.idx file)", rich_help_panel=PANEL_DATABASES, ), virus_index_url: str | None = typer.Option( None, "--virus-index-url", - help="URL to download a prebuilt vertebrate-infecting virus deacon index", + help="URL to download a prebuilt Deacon target-enrichment index", rich_help_panel=PANEL_DATABASES, ), virus_reference_fasta: Path | None = typer.Option( None, "--virus-reference-fasta", - help="Custom vertebrate-infecting virus FASTA for building an enrichment index", + help="Custom target FASTA for building a Deacon enrichment index", rich_help_panel=PANEL_DATABASES, ), no_enrichment: bool | None = typer.Option( @@ -231,25 +231,25 @@ def run( virus_kmer_size: int | None = typer.Option( None, "--virus-kmer-size", - help="K-mer size for building a custom virus enrichment index (default: 31)", + help="K-mer size for building a custom target-enrichment index (default: 31)", rich_help_panel=PANEL_DATABASES, ), virus_window_size: int | None = typer.Option( None, "--virus-window-size", - help="Minimizer window size for building a custom virus enrichment index (default: 1)", + help="Minimizer window size for building a custom target-enrichment index (default: 1)", rich_help_panel=PANEL_DATABASES, ), virus_abs_threshold: int | None = typer.Option( None, "--virus-abs-threshold", - help="Minimum absolute minimizer hits for virus read enrichment (default: 1)", + help="Minimum absolute minimizer hits for target enrichment (default: 1)", rich_help_panel=PANEL_DATABASES, ), virus_rel_threshold: float | None = typer.Option( None, "--virus-rel-threshold", - help="Minimum relative minimizer proportion for virus read enrichment (default: 0.0)", + help="Minimum relative minimizer proportion for target enrichment (default: 0.0)", rich_help_panel=PANEL_DATABASES, ), # ------------------------------------------------------------------------- diff --git a/lib/py_nvd/models.py b/lib/py_nvd/models.py index b50c323d..99c7e592 100644 --- a/lib/py_nvd/models.py +++ b/lib/py_nvd/models.py @@ -177,7 +177,7 @@ def reject_removed_params(cls, data: Any) -> Any: # noqa: ANN401 # pydantic ho ) virus_index_version: str | None = Field( None, - description="Virus enrichment index version", + description="Target-enrichment index version", json_schema_extra={"category": "Databases"}, ) @@ -198,42 +198,42 @@ def reject_removed_params(cls, data: Any) -> Any: # noqa: ANN401 # pydantic ho ) virus_index: Path | None = Field( None, - description="Path to prebuilt vertebrate-infecting virus deacon index (.idx file)", + description="Path to a prebuilt Deacon target-enrichment index (.idx file)", json_schema_extra={"category": "Databases"}, ) virus_index_url: str | None = Field( None, - description="URL to download a prebuilt vertebrate-infecting virus deacon index", + description="URL to download a prebuilt Deacon target-enrichment index", json_schema_extra={"category": "Databases"}, ) virus_reference_fasta: Path | None = Field( None, - description="Custom vertebrate-infecting virus FASTA for building an enrichment index", + description="Custom target FASTA for building a Deacon enrichment index", json_schema_extra={"category": "Databases"}, ) no_enrichment: bool = Field( default=False, - description="Disable target enrichment even when a virus index source is provided", + description="Disable target enrichment even when an index source is provided", json_schema_extra={"category": "Databases"}, ) virus_kmer_size: int = Field( 31, - description="K-mer size for building a custom virus enrichment index", + description="K-mer size for building a custom target-enrichment index", json_schema_extra={"category": "Databases"}, ) virus_window_size: int = Field( 1, - description="Minimizer window size for building a custom virus enrichment index", + description="Minimizer window size for building a custom target-enrichment index", json_schema_extra={"category": "Databases"}, ) virus_abs_threshold: int = Field( 1, - description="Minimum absolute minimizer hits for virus read enrichment", + description="Minimum absolute minimizer hits for target enrichment", json_schema_extra={"category": "Databases"}, ) virus_rel_threshold: float = Field( 0.0, - description="Minimum relative proportion of minimizers for virus read enrichment (0.0-1.0)", + description="Minimum relative proportion of minimizers for target enrichment (0.0-1.0)", json_schema_extra={"category": "Databases"}, ) sourmash_ksize: int = Field( diff --git a/modules/blast.nf b/modules/blast.nf index 55cded60..83a825f5 100644 --- a/modules/blast.nf +++ b/modules/blast.nf @@ -149,7 +149,7 @@ process MEGABLAST { """ } -//Create file with megablast results annotated with megablast task and full taxonomic rank. +// Create a TSV with MEGABLAST task metadata and full taxonomic lineages. process ANNOTATE_MEGABLAST_RESULTS { tag "${sample_id}, ${query_class}" @@ -163,7 +163,7 @@ process ANNOTATE_MEGABLAST_RESULTS { val taxonomy_dir output: - tuple val(sample_id), val(query_class), path("${sample_id}.${query_class}.annotated_megablast.txt"), emit: hits + tuple val(sample_id), val(query_class), path("${sample_id}.${query_class}.annotated_megablast.tsv"), emit: hits script: def taxonomy_dir_arg = taxonomy_dir ? "--taxonomy-dir '${taxonomy_dir}'" : "" @@ -172,7 +172,7 @@ process ANNOTATE_MEGABLAST_RESULTS { """ annotate_blast_results.py \ --input_file ${blast_txt} \ - --output_file ${sample_id}.${query_class}.annotated_megablast.txt \ + --output_file ${sample_id}.${query_class}.annotated_megablast.tsv \ --sample_name ${sample_id} \ --task 'megablast' \ ${taxonomy_dir_arg} \ @@ -203,43 +203,6 @@ process SELECT_TOP_BLAST_HITS { --blast-retention-count ${params.blast_retention_count} """ } -/* - -Remove any query groups that do not have at least one hit corresponding to viruses. - -This rule handles cases where the input file is empty and filters out -blast query read groups where none of the entries are from superkingdom: Viruses -or where the only viral hits are phages. This will leave in any query seq that has things that are -non virus as long as one of the qseqid reads has a viral hit. -*/ -process FILTER_NON_VIRUS_MEGABLAST_NODES { - - tag "${sample_id}, ${query_class}" - label "low" - - errorStrategy { task.attempt < 3 ? 'retry' : 'ignore' } - maxRetries 2 - - input: - tuple val(sample_id), val(query_class), path(annotated_blast) - - output: - tuple val(sample_id), val(query_class), path("${sample_id}.${query_class}.mb_virus_only.txt"), emit: hits - tuple val(sample_id), val(query_class), path("${sample_id}.${query_class}.megablast_query_filtering.tsv"), emit: decisions - - script: - """ - filter_non_virus_blast_nodes.py ${annotated_blast} ${sample_id}.${query_class}.mb_virus_only.txt - summarize_blast_filtering.py \ - --sample-id '${sample_id}' \ - --query-class '${query_class}' \ - --stage megablast_virus_filter \ - --input ${annotated_blast} \ - --output ${sample_id}.${query_class}.mb_virus_only.txt \ - --summary ${sample_id}.${query_class}.megablast_query_filtering.tsv - """ -} - /* Partition queries according to whether MEGABLAST produced any hit. @@ -310,7 +273,7 @@ process BLASTN_CLASSIFY { """ } -// Create file with blastn results annotated with blastn task and full taxonomic rank. +// Create a TSV with BLASTN task metadata and full taxonomic lineages. process ANNOTATE_BLASTN_RESULTS { tag "${sample_id}, ${query_class}" @@ -324,7 +287,7 @@ process ANNOTATE_BLASTN_RESULTS { val taxonomy_dir output: - tuple val(sample_id), val(query_class), path("${sample_id}.${query_class}.annotated_blastn.txt") + tuple val(sample_id), val(query_class), path("${sample_id}.${query_class}.annotated_blastn.tsv") script: def taxonomy_dir_arg = taxonomy_dir ? "--taxonomy-dir '${taxonomy_dir}'" : "" @@ -334,7 +297,7 @@ process ANNOTATE_BLASTN_RESULTS { annotate_blast_results.py \ --sample_name ${sample_id} \ --input_file ${blastn_txt} \ - --output_file ${sample_id}.${query_class}.annotated_blastn.txt \ + --output_file ${sample_id}.${query_class}.annotated_blastn.tsv \ --task 'blastn' \ ${taxonomy_dir_arg} \ ${taxonomy_mode_arg} \ @@ -342,41 +305,6 @@ process ANNOTATE_BLASTN_RESULTS { """ } -/* -Remove any query groups that do not have at least one hit corresponding to viruses. - -This rule handles cases where the input file is empty and filters out -groups where none of the entries are from superkingdom: Viruses -or where the only viral hits are phages. -*/ -process FILTER_NON_VIRUS_BLASTN_NODES { - - tag "${sample_id}, ${query_class}" - label "low" - - errorStrategy { task.attempt < 3 ? 'retry' : 'ignore' } - maxRetries 2 - - input: - tuple val(sample_id), val(query_class), path(annotated_blast) - - output: - tuple val(sample_id), val(query_class), path("${sample_id}.${query_class}.nt_virus_only.txt"), emit: hits - tuple val(sample_id), val(query_class), path("${sample_id}.${query_class}.blastn_query_filtering.tsv"), emit: decisions - - script: - """ - filter_non_virus_blast_nodes.py ${annotated_blast} ${sample_id}.${query_class}.nt_virus_only.txt - summarize_blast_filtering.py \ - --sample-id '${sample_id}' \ - --query-class '${query_class}' \ - --stage blastn_virus_filter \ - --input ${annotated_blast} \ - --output ${sample_id}.${query_class}.nt_virus_only.txt \ - --summary ${sample_id}.${query_class}.blastn_query_filtering.tsv - """ -} - // Combine MEGABLAST and optional BLASTN search hits for one query-class batch. process COMBINE_BATCH_SEARCH_HITS { diff --git a/modules/reporting.nf b/modules/reporting.nf index 3f0560a2..841cb7ec 100644 --- a/modules/reporting.nf +++ b/modules/reporting.nf @@ -9,7 +9,8 @@ process BUILD_SEQUENCE_FLOW { path "sequence_flow.tsv", emit: sequence_flow script: - def input_args = input_files.collect { input -> "--input '${input}'" }.join(" ") + def files = input_files instanceof List ? input_files : [input_files] + def input_args = files.collect { input -> "--input '${input}'" }.join(" ") """ build_sequence_flow.py \ ${input_args} \ diff --git a/nextflow.config b/nextflow.config index 1892b832..c33f2634 100644 --- a/nextflow.config +++ b/nextflow.config @@ -64,8 +64,8 @@ params { blast_db = null blast_db_prefix = null - // Vertebrate-infecting virus enrichment with deacon. Enabled automatically - // when any virus index source is provided; set no_enrichment to disable it. + // Target enrichment with Deacon. Enabled automatically when any target index + // source is provided; set no_enrichment to disable it. no_enrichment = false virus_index = null virus_index_url = null diff --git a/subworkflows/classify_with_blastn.nf b/subworkflows/classify_with_blastn.nf index 9c52d623..3bc3059e 100644 --- a/subworkflows/classify_with_blastn.nf +++ b/subworkflows/classify_with_blastn.nf @@ -1,7 +1,6 @@ include { BLASTN_CLASSIFY ; ANNOTATE_BLASTN_RESULTS ; - FILTER_NON_VIRUS_BLASTN_NODES ; COMBINE_BATCH_SEARCH_HITS ; SELECT_TOP_BLAST_HITS } from "../modules/blast" @@ -9,7 +8,7 @@ include { ANNOTATE_LEAST_COMMON_ANCESTORS } from "../modules/utils" workflow CLASSIFY_WITH_BLASTN { take: - ch_filtered_megablast + ch_megablast_hits ch_megablast_query_partition ch_blast_db_files ch_taxonomy_dir // value channel: taxonomy directory path for taxonomy lookups @@ -19,7 +18,7 @@ workflow CLASSIFY_WITH_BLASTN { // this subworkflow. BLASTN is still skipped when MEGABLAST leaves no candidate // query sequences, but the skipped case remains a channel emission instead // of disappearing. - ch_samples_after_megablast = ch_filtered_megablast + ch_samples_after_megablast = ch_megablast_hits .join(ch_megablast_query_partition, by: [0, 1]) .map { sample_id, query_class, megablast_hits, accounted_query_ids, blastn_candidate_fasta, _partition_summary -> def needs_blastn = file(blastn_candidate_fasta).size() > 0 @@ -46,11 +45,7 @@ workflow CLASSIFY_WITH_BLASTN { ch_taxonomy_dir, ) - FILTER_NON_VIRUS_BLASTN_NODES( - ANNOTATE_BLASTN_RESULTS.out - ) - - ch_blastn_terminal = FILTER_NON_VIRUS_BLASTN_NODES.out.hits + ch_blastn_hits = ANNOTATE_BLASTN_RESULTS.out .join(ch_blastn_context, by: [0, 1]) .map { sample_id, query_class, blastn_hits, _meta, megablast_hits -> tuple(sample_id, query_class, [megablast_hits, blastn_hits]) @@ -62,7 +57,7 @@ workflow CLASSIFY_WITH_BLASTN { tuple(meta.id, meta.query_class, [megablast_hits]) } - ch_merged_input = ch_blastn_terminal.mix(ch_samples_skipping_blastn) + ch_merged_input = ch_blastn_hits.mix(ch_samples_skipping_blastn) COMBINE_BATCH_SEARCH_HITS(ch_merged_input) @@ -70,5 +65,4 @@ workflow CLASSIFY_WITH_BLASTN { emit: merged_results = ANNOTATE_LEAST_COMMON_ANCESTORS.out // tuple(sample_id, query_class, batch_lca_tsv) - filter_decisions = FILTER_NON_VIRUS_BLASTN_NODES.out.decisions } diff --git a/subworkflows/classify_with_megablast.nf b/subworkflows/classify_with_megablast.nf index b4875426..36ca650c 100644 --- a/subworkflows/classify_with_megablast.nf +++ b/subworkflows/classify_with_megablast.nf @@ -1,7 +1,6 @@ include { MEGABLAST ; ANNOTATE_MEGABLAST_RESULTS ; - FILTER_NON_VIRUS_MEGABLAST_NODES ; PARTITION_MEGABLAST_QUERIES ; SELECT_TOP_BLAST_HITS } from "../modules/blast" @@ -38,18 +37,12 @@ workflow CLASSIFY_WITH_MEGABLAST { ch_taxonomy_dir ) - // Capture this output for the LabKey table - FILTER_NON_VIRUS_MEGABLAST_NODES( - ANNOTATE_MEGABLAST_RESULTS.out.hits - ) - PARTITION_MEGABLAST_QUERIES( MEGABLAST.out.join(ch_megablast_candidates.for_partition, by: [0, 1]) ) emit: - filtered_megablast = FILTER_NON_VIRUS_MEGABLAST_NODES.out.hits - filter_decisions = FILTER_NON_VIRUS_MEGABLAST_NODES.out.decisions + annotated_hits = ANNOTATE_MEGABLAST_RESULTS.out.hits megablast_query_partition = PARTITION_MEGABLAST_QUERIES.out megablast = SELECT_TOP_BLAST_HITS.out } diff --git a/workflows/nvd_main.nf b/workflows/nvd_main.nf index 504a2b07..2627c77e 100644 --- a/workflows/nvd_main.nf +++ b/workflows/nvd_main.nf @@ -117,7 +117,7 @@ workflow NVD_MAIN { ) CLASSIFY_WITH_BLASTN( - CLASSIFY_WITH_MEGABLAST.out.filtered_megablast, + CLASSIFY_WITH_MEGABLAST.out.annotated_hits, CLASSIFY_WITH_MEGABLAST.out.megablast_query_partition, ch_blast_db_files, ch_taxonomy_dir, @@ -130,8 +130,6 @@ workflow NVD_MAIN { .mix(PREPARE_BLAST_QUERIES.out.mapback_count_files.map { _sample_id, counts -> counts }) .mix(PREPARE_BLAST_QUERIES.out.blast_query_summaries.map { _sample_id, summary -> summary }) .mix(CLASSIFY_WITH_MEGABLAST.out.megablast_query_partition.map { _sample_id, _query_class, _accounted_ids, _blastn_candidates, summary -> summary }) - .mix(CLASSIFY_WITH_MEGABLAST.out.filter_decisions.map { _sample_id, _query_class, decision -> decision }) - .mix(CLASSIFY_WITH_BLASTN.out.filter_decisions.map { _sample_id, _query_class, decision -> decision }) REPORTING( CLASSIFY_WITH_BLASTN.out.merged_results, From eaa8ab10d73f9ee7ad2d5b3e2fe3f5bf8f9b4d73 Mon Sep 17 00:00:00 2001 From: Nick Minor Date: Fri, 31 Jul 2026 12:23:50 -0500 Subject: [PATCH 02/10] docs: point repository override example at v3.4 Use the v3.4.0 checkout in the NVD_REPO installation example after the v3.4 release line, while preserving the existing launcher and compatibility guidance. --- docs/INSTALLATION.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/INSTALLATION.md b/docs/INSTALLATION.md index 46f2a1be..5bb61ac3 100644 --- a/docs/INSTALLATION.md +++ b/docs/INSTALLATION.md @@ -102,7 +102,7 @@ source ~/.bashrc The generated `~/.local/bin/nvd` launcher uses the repository selected during setup. To run another installed version, set `NVD_REPO` to that checkout root; the launcher then uses its CLI, Pixi environment, and Nextflow sources together: ```bash -NVD_REPO=/home/you/.nvd/v3.2.1 nvd run --samplesheet samples.csv +NVD_REPO=/home/you/.nvd/v3.4.0 nvd run --samplesheet samples.csv ``` `NVD_PIPELINE_ROOT` remains a deprecated compatibility alias during NVD v3. If both variables are set, they must resolve to the same checkout. From d46a0d78b741189b9d7a7db87eaf748630aaa1a8 Mon Sep 17 00:00:00 2001 From: Nick Minor Date: Thu, 30 Jul 2026 18:43:52 -0500 Subject: [PATCH 03/10] build: manage development tools with Mise --- .envrc | 2 - .gitignore | 4 +- Containerfile | 2 +- docs/contributor_guide.md | 24 +++---- flake.nix | 90 -------------------------- mise.lock | 132 ++++++++++++++++++++++++++++++++++++++ mise.toml | 12 ++++ pyproject.toml | 20 +++--- 8 files changed, 166 insertions(+), 120 deletions(-) delete mode 100644 .envrc delete mode 100644 flake.nix create mode 100644 mise.lock create mode 100644 mise.toml diff --git a/.envrc b/.envrc deleted file mode 100644 index bc62fdc5..00000000 --- a/.envrc +++ /dev/null @@ -1,2 +0,0 @@ -use flake . --quiet - diff --git a/.gitignore b/.gitignore index 84214b51..1c08835f 100644 --- a/.gitignore +++ b/.gitignore @@ -15,8 +15,8 @@ !pixi.lock !pyproject.toml !uv.lock -!.envrc -!flake.nix +!mise.lock +!mise.toml # workflows !/workflows diff --git a/Containerfile b/Containerfile index ac5a6fda..e071157e 100644 --- a/Containerfile +++ b/Containerfile @@ -40,7 +40,7 @@ ENV PATH="${HOME}/.pixi/bin:${PATH}" # environment, so skip those direct packages and their dependency subtrees. RUN cd $HOME && \ pixi install --frozen \ - --skip-with-deps apptainer \ + --skip apptainer \ --skip-with-deps rust-script \ --skip-with-deps rust \ --skip-with-deps compilers \ diff --git a/docs/contributor_guide.md b/docs/contributor_guide.md index 39529c76..f9ac4508 100644 --- a/docs/contributor_guide.md +++ b/docs/contributor_guide.md @@ -4,33 +4,27 @@ Welcome to the NVD development guide. This document is for people changing the p ## Quick start -For full-stack development, use Pixi. It installs the Python package, Nextflow, and the bioinformatics tools used by the workflow. +For full-stack development, use Mise to install the repository's pinned Pixi, uv, and just versions. Install Mise and [activate it for your shell](https://mise.jdx.dev/getting-started.html#activate-mise) first. Pixi then installs the Python package, Nextflow, and the bioinformatics tools used by the workflow. ```bash -curl -fsSL https://pixi.sh/install.sh | bash git clone https://github.com/dholab/nvd.git cd nvd -pixi install --frozen +mise trust +mise install +just setup pixi shell ``` -For Python-only work, `uv` is faster and usually enough for CLI, params, preset, setup, taxonomy, and helper-script changes. - -```bash -curl -LsSf https://astral.sh/uv/install.sh | sh -git clone https://github.com/dholab/nvd.git -cd nvd -uv sync -``` +Mise adds `.pixi/envs/default/bin` to `PATH` when its shell activation is enabled. It does not install the Pixi environment merely because you enter the repository; run `just setup` after cloning or when the lockfiles change. -The Nix/direnv setup remains available for contributors who want a more portable system shell: +For Python-only work, the Mise-provided `uv` is faster and usually enough for CLI, params, preset, setup, taxonomy, and helper-script changes. ```bash -curl -L https://nixos.org/nix/install | sh -curl -sfL https://direnv.net/install.sh | bash git clone https://github.com/dholab/nvd.git cd nvd -direnv allow +mise trust +mise install +uv sync ``` ## Current usage model diff --git a/flake.nix b/flake.nix deleted file mode 100644 index b870edf7..00000000 --- a/flake.nix +++ /dev/null @@ -1,90 +0,0 @@ -{ - description = "Reproducible dev shell for the `NVD` bioinformatic processing pipeline"; - - inputs = { - nixpkgs.url = "github:NixOS/nixpkgs/nixpkgs-unstable"; - flake-utils.url = "github:numtide/flake-utils"; - }; - - outputs = - { - self, - nixpkgs, - flake-utils, - ... - }: - flake-utils.lib.eachDefaultSystem ( - system: - let - pkgs = import nixpkgs { - inherit system; - }; - - pixiVersion = "0.74.0"; - pixiRelease = { - aarch64-darwin = { - target = "aarch64-apple-darwin"; - hash = "sha256-t8kqwVMXHSXEanOHJyvaGs5tn530b9NLPJc3MAsyJBU="; - }; - x86_64-darwin = { - target = "x86_64-apple-darwin"; - hash = "sha256-II5HVbzfrIqWxTAU/WOe/nUy5qswGwNgq03vK+Max98="; - }; - aarch64-linux = { - target = "aarch64-unknown-linux-musl"; - hash = "sha256-hJ3dnaP82/yZpZvi8zlzEX0ZrVOjQIYewrmahqm2F98="; - }; - x86_64-linux = { - target = "x86_64-unknown-linux-musl"; - hash = "sha256-BuMYXJdAr5/9NFYQHXvw6uw7KUrxgIkijj8lnDneC2Q="; - }; - }.${system}; - pixi = pkgs.stdenvNoCC.mkDerivation { - pname = "pixi"; - version = pixiVersion; - src = pkgs.fetchurl { - url = "https://github.com/prefix-dev/pixi/releases/download/v${pixiVersion}/pixi-${pixiRelease.target}"; - hash = pixiRelease.hash; - }; - dontUnpack = true; - installPhase = '' - runHook preInstall - install -Dm755 "$src" "$out/bin/pixi" - runHook postInstall - ''; - }; - - in - { - devShells.default = pkgs.mkShell { - name = "NVD"; - - buildInputs = [ - pkgs.stdenv - pkgs.gcc - pkgs.curl - pkgs.wget - pkgs.openjdk - pkgs.git - pkgs.cmake - pkgs.libxml2 - pkgs.libxslt - pkgs.libffi - pixi - pkgs.graphviz - ]; - - shellHook = '' - echo "🔧 Entering NVD dev shell" - export PS1="(NVD) $PS1" - if [ ! -d .pixi/envs/default ]; then - echo "Pixi env not found. Running install..." - pixi install --frozen - fi - - export PATH="$PWD/.pixi/envs/default/bin:$PATH" - ''; - }; - } - ); -} diff --git a/mise.lock b/mise.lock new file mode 100644 index 00000000..e85f9f1e --- /dev/null +++ b/mise.lock @@ -0,0 +1,132 @@ +# @generated - this file is auto-generated by `mise lock` https://mise.jdx.dev/dev-tools/mise-lock.html + +[[tools."github:prefix-dev/pixi"]] +version = "0.74.0" +backend = "github:prefix-dev/pixi" + +[tools."github:prefix-dev/pixi"."platforms.linux-arm64"] +checksum = "sha256:86b9eaace45f28e2fa599e5cb1ca7b123cb3db6c7a9175cda6cb139815724d75" +url = "https://github.com/prefix-dev/pixi/releases/download/v0.74.0/pixi-aarch64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/prefix-dev/pixi/releases/assets/491454287" +provenance = "github-attestations" + +[tools."github:prefix-dev/pixi"."platforms.linux-arm64-musl"] +checksum = "sha256:86b9eaace45f28e2fa599e5cb1ca7b123cb3db6c7a9175cda6cb139815724d75" +url = "https://github.com/prefix-dev/pixi/releases/download/v0.74.0/pixi-aarch64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/prefix-dev/pixi/releases/assets/491454287" +provenance = "github-attestations" + +[tools."github:prefix-dev/pixi"."platforms.linux-x64"] +checksum = "sha256:5729e99fb569a25f0b65dd6f5d1f32fcf3b967b36c1e89b9bab6a3bf845cc655" +url = "https://github.com/prefix-dev/pixi/releases/download/v0.74.0/pixi-x86_64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/prefix-dev/pixi/releases/assets/491454402" +provenance = "github-attestations" + +[tools."github:prefix-dev/pixi"."platforms.linux-x64-musl"] +checksum = "sha256:5729e99fb569a25f0b65dd6f5d1f32fcf3b967b36c1e89b9bab6a3bf845cc655" +url = "https://github.com/prefix-dev/pixi/releases/download/v0.74.0/pixi-x86_64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/prefix-dev/pixi/releases/assets/491454402" +provenance = "github-attestations" + +[tools."github:prefix-dev/pixi"."platforms.macos-arm64"] +checksum = "sha256:abe8daa74933609e65f5cb7a2ae99a1fe292289854d731d1d46b44606fec2acc" +url = "https://github.com/prefix-dev/pixi/releases/download/v0.74.0/pixi-aarch64-apple-darwin.tar.gz" +url_api = "https://api.github.com/repos/prefix-dev/pixi/releases/assets/491454236" +provenance = "github-attestations" + +[tools."github:prefix-dev/pixi"."platforms.macos-x64"] +checksum = "sha256:42565c447d77655960946749ce1d1a5d463db3fceba216ed24d6ced3a7a33b8a" +url = "https://github.com/prefix-dev/pixi/releases/download/v0.74.0/pixi-x86_64-apple-darwin.tar.gz" +url_api = "https://api.github.com/repos/prefix-dev/pixi/releases/assets/491454345" +provenance = "github-attestations" + +[tools."github:prefix-dev/pixi"."platforms.windows-x64"] +checksum = "sha256:5770edf521b1d2be114d93121cb1fa366e8a201ba01c188179d71632804480cf" +url = "https://github.com/prefix-dev/pixi/releases/download/v0.74.0/pixi-x86_64-pc-windows-msvc.zip" +url_api = "https://api.github.com/repos/prefix-dev/pixi/releases/assets/491454373" +provenance = "github-attestations" + +[[tools.just]] +version = "1.57.0" +backend = "aqua:casey/just" + +[tools.just."platforms.linux-arm64"] +checksum = "sha256:f225044a81adea6e0b3a8b9370aaf374e6af76c8735ae263ac993df55fd137ec" +url = "https://github.com/casey/just/releases/download/1.57.0/just-1.57.0-aarch64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/casey/just/releases/assets/481976912" + +[tools.just."platforms.linux-arm64-musl"] +checksum = "sha256:f225044a81adea6e0b3a8b9370aaf374e6af76c8735ae263ac993df55fd137ec" +url = "https://github.com/casey/just/releases/download/1.57.0/just-1.57.0-aarch64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/casey/just/releases/assets/481976912" + +[tools.just."platforms.linux-x64"] +checksum = "sha256:45b548094283cb9739af8f13273b8cddeee869f5b4ef2bb631b1f311cb566155" +url = "https://github.com/casey/just/releases/download/1.57.0/just-1.57.0-x86_64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/casey/just/releases/assets/481976790" + +[tools.just."platforms.linux-x64-musl"] +checksum = "sha256:45b548094283cb9739af8f13273b8cddeee869f5b4ef2bb631b1f311cb566155" +url = "https://github.com/casey/just/releases/download/1.57.0/just-1.57.0-x86_64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/casey/just/releases/assets/481976790" + +[tools.just."platforms.macos-arm64"] +checksum = "sha256:0381db216c2f97ce31d838a1562c1064dfbfa73f5a8a81581338a2cd9217df47" +url = "https://github.com/casey/just/releases/download/1.57.0/just-1.57.0-aarch64-apple-darwin.tar.gz" +url_api = "https://api.github.com/repos/casey/just/releases/assets/481976803" + +[tools.just."platforms.macos-x64"] +checksum = "sha256:5e6ade3698095576274b2b32cc9e5d467185e8e40b04949004c04cc3d7e962dc" +url = "https://github.com/casey/just/releases/download/1.57.0/just-1.57.0-x86_64-apple-darwin.tar.gz" +url_api = "https://api.github.com/repos/casey/just/releases/assets/481976824" + +[tools.just."platforms.windows-x64"] +checksum = "sha256:4c7391d17cb1d17b758b52004ee6411372b8a13ff37c3c9b9031625cb6026e09" +url = "https://github.com/casey/just/releases/download/1.57.0/just-1.57.0-x86_64-pc-windows-msvc.zip" +url_api = "https://api.github.com/repos/casey/just/releases/assets/481977734" + +[[tools.uv]] +version = "0.9.18" +backend = "aqua:astral-sh/uv" + +[tools.uv."platforms.linux-arm64"] +checksum = "sha256:b710ceb9889276cbd7ce04e2ca06b5bd3e288da465bd38f7dd17955c4e703a65" +url = "https://github.com/astral-sh/uv/releases/download/0.9.18/uv-aarch64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/astral-sh/uv/releases/assets/329397626" +provenance = "github-attestations" + +[tools.uv."platforms.linux-arm64-musl"] +checksum = "sha256:b710ceb9889276cbd7ce04e2ca06b5bd3e288da465bd38f7dd17955c4e703a65" +url = "https://github.com/astral-sh/uv/releases/download/0.9.18/uv-aarch64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/astral-sh/uv/releases/assets/329397626" +provenance = "github-attestations" + +[tools.uv."platforms.linux-x64"] +checksum = "sha256:a55ae2d0d53c8f6541bb4d6afc95857ff33a97de8f1d23e9d09acdcb865c4a00" +url = "https://github.com/astral-sh/uv/releases/download/0.9.18/uv-x86_64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/astral-sh/uv/releases/assets/329397666" +provenance = "github-attestations" + +[tools.uv."platforms.linux-x64-musl"] +checksum = "sha256:a55ae2d0d53c8f6541bb4d6afc95857ff33a97de8f1d23e9d09acdcb865c4a00" +url = "https://github.com/astral-sh/uv/releases/download/0.9.18/uv-x86_64-unknown-linux-musl.tar.gz" +url_api = "https://api.github.com/repos/astral-sh/uv/releases/assets/329397666" +provenance = "github-attestations" + +[tools.uv."platforms.macos-arm64"] +checksum = "sha256:dc3bee4abbb3bac267a3985a23ea7617d19d41ff381dbaf560ba415ad65af68f" +url = "https://github.com/astral-sh/uv/releases/download/0.9.18/uv-aarch64-apple-darwin.tar.gz" +url_api = "https://api.github.com/repos/astral-sh/uv/releases/assets/329397618" +provenance = "github-attestations" + +[tools.uv."platforms.macos-x64"] +checksum = "sha256:f86836c637333c65bbc7902acc9c49888eef9fbd15dccbc1946b10e30b041073" +url = "https://github.com/astral-sh/uv/releases/download/0.9.18/uv-x86_64-apple-darwin.tar.gz" +url_api = "https://api.github.com/repos/astral-sh/uv/releases/assets/329397657" +provenance = "github-attestations" + +[tools.uv."platforms.windows-x64"] +checksum = "sha256:28cbe5d30907a774bfe27a517a39b494ec6f7d3816bda8bbf6f9645490449182" +url = "https://github.com/astral-sh/uv/releases/download/0.9.18/uv-x86_64-pc-windows-msvc.zip" +url_api = "https://api.github.com/repos/astral-sh/uv/releases/assets/329397660" +provenance = "github-attestations" diff --git a/mise.toml b/mise.toml new file mode 100644 index 00000000..64b5be41 --- /dev/null +++ b/mise.toml @@ -0,0 +1,12 @@ +min_version = "2026.7.13" + +[tools] +"github:prefix-dev/pixi" = "0.74.0" +just = "1.57.0" +uv = "0.9.18" + +[env] +_.path = [".pixi/envs/default/bin", "bin"] + +[settings] +lockfile = true diff --git a/pyproject.toml b/pyproject.toml index da2d43f7..5c5eec88 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -37,38 +37,38 @@ platforms = ["linux-64", "linux-aarch64", "osx-64", "osx-arm64"] [tool.pixi.dependencies] # Bioinformatics ToolKits # ----------------------------------------------------------------------------------- # +alignoth = ">=1.6.0,<2" bbmap = ">=39.19,<40" blast = ">=2.16.0,<3" deacon = ">=0.15.0,<0.16" fastp = ">=1.0.1,<2" +fastqc = "==0.12.1" flye = ">=2.9.6,<3" metamdbg = ">=1.2,<2" minimap2 = ">=2.22,<3" mmseqs2 = ">=18.8cc5c,<19" +multiqc = "==1.35" myloasm = ">=0.5.1,<0.7" +rust = ">=1.94.0,<1.95" samtools = ">=1.18,<2" seqkit = ">=2.10.0,<3" spades = ">=4.1.0,<5" -rust = ">=1.94.0,<1.95" +sracha = { version = ">=0.3.10,<0.4", channel = "bioconda" } +sourmash = ">=4.9.4,<5" +sourmash_plugin_branchwater = ">=0.9.15,<0.10" +taxburst = ">=0.3.2,<0.4" +vsearch = ">=2.31.0,<3" # System Dependencies # ----------------------------------------------------------------------------------- # compilers = ">=1.9.0,<2" +nextflow = ">=25.10.2,<26" openssl = ">=3.5.0,<4" pkg-config = ">=0.29.2,<0.30" sqlite = ">=3.49.1,<4" unzip = ">=6.0,<7" yaml = ">=0.2.5,<0.3" zstd = ">=1.5.7,<2" -nextflow = ">=25.10.2,<26" -sracha = { version = ">=0.3.10,<0.4", channel = "bioconda" } -sourmash = ">=4.9.4,<5" -sourmash_plugin_branchwater = ">=0.9.15,<0.10" -taxburst = ">=0.3.2,<0.4" -vsearch = ">=2.31.0,<3" -alignoth = ">=1.6.0,<2" -fastqc = "==0.12.1" -multiqc = "==1.35" # ----------------------------------------------------------------------------------- # From 2baa60e4e582b5f7cb5e28515ca1f6c5b2fc51cf Mon Sep 17 00:00:00 2001 From: Nick Minor Date: Thu, 30 Jul 2026 18:50:24 -0500 Subject: [PATCH 04/10] feat: stream SRA reads through sracha --- README.md | 2 + conf/results.config | 2 +- docs/nvd_cli_guide.md | 2 + modules/deacon.nf | 101 +++++++++++++++++++++++++++++++ pixi.lock | 64 ++++++++++---------- pyproject.toml | 2 +- subworkflows/gather_reads.nf | 32 ++-------- subworkflows/preprocess_reads.nf | 24 +++++++- tests/test_read_filtering.py | 23 ++++--- tests/test_target_enrichment.py | 24 +++++--- workflows/nvd_main.nf | 5 +- 11 files changed, 197 insertions(+), 84 deletions(-) diff --git a/README.md b/README.md index 6da44fd5..55aa6592 100644 --- a/README.md +++ b/README.md @@ -212,6 +212,8 @@ Then generate a samplesheet from that accession list: nvd samplesheet generate --from-sra accessions.txt --platform illumina --output samplesheet.csv ``` +NVD decodes each SRA run as a stream through target enrichment rather than materializing decoded raw FASTQ files. Raw-read FastQC therefore runs only for local FASTQ inputs; the `skip_fastqc` setting controls those local tasks. SRA reads rejoin local reads immediately after target enrichment and receive the same subsequent preprocessing. + If you want to inspect what NVD would write before touching the filesystem, use dry-run mode: ```bash diff --git a/conf/results.config b/conf/results.config index 7f8bc390..00b3306b 100644 --- a/conf/results.config +++ b/conf/results.config @@ -109,7 +109,7 @@ process { // ========================================================================= // 01 — Target-enriched reads // ========================================================================= - withName: 'DEACON_ENRICH_TARGET_READS' { + withName: '^(DEACON_ENRICH_TARGET_READS|DEACON_ENRICH_SRA_READS)$' { publishDir = [ [ path: { params.target_enriched_reads }, diff --git a/docs/nvd_cli_guide.md b/docs/nvd_cli_guide.md index 0ec599db..d5bf0db4 100644 --- a/docs/nvd_cli_guide.md +++ b/docs/nvd_cli_guide.md @@ -187,6 +187,8 @@ Generate from SRA accessions, one accession per line: nvd samplesheet generate --from-sra accessions.txt --platform illumina --output samplesheet.csv ``` +NVD decodes each SRA run as a stream through target enrichment rather than materializing decoded raw FASTQ files. Raw-read FastQC therefore runs only for local FASTQ inputs; the `skip_fastqc` setting controls those local tasks. SRA reads rejoin local reads immediately after target enrichment and receive the same subsequent preprocessing. + Preview without writing: ```bash diff --git a/modules/deacon.nf b/modules/deacon.nf index 4a7a9243..ecb111ec 100644 --- a/modules/deacon.nf +++ b/modules/deacon.nf @@ -236,6 +236,107 @@ process DEACON_ENRICH_TARGET_READS { """ } +process DEACON_ENRICH_SRA_READS { + /* Stream one resolved SRA run through deacon without materializing decoded FASTQ files. */ + + tag "${id}, ${run_accession}" + label "medium" + + errorStrategy { task.attempt < 3 ? 'retry' : 'ignore' } + maxRetries 2 + maxForks params.max_concurrent_downloads + + input: + tuple val(id), val(platform), val(run_accession), path(deacon_idx), val(target_enrichment_enabled) + + output: + tuple val(id), val(platform), path("${id}.sra_read_structure.txt"), path("${id}.target_enriched.fastq.gz"), emit: reads + tuple val(id), path("${id}.deacon_filter.json"), emit: stats + + script: + def cpus = task.cpus as int + def sracha_threads = Math.max(1, cpus.intdiv(2)) + def deacon_threads = Math.max(1, cpus - sracha_threads) + def deplete_arg = target_enrichment_enabled ? "" : "--deplete" + """ + set -euo pipefail + + metadata_file='${id}.sracha_info.tsv' + sracha info --format tsv '${run_accession}' > "\${metadata_file}" + + metadata_lines=() + while IFS= read -r line || [[ -n "\${line}" ]]; do + metadata_lines+=("\${line}") + done < "\${metadata_file}" + + if [[ \${#metadata_lines[@]} -ne 2 ]]; then + printf 'Expected one sracha metadata row for %s, received %s lines\n' \ + '${run_accession}' "\${#metadata_lines[@]}" >&2 + exit 1 + fi + + expected_header=\$'accession\tarchive_type\tlayout\tnreads\tspots\tsize_bytes\tplatform\tmd5' + if [[ "\${metadata_lines[0]}" != "\${expected_header}" ]]; then + printf 'Unexpected sracha metadata header for %s: %s\n' \ + '${run_accession}' "\${metadata_lines[0]}" >&2 + exit 1 + fi + + IFS=\$'\t' read -r -a metadata_fields <<< "\${metadata_lines[1]}" + if [[ \${#metadata_fields[@]} -ne 8 ]]; then + printf 'Expected eight sracha metadata fields for %s, received %s\n' \ + '${run_accession}' "\${#metadata_fields[@]}" >&2 + exit 1 + fi + + resolved_accession="\${metadata_fields[0]}" + layout="\${metadata_fields[2]}" + nreads="\${metadata_fields[3]}" + + if [[ "\${resolved_accession}" != '${run_accession}' ]]; then + printf 'Sracha resolved %s while %s was requested\n' \ + "\${resolved_accession}" '${run_accession}' >&2 + exit 1 + fi + + case "\${layout}/\${nreads}" in + SINGLE/1) + read_structure='single' + deacon_inputs=(-) + ;; + PAIRED/2) + read_structure='interleaved' + deacon_inputs=(- -) + ;; + *) + printf 'Unsupported sracha read layout for %s: layout=%s nreads=%s\n' \ + '${run_accession}' "\${layout}" "\${nreads}" >&2 + exit 1 + ;; + esac + + printf '%s\n' "\${read_structure}" > '${id}.sra_read_structure.txt' + + sracha get \ + --output-dir . \ + --stdout \ + --split interleaved \ + --threads ${sracha_threads} \ + --no-progress \ + --yes \ + '${run_accession}' \ + | deacon filter \ + ${deplete_arg} \ + --threads ${deacon_threads} \ + --abs-threshold ${params.virus_abs_threshold} \ + --rel-threshold ${params.virus_rel_threshold} \ + --summary '${id}.deacon_filter.json' \ + --output '${id}.target_enriched.fastq.gz' \ + ${deacon_idx} \ + "\${deacon_inputs[@]}" + """ +} + process DEACON_FILTER_CONTIGS { /* * Retain target contigs using deacon filter on assembled FASTA. diff --git a/pixi.lock b/pixi.lock index ab717ea0..49c3226c 100644 --- a/pixi.lock +++ b/pixi.lock @@ -52,7 +52,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/linux-64/samtools-1.24-h9dcdb79_1.conda - conda: https://conda.anaconda.org/bioconda/linux-64/seqkit-2.13.0-he881be0_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/spades-4.3.0-hde4eca7_1.conda - - conda: https://conda.anaconda.org/bioconda/linux-64/sracha-0.3.11-h54198d6_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-64/sracha-0.5.0-h54198d6_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/vsearch-2.31.0-hd2be7a0_0.conda - conda: https://conda.anaconda.org/bioconda/noarch/fastqc-0.12.1-hdfd78af_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/noarch/multiqc-1.35-pyhdfd78af_1.conda @@ -366,7 +366,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/linux-aarch64/samtools-1.24-h391949c_1.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/seqkit-2.13.0-h8865c2f_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/spades-4.3.0-h07e4b28_1.conda - - conda: https://conda.anaconda.org/bioconda/linux-aarch64/sracha-0.3.11-h389c417_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-aarch64/sracha-0.5.0-h389c417_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/vsearch-2.31.0-hd0c913a_0.conda - conda: https://conda.anaconda.org/bioconda/noarch/fastqc-0.12.1-hdfd78af_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/noarch/multiqc-1.35-pyhdfd78af_1.conda @@ -686,7 +686,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-64/samtools-1.24-he9a8aba_1.conda - conda: https://conda.anaconda.org/bioconda/osx-64/seqkit-2.13.0-h193322a_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/spades-4.3.0-hb11480c_1.conda - - conda: https://conda.anaconda.org/bioconda/osx-64/sracha-0.3.11-h518e4cb_0.conda + - conda: https://conda.anaconda.org/bioconda/osx-64/sracha-0.5.0-h518e4cb_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/vsearch-2.31.0-h7aab66d_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/annotated-types-0.8.0-pyhd8ed1ab_0.conda @@ -973,7 +973,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-arm64/samtools-1.24-h36b3a25_1.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/seqkit-2.13.0-hd5f1084_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/spades-4.3.0-hd468e49_1.conda - - conda: https://conda.anaconda.org/bioconda/osx-arm64/sracha-0.3.11-h1a4e967_0.conda + - conda: https://conda.anaconda.org/bioconda/osx-arm64/sracha-0.5.0-h1a4e967_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/vsearch-2.31.0-h0448ff9_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/annotated-types-0.8.0-pyhd8ed1ab_0.conda @@ -1258,7 +1258,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/linux-64/samtools-1.24-h9dcdb79_1.conda - conda: https://conda.anaconda.org/bioconda/linux-64/seqkit-2.13.0-he881be0_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/spades-4.3.0-hde4eca7_1.conda - - conda: https://conda.anaconda.org/bioconda/linux-64/sracha-0.3.11-h54198d6_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-64/sracha-0.5.0-h54198d6_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/vsearch-2.31.0-hd2be7a0_0.conda - conda: https://conda.anaconda.org/bioconda/noarch/fastqc-0.12.1-hdfd78af_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/noarch/multiqc-1.35-pyhdfd78af_1.conda @@ -1670,7 +1670,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/linux-aarch64/samtools-1.24-h391949c_1.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/seqkit-2.13.0-h8865c2f_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/spades-4.3.0-h07e4b28_1.conda - - conda: https://conda.anaconda.org/bioconda/linux-aarch64/sracha-0.3.11-h389c417_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-aarch64/sracha-0.5.0-h389c417_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/vsearch-2.31.0-hd0c913a_0.conda - conda: https://conda.anaconda.org/bioconda/noarch/fastqc-0.12.1-hdfd78af_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/noarch/multiqc-1.35-pyhdfd78af_1.conda @@ -2088,7 +2088,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-64/samtools-1.24-he9a8aba_1.conda - conda: https://conda.anaconda.org/bioconda/osx-64/seqkit-2.13.0-h193322a_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/spades-4.3.0-hb11480c_1.conda - - conda: https://conda.anaconda.org/bioconda/osx-64/sracha-0.3.11-h518e4cb_0.conda + - conda: https://conda.anaconda.org/bioconda/osx-64/sracha-0.5.0-h518e4cb_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/vsearch-2.31.0-h7aab66d_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/annotated-types-0.8.0-pyhd8ed1ab_0.conda @@ -2474,7 +2474,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-arm64/samtools-1.24-h36b3a25_1.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/seqkit-2.13.0-hd5f1084_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/spades-4.3.0-hd468e49_1.conda - - conda: https://conda.anaconda.org/bioconda/osx-arm64/sracha-0.3.11-h1a4e967_0.conda + - conda: https://conda.anaconda.org/bioconda/osx-arm64/sracha-0.5.0-h1a4e967_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/vsearch-2.31.0-h0448ff9_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_3.conda - conda: https://conda.anaconda.org/conda-forge/noarch/annotated-types-0.8.0-pyhd8ed1ab_0.conda @@ -3163,9 +3163,9 @@ packages: - spades >=4.3.0,<5.0a0 size: 18750876 timestamp: 1785589300167 -- conda: https://conda.anaconda.org/bioconda/linux-64/sracha-0.3.11-h54198d6_0.conda - sha256: 3d2afd00513bb756cf7ae10774125f2e5ccdecc9eb16e38b0eb8eb28b4b42beb - md5: 3c478007d246ca08bbe3dbc8bb082b31 +- conda: https://conda.anaconda.org/bioconda/linux-64/sracha-0.5.0-h54198d6_0.conda + sha256: caabfd53f72dd500a866279703dbf724819d6a667f46b59f8ed642ff10842175 + md5: 23417d39f3701e397ff7bcfd93227134 depends: - libgcc >=14 constrains: @@ -3174,9 +3174,9 @@ packages: purls: [] run_exports: weak: - - sracha >=0.3.11,<0.4.0a0 - size: 3025009 - timestamp: 1783873616890 + - sracha >=0.5.0,<0.6.0a0 + size: 3061014 + timestamp: 1785719882966 - conda: https://conda.anaconda.org/bioconda/linux-64/vsearch-2.31.0-hd2be7a0_0.conda sha256: 2cb6740460d605597869ba3223b4ea337f19256d5fc9de6e14510f3f6171ac0b md5: 1688be9a6aa97ea509e853055f57b03d @@ -3537,9 +3537,9 @@ packages: - spades >=4.3.0,<5.0a0 size: 34686930 timestamp: 1785588815442 -- conda: https://conda.anaconda.org/bioconda/linux-aarch64/sracha-0.3.11-h389c417_0.conda - sha256: 883b7a796f9978075c97af1949d5f83728690e6f0edaca2d26b4778ce570f70d - md5: e025f830c42305327f4380f32d5d5480 +- conda: https://conda.anaconda.org/bioconda/linux-aarch64/sracha-0.5.0-h389c417_0.conda + sha256: 582d2a97d3977f64cf499272cbcda26d83388478bd2626b1ac3c2d43fcd3d846 + md5: 0179cdf699064f0dde1f9e08a45627e6 depends: - libgcc >=14 constrains: @@ -3548,9 +3548,9 @@ packages: purls: [] run_exports: weak: - - sracha >=0.3.11,<0.4.0a0 - size: 2947778 - timestamp: 1783873416670 + - sracha >=0.5.0,<0.6.0a0 + size: 2987330 + timestamp: 1785719661546 - conda: https://conda.anaconda.org/bioconda/linux-aarch64/vsearch-2.31.0-hd0c913a_0.conda sha256: e70bab441117c0420a4d77802ccf9123bbec3a76070956295e47688967197b39 md5: 549678a7db61a3883543250b65d7c1a2 @@ -4059,18 +4059,18 @@ packages: - spades >=4.3.0,<5.0a0 size: 26944692 timestamp: 1785592045332 -- conda: https://conda.anaconda.org/bioconda/osx-64/sracha-0.3.11-h518e4cb_0.conda - sha256: a391edf05d0f74060a22a837c2c8a3db3ce5ba01b3d4d28a25c38a645d7a1c87 - md5: 1e8d9b1483f9fee8e9a8bc49ed3f8570 +- conda: https://conda.anaconda.org/bioconda/osx-64/sracha-0.5.0-h518e4cb_0.conda + sha256: 0a2f298751071c80d99bc506b043f169ae4fc3586d30a785442220d3ac144a48 + md5: 43e83d4f2b6e8a1da473f715361dac44 constrains: - __osx >=11.0 license: MIT purls: [] run_exports: weak: - - sracha >=0.3.11,<0.4.0a0 - size: 2937957 - timestamp: 1783875310358 + - sracha >=0.5.0,<0.6.0a0 + size: 2972246 + timestamp: 1785720990903 - conda: https://conda.anaconda.org/bioconda/osx-64/vsearch-2.31.0-h7aab66d_0.conda sha256: 47236db3b8c996a26b3fbf6bfade7c209e3c222a799e53e5c8bb414a755a754f md5: cada888bf20c710bf62e724b6f01dcd7 @@ -4434,18 +4434,18 @@ packages: - spades >=4.3.0,<5.0a0 size: 22442465 timestamp: 1785588192881 -- conda: https://conda.anaconda.org/bioconda/osx-arm64/sracha-0.3.11-h1a4e967_0.conda - sha256: 3f89f535599b7db8e84f08a33026160cb9d6949fb980c5a0c12d07e4fa6f7381 - md5: 1d4da04d4876deb89f431c89163b4b70 +- conda: https://conda.anaconda.org/bioconda/osx-arm64/sracha-0.5.0-h1a4e967_0.conda + sha256: af0bb6a50bf6940ead061f998992598982c2317c93aef9d00a40a5a86d457b15 + md5: 3c48083993c15db8ed294596c2ca27a2 constrains: - __osx >=11.0 license: MIT purls: [] run_exports: weak: - - sracha >=0.3.11,<0.4.0a0 - size: 2649966 - timestamp: 1783873350576 + - sracha >=0.5.0,<0.6.0a0 + size: 2686424 + timestamp: 1785719457985 - conda: https://conda.anaconda.org/bioconda/osx-arm64/vsearch-2.31.0-h0448ff9_0.conda sha256: 992bb9bf6858d0c53f502f11575fe45526862b47277c5d2bbd3de008136af7f0 md5: 003a2cc62a2f4350d19d549a6db39e25 diff --git a/pyproject.toml b/pyproject.toml index 5c5eec88..a2e871f8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -53,7 +53,7 @@ rust = ">=1.94.0,<1.95" samtools = ">=1.18,<2" seqkit = ">=2.10.0,<3" spades = ">=4.1.0,<5" -sracha = { version = ">=0.3.10,<0.4", channel = "bioconda" } +sracha = { version = ">=0.5.0", channel = "bioconda" } sourmash = ">=4.9.4,<5" sourmash_plugin_branchwater = ">=0.9.15,<0.10" taxburst = ">=0.3.2,<0.4" diff --git a/subworkflows/gather_reads.nf b/subworkflows/gather_reads.nf index 61a4f76f..6499e81f 100644 --- a/subworkflows/gather_reads.nf +++ b/subworkflows/gather_reads.nf @@ -1,5 +1,3 @@ -include { FETCH_FASTQ } from "../modules/sratools" - process RESOLVE_READ_INPUTS { label "low" @@ -35,7 +33,10 @@ workflow GATHER_READS { ch_sra_accessions = ch_resolved_reads .filter { rec -> rec.source == "sra" } - .map { rec -> tuple(rec.sample_id, rec.platform, rec.srr) } + .map { rec -> + def accession = rec.srr.toString().toUpperCase(java.util.Locale.ROOT) + tuple(rec.sample_id, rec.platform, accession) + } ch_local_bundles = ch_resolved_reads .filter { rec -> rec.source != "sra" } @@ -58,30 +59,9 @@ workflow GATHER_READS { tuple(meta, r1 + r2) } - FETCH_FASTQ(ch_sra_accessions) - - ch_sra_bundles = FETCH_FASTQ.out - .map { sample_id, platform, fastq_files -> - def files = (fastq_files instanceof List ? fastq_files : [fastq_files]) - .collect { file(it) } - .sort { a, b -> a.getName() <=> b.getName() } - def read1 = files.find { path -> path.getName().endsWith("1.fastq") || path.getName().endsWith("1.fastq.gz") || path.getName().contains("_R1_") } - def read2 = files.find { path -> path.getName().endsWith("2.fastq") || path.getName().endsWith("2.fastq.gz") || path.getName().contains("_R2_") } - if (read1 && read2) { - def meta = [id: sample_id, platform: platform, source: "sra", read_mode: "paired", r1_count: 1, deacon_read_structure: "interleaved"] - return tuple(meta, [read1, read2]) - } - if (files.size() == 1) { - def meta = [id: sample_id, platform: platform, source: "sra", read_mode: "single", r1_count: 1, deacon_read_structure: "single"] - return tuple(meta, files) - } - throw new IllegalArgumentException("Could not determine SRA FASTQ pairing for ${sample_id}: ${files*.getName()}") - } - - ch_raw_reads = ch_local_bundles.mix(ch_sra_bundles) - emit: - reads = ch_raw_reads + reads = ch_local_bundles + sra_accessions = ch_sra_accessions resolved_reads = RESOLVE_READ_INPUTS.out.jsonl } diff --git a/subworkflows/preprocess_reads.nf b/subworkflows/preprocess_reads.nf index a5aa5827..0184f8f4 100644 --- a/subworkflows/preprocess_reads.nf +++ b/subworkflows/preprocess_reads.nf @@ -4,6 +4,7 @@ include { DEACON_BUILD_TARGET_INDEX_FROM_FASTA } from "../modules/dea include { DEACON_BUILD_INDEX_FROM_FASTA } from "../modules/deacon" include { DEACON_UNION_INDEXES } from "../modules/deacon" include { DEACON_ENRICH_TARGET_READS } from "../modules/deacon" +include { DEACON_ENRICH_SRA_READS } from "../modules/deacon" include { DEACON_DEPLETE } from "../modules/deacon" include { MERGE_PAIRS ; DEDUP_WITH_CLUMPIFY ; TRIM_ADAPTERS ; FILTER_READS } from "../modules/bbmap" include { PROFILE_FASTX as PROFILE_READS ; PLOT_FASTX_LENGTH_PROFILE as PLOT_READ_LENGTH_PROFILES ; PLOT_FASTX_QUALITY_PROFILE as PLOT_READ_QUALITY_PROFILES } from "../modules/fastx" @@ -12,6 +13,7 @@ include { FASTQC_RAW } from "../modules/fastqc" workflow PREPROCESS_READS { take: ch_read_bundles // tuple(meta, read_files) from GATHER_READS + ch_sra_accessions // tuple(sample_id, platform, accession) from GATHER_READS main: @@ -60,13 +62,29 @@ workflow PREPROCESS_READS { ch_read_bundles.combine(ch_target_index) .combine(ch_target_enrichment_enabled) ) + DEACON_ENRICH_SRA_READS( + ch_sra_accessions.combine(ch_target_index) + .combine(ch_target_enrichment_enabled) + ) // ------------------------------------------------------------------------- // Step 2: Inlined preprocessing on target-enriched reads // ------------------------------------------------------------------------- + ch_sra_target_reads = DEACON_ENRICH_SRA_READS.out.reads + .map { sample_id, platform, read_structure_file, reads -> + def read_structure = read_structure_file.toFile().text.trim() + tuple(sample_id, platform, read_structure, reads) + } + + ch_target_reads = DEACON_ENRICH_TARGET_READS.out.reads + .mix(ch_sra_target_reads) + + ch_target_enrichment_stats = DEACON_ENRICH_TARGET_READS.out.stats + .mix(DEACON_ENRICH_SRA_READS.out.stats) + // Extract input and retained counts from Deacon's summary. A malformed // summary is a process failure, never an empty-sample completion. - ch_enrichment_counts = DEACON_ENRICH_TARGET_READS.out.stats + ch_enrichment_counts = ch_target_enrichment_stats .map { sample_id, json_file -> def summary = new groovy.json.JsonSlurper().parse(json_file.toFile()) def input_count = summary.seqs_in @@ -79,7 +97,7 @@ workflow PREPROCESS_READS { ch_read_counts = ch_enrichment_counts .map { sample_id, input_count, _retained_count -> tuple(sample_id, input_count) } - ch_target_reads_by_retention = DEACON_ENRICH_TARGET_READS.out.reads + ch_target_reads_by_retention = ch_target_reads .join(ch_enrichment_counts, by: 0) .branch { _sample_id, _platform, _read_structure, _reads, _input_count, retained_count -> retained: retained_count > 0 @@ -414,7 +432,7 @@ workflow PREPROCESS_READS { single_reads_for_mapback = ch_single_reads_for_mapback read_counts = ch_read_counts complete_empty_samples = ch_complete_empty_samples - target_enrichment_stats = DEACON_ENRICH_TARGET_READS.out.stats + target_enrichment_stats = ch_target_enrichment_stats depletion_stats = ch_depletion_stats raw_fastqc_packages = FASTQC_RAW.out.packages raw_fastqc_zips = FASTQC_RAW.out.zips diff --git a/tests/test_read_filtering.py b/tests/test_read_filtering.py index a2c2f2e3..33ccd83f 100644 --- a/tests/test_read_filtering.py +++ b/tests/test_read_filtering.py @@ -306,16 +306,19 @@ def test_preprocessing_deduplicates_after_adapter_trimming(tmp_path: Path) -> No include {{ PREPROCESS_READS }} from '{PREPROCESS_READS}' workflow {{ - PREPROCESS_READS(Channel.of(tuple( - [ - id: 'sample_A', - platform: 'illumina', - read_mode: 'single', - r1_count: 1, - deacon_read_structure: 'single', - ], - [file('{reads}')], - ))) + PREPROCESS_READS( + Channel.of(tuple( + [ + id: 'sample_A', + platform: 'illumina', + read_mode: 'single', + r1_count: 1, + deacon_read_structure: 'single', + ], + [file('{reads}')], + )), + Channel.empty(), + ) PREPROCESS_READS.out.read_batches.view {{ _id, _platform, _structure, _query_class, output -> "FINAL_READS:${{output.name}}" diff --git a/tests/test_target_enrichment.py b/tests/test_target_enrichment.py index 69389e8d..3bb005ea 100644 --- a/tests/test_target_enrichment.py +++ b/tests/test_target_enrichment.py @@ -156,6 +156,7 @@ def test_zero_target_reads_stop_before_preprocessing(tmp_path: Path) -> None: params.no_enrichment = false params.skip_fastqc = true +params.max_concurrent_downloads = 1 params.merge_pairs = false params.virus_index = '{target_index}' params.virus_index_url = null @@ -178,16 +179,19 @@ def test_zero_target_reads_stop_before_preprocessing(tmp_path: Path) -> None: include {{ PREPROCESS_READS }} from '{PREPROCESS_READS}' workflow {{ - PREPROCESS_READS(Channel.of(tuple( - [ - id: 'sample_A', - platform: 'illumina', - read_mode: 'single', - r1_count: 1, - deacon_read_structure: 'single', - ], - [file('{reads}')], - ))) + PREPROCESS_READS( + Channel.of(tuple( + [ + id: 'sample_A', + platform: 'illumina', + read_mode: 'single', + r1_count: 1, + deacon_read_structure: 'single', + ], + [file('{reads}')], + )), + Channel.empty(), + ) PREPROCESS_READS.out.reads.view {{ sample_id, _platform, _structure, output -> "READ_OUTPUT:${{sample_id}}:${{output.size()}}" diff --git a/workflows/nvd_main.nf b/workflows/nvd_main.nf index 2627c77e..17201320 100644 --- a/workflows/nvd_main.nf +++ b/workflows/nvd_main.nf @@ -62,7 +62,10 @@ workflow NVD_MAIN { GATHER_READS(ch_samplesheet) - PREPROCESS_READS(GATHER_READS.out.reads) + PREPROCESS_READS( + GATHER_READS.out.reads, + GATHER_READS.out.sra_accessions, + ) ch_risk_group_lookup = Channel.value(file("${projectDir}/assets/human_virus_risk_group_lookup.tsv")) From c30c214576d38fa9f652f790c9783047414dda8a Mon Sep 17 00:00:00 2001 From: Nick Minor Date: Mon, 27 Jul 2026 15:53:35 -0500 Subject: [PATCH 05/10] version 3.5.0 minor bump release candidate Bump __version__, the Nextflow manifest, pyproject, and uv.lock to 3.5.0, and publish schemas/nvd-params.v3.5.0.schema.json as the parameter contract for the release. The rolling nvd-params.latest.schema.json symlink and params.SCHEMA_URL both move to it, so IDE completion and generated templates describe 3.5. The v3.5.0 schema is a copy of the v3.4.0 contract with only $id changed. That contract already reflects the v3.4.0 merge: the six sourmash reference and threshold params are gone, labkey_insert_batch_size and skip_unassembled_read_queries are present, and merge_pairs defaults to true. The rapid-screen skip option that originally introduced this schema was dropped during the rebase onto the v3.4.0 merge, since sourmash rapid screening no longer exists to skip. Only its version plumbing survives here. Also drops bin/filter_non_virus_blast_nodes.py, bin/summarize_blast_filtering.py, and bin/test_filter_non_virus_blast_nodes.py. "feat: retain BLAST hits across taxa" deleted all three, a bad merge resurrected them here (the summarize_blast_filtering.py copy predating its qseqid column check), and nothing in the pipeline referenced them. --- lib/py_nvd/__init__.py | 2 +- lib/py_nvd/params.py | 2 +- lib/py_nvd/test_version_sync.py | 8 +- nextflow.config | 2 +- pyproject.toml | 2 +- schemas/README.md | 3 +- schemas/nvd-params.latest.schema.json | 2 +- schemas/nvd-params.v3.5.0.schema.json | 453 ++++++++++++++++++++++++++ uv.lock | 2 +- 9 files changed, 465 insertions(+), 11 deletions(-) create mode 100644 schemas/nvd-params.v3.5.0.schema.json diff --git a/lib/py_nvd/__init__.py b/lib/py_nvd/__init__.py index 8deedb1f..41304222 100644 --- a/lib/py_nvd/__init__.py +++ b/lib/py_nvd/__init__.py @@ -1,6 +1,6 @@ """NVD CLI and pipeline helper library.""" -__version__ = "3.4.1" +__version__ = "3.5.0" # Re-export key modules for convenient access. from py_nvd import models, params, paths, presets, taxonomy diff --git a/lib/py_nvd/params.py b/lib/py_nvd/params.py index f7719d16..c146739b 100644 --- a/lib/py_nvd/params.py +++ b/lib/py_nvd/params.py @@ -22,7 +22,7 @@ SCHEMA_FILENAME = "nvd-params.latest.schema.json" # GitHub raw URL for schema (fallback and for generated templates) -SCHEMA_URL = "https://raw.githubusercontent.com/dholab/nvd/main/schemas/nvd-params.v3.4.0.schema.json" +SCHEMA_URL = "https://raw.githubusercontent.com/dholab/nvd/main/schemas/nvd-params.v3.5.0.schema.json" def _find_schema_path() -> Path: diff --git a/lib/py_nvd/test_version_sync.py b/lib/py_nvd/test_version_sync.py index 386c0fd4..1024c524 100644 --- a/lib/py_nvd/test_version_sync.py +++ b/lib/py_nvd/test_version_sync.py @@ -112,13 +112,13 @@ def test_preprocess_param_is_gone() -> None: assert "preprocess" not in NvdParams.model_fields -def test_latest_params_schema_points_to_v3_4() -> None: - """The rolling schema link should expose the v3.4 parameter contract.""" +def test_latest_params_schema_points_to_v3_5() -> None: + """The rolling schema link should expose the v3.5 parameter contract.""" latest_schema = ROOT / "schemas" / "nvd-params.latest.schema.json" assert latest_schema.is_symlink() - assert latest_schema.readlink() == Path("nvd-params.v3.4.0.schema.json") - assert SCHEMA_URL.endswith("/nvd-params.v3.4.0.schema.json") + assert latest_schema.readlink() == Path("nvd-params.v3.5.0.schema.json") + assert SCHEMA_URL.endswith("/nvd-params.v3.5.0.schema.json") def test_v3_3_2_schema_corrects_only_the_read_entropy_default() -> None: diff --git a/nextflow.config b/nextflow.config index c33f2634..9219943b 100644 --- a/nextflow.config +++ b/nextflow.config @@ -338,7 +338,7 @@ manifest { homePage = '' mainScript = 'main.nf' defaultBranch = 'main' - version = '3.4.1' + version = '3.5.0' description = '' author = '' } diff --git a/pyproject.toml b/pyproject.toml index a2e871f8..5149193f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ authors = [ { name = "William Gardner", email = "wkgardner@wisc.edu" }, ] requires-python = ">= 3.11, < 3.14" -version = "3.4.1" +version = "3.5.0" dependencies = [ "altair>=6.1.0", "biopython>=1.85", diff --git a/schemas/README.md b/schemas/README.md index 1937b3d1..6218f3a3 100644 --- a/schemas/README.md +++ b/schemas/README.md @@ -12,6 +12,7 @@ This directory contains JSON Schema definitions for the NVD pipeline. | `nvd-params.v3.3.0.schema.json` | Pipeline parameters schema (version 3.3.0) | | `nvd-params.v3.3.2.schema.json` | Pipeline parameters schema (version 3.3.2) | | `nvd-params.v3.4.0.schema.json` | Pipeline parameters schema (version 3.4.0) | +| `nvd-params.v3.5.0.schema.json` | Pipeline parameters schema (version 3.5.0) | | `nvd-params.latest.schema.json` | Symlink to the current version | ## Usage @@ -74,7 +75,7 @@ points to the current version for users who want to track updates automatically. For reproducibility, you can reference a specific version: ```yaml -# yaml-language-server: $schema=https://raw.githubusercontent.com/dholab/nvd/main/schemas/nvd-params.v3.4.0.schema.json +# yaml-language-server: $schema=https://raw.githubusercontent.com/dholab/nvd/main/schemas/nvd-params.v3.5.0.schema.json ``` ## Validation diff --git a/schemas/nvd-params.latest.schema.json b/schemas/nvd-params.latest.schema.json index dedd9b38..ad334926 120000 --- a/schemas/nvd-params.latest.schema.json +++ b/schemas/nvd-params.latest.schema.json @@ -1 +1 @@ -nvd-params.v3.4.0.schema.json \ No newline at end of file +nvd-params.v3.5.0.schema.json \ No newline at end of file diff --git a/schemas/nvd-params.v3.5.0.schema.json b/schemas/nvd-params.v3.5.0.schema.json new file mode 100644 index 00000000..127bea39 --- /dev/null +++ b/schemas/nvd-params.v3.5.0.schema.json @@ -0,0 +1,453 @@ +{ + "$schema": "https://json-schema.org/draft/2020-12/schema", + "$id": "https://raw.githubusercontent.com/dholab/nvd/main/schemas/nvd-params.v3.5.0.schema.json", + "title": "NVD Pipeline Parameters", + "description": "Parameters for the Novel Virus Detection (NVD) metagenomics pipeline. Use with Nextflow's -params-file option or register as a preset with nvd preset register.", + "type": "object", + "properties": { + "labkey": { + "type": "boolean", + "default": false, + "description": "Enable LabKey integration for result uploads. When enabled, all labkey_* params must be set." + }, + "labkey_server": { + "type": "string", + "description": "LabKey server URL" + }, + "labkey_project_name": { + "type": "string", + "description": "LabKey project name" + }, + "labkey_webdav": { + "type": "string", + "description": "LabKey WebDAV endpoint URL" + }, + "labkey_schema": { + "type": "string", + "description": "LabKey schema name" + }, + "labkey_blast_meta_hits_list": { + "type": "string", + "description": "LabKey list name for BLAST meta hits" + }, + "labkey_blast_fasta_list": { + "type": "string", + "description": "LabKey list name for BLAST FASTA uploads" + }, + "labkey_insert_batch_size": { + "type": "integer", + "minimum": 1, + "default": 1000, + "description": "Rows per LabKey insert call. Read-derived query classes produce payloads large enough that a single call can hang the server." + }, + "samplesheet": { + "type": "string", + "format": "path", + "description": "Path to samplesheet CSV with columns: sample_id, srr, platform, fastq1, fastq2" + }, + "results": { + "type": "string", + "format": "path", + "description": "Directory for pipeline output files" + }, + "experiment_id": { + "type": "string", + "description": "Experiment identifier for tracking and LabKey integration" + }, + "max_concurrent_downloads": { + "type": "integer", + "default": 3, + "minimum": 1, + "description": "Maximum number of concurrent SRA downloads" + }, + "cleanup": { + "type": "boolean", + "description": "Whether to empty the work directory after successful completion" + }, + "work_dir": { + "type": [ + "string", + "null" + ], + "format": "path", + "description": "Nextflow work directory for intermediate files" + }, + "experimental": { + "type": "boolean", + "default": false, + "description": "Enable experimental release-candidate features that may change or be removed." + }, + "skip_assembly": { + "type": "boolean", + "default": false, + "description": "Skip SPAdes assembly and all downstream contig classification." + }, + "skip_blast": { + "type": "boolean", + "default": false, + "description": "Skip MEGABLAST and BLASTN contig search." + }, + "skip_fastqc": { + "type": "boolean", + "default": false, + "description": "Skip per-file raw-read FastQC." + }, + "skip_unassembled_read_queries": { + "type": "boolean", + "default": false, + "description": "Skip BLAST querying of unassembled reads (overlap-merged pairs and single reads), leaving assembly contigs as the only query classes." + }, + "blast_db_version": { + "type": "string", + "description": "BLAST database version identifier" + }, + "virus_index_version": { + "type": "string", + "description": "Virus enrichment index version" + }, + "nvd_files": { + "type": "string", + "format": "path", + "description": "Path to NVD resource files directory" + }, + "blast_db": { + "type": "string", + "format": "path", + "description": "Path to BLAST database directory" + }, + "blast_db_prefix": { + "type": "string", + "description": "BLAST database name prefix" + }, + "virus_index": { + "type": [ + "string", + "null" + ], + "default": null, + "description": "Path to a prebuilt vertebrate-infecting virus deacon index (.idx file). Providing this enables virus read enrichment." + }, + "virus_index_url": { + "type": [ + "string", + "null" + ], + "default": null, + "description": "URL to download a prebuilt vertebrate-infecting virus deacon index. Providing this enables virus read enrichment when virus_index is not set." + }, + "virus_reference_fasta": { + "type": [ + "string", + "null" + ], + "default": null, + "description": "Custom vertebrate-infecting virus FASTA to build and union with other virus indexes. Providing this enables virus read enrichment." + }, + "no_enrichment": { + "type": "boolean", + "default": false, + "description": "Disable target enrichment even when virus_index, virus_index_url, or virus_reference_fasta is set." + }, + "virus_kmer_size": { + "type": "integer", + "default": 31, + "minimum": 1, + "description": "K-mer size for building a custom virus enrichment index" + }, + "virus_window_size": { + "type": "integer", + "default": 1, + "minimum": 1, + "description": "Minimizer window size for building a custom virus enrichment index" + }, + "virus_abs_threshold": { + "type": "integer", + "default": 1, + "minimum": 1, + "description": "Minimum absolute minimizer hits for virus read enrichment" + }, + "virus_rel_threshold": { + "type": "number", + "default": 0.0, + "minimum": 0, + "maximum": 1, + "description": "Minimum relative proportion of minimizers for virus read enrichment (0.0-1.0)" + }, + "sourmash_ksize": { + "type": "integer", + "default": 31, + "minimum": 1, + "description": "K-mer size for experimental sourmash sketching" + }, + "sourmash_scaled": { + "type": "integer", + "default": 50, + "minimum": 1, + "description": "Scaled value for experimental sourmash sketching" + }, + "merge_pairs": { + "type": "boolean", + "default": true, + "description": "Merge overlapping paired-end reads before contig mapback. On by default; disable with --no-merge-pairs." + }, + "dedup": { + "type": "boolean", + "default": false, + "description": "Deduplicate sequencing reads (umbrella: enables both dedup_seq and dedup_pos)" + }, + "dedup_seq": { + "type": "boolean", + "default": false, + "description": "Sequence-based deduplication with clumpify (preprocessing)" + }, + "dedup_pos": { + "type": "boolean", + "default": false, + "description": "Positional deduplication with samtools markdup (after alignment)" + }, + "trim_adapters": { + "type": "boolean", + "description": "Trim Illumina adapters from reads" + }, + "filter_reads": { + "type": [ + "boolean", + "null" + ], + "default": null, + "description": "Apply quality and length filtering to reads" + }, + "filter_low_complexity_reads": { + "type": "boolean", + "default": false, + "description": "Filter reads below the minimum normalized 5-mer entropy" + }, + "min_read_quality_illumina": { + "type": "integer", + "default": 20, + "minimum": 0, + "description": "Minimum average quality score for Illumina reads" + }, + "min_read_quality_nanopore": { + "type": "integer", + "default": 12, + "minimum": 0, + "description": "Minimum average quality score for Nanopore reads" + }, + "min_read_length": { + "type": "integer", + "default": 50, + "minimum": 1, + "description": "Minimum read length to retain" + }, + "max_read_length": { + "type": [ + "integer", + "null" + ], + "default": null, + "minimum": 1, + "description": "Maximum read length to retain (no limit if not specified)" + }, + "min_read_entropy": { + "type": "number", + "default": 0.5, + "minimum": 0, + "maximum": 1, + "description": "Minimum normalized 5-mer entropy over 50-base windows" + }, + "host_index": { + "type": [ + "string", + "null" + ], + "default": null, + "description": "Path to a prebuilt host/contaminant index (.idx file). Providing this enables host depletion." + }, + "host_index_url": { + "type": [ + "string", + "null" + ], + "default": null, + "description": "URL to download a prebuilt host/contaminant index. Providing this enables host depletion when host_index is not set." + }, + "host_contaminants_fasta": { + "type": [ + "string", + "null" + ], + "default": null, + "description": "Custom contaminant FASTA to build and union with other host indexes. Providing this enables host depletion." + }, + "host_kmer_size": { + "type": "integer", + "default": 31, + "minimum": 1, + "description": "K-mer size for building a custom host/contaminant index" + }, + "host_window_size": { + "type": "integer", + "default": 15, + "minimum": 1, + "description": "Minimizer window size for building a custom host/contaminant index" + }, + "host_abs_threshold": { + "type": "integer", + "default": 2, + "minimum": 1, + "description": "Minimum absolute minimizer hits to classify as contaminant" + }, + "host_rel_threshold": { + "type": "number", + "default": 0.01, + "minimum": 0, + "maximum": 1, + "description": "Minimum relative proportion of minimizers (0.0-1.0)" + }, + "cutoff_percent": { + "type": "number", + "default": 0.001, + "minimum": 0, + "maximum": 1, + "description": "Minimum abundance threshold for reporting taxa (0-1)" + }, + "entropy": { + "type": "number", + "default": 0.9, + "minimum": 0, + "maximum": 1, + "description": "Entropy threshold for sequence complexity filtering (0-1)" + }, + "min_consecutive_bases": { + "type": "integer", + "default": 200, + "minimum": 1, + "description": "Minimum number of consecutive bases required" + }, + "qtrim": { + "type": "string", + "default": "t", + "description": "Quality trimming mode" + }, + "tax_stringency": { + "type": "number", + "default": 0.7, + "minimum": 0, + "maximum": 1, + "description": "Stringency for taxonomic classification (0-1)" + }, + "include_children": { + "type": "boolean", + "default": true, + "description": "Include child taxa in taxonomic analysis" + }, + "human_virus_families": { + "type": "array", + "items": { + "type": "string" + }, + "default": [ + "Adenoviridae", + "Anelloviridae", + "Arenaviridae", + "Arteriviridae", + "Astroviridae", + "Bornaviridae", + "Peribunyaviridae", + "Caliciviridae", + "Coronaviridae", + "Filoviridae", + "Flaviviridae", + "Hepadnaviridae", + "Hepeviridae", + "Orthoherpesviridae", + "Orthomyxoviridae", + "Papillomaviridae", + "Paramyxoviridae", + "Parvoviridae", + "Picobirnaviridae", + "Picornaviridae", + "Pneumoviridae", + "Polyomaviridae", + "Poxviridae", + "Sedoreoviridae", + "Retroviridae", + "Rhabdoviridae", + "Togaviridae", + "Kolmioviridae" + ], + "description": "List of virus family names to include in human virus analysis" + }, + "max_blast_targets": { + "type": "integer", + "default": 100, + "minimum": 1, + "description": "Maximum number of BLAST hits to consider before calling hits" + }, + "blast_retention_count": { + "type": "integer", + "default": 5, + "minimum": 1, + "description": "Number of top BLAST hits to retain" + }, + "refman_registry": { + "type": "string", + "format": "path", + "description": "Path to refman registry file for reference management" + }, + "monoimage": { + "type": "string", + "default": "nrminor/nvd:latest", + "description": "Container image to use for pipeline execution" + }, + "slack_enabled": { + "type": "boolean", + "default": false, + "description": "Enable Slack notifications for run completion" + }, + "slack_channel": { + "type": "string", + "pattern": "^C[A-Z0-9]+$", + "description": "Slack channel ID for notifications (e.g., 'C0123456789')" + }, + "taxonomy_dir": { + "type": [ + "string", + "null" + ], + "format": "path", + "description": "Explicit taxonomy database directory containing NCBI taxdump files. When unset, Python taxonomy helpers resolve the cache location." + }, + "taxonomy_mode": { + "type": [ + "string", + "null" + ], + "enum": [ + "read_only", + "missing", + null + ], + "default": null, + "description": "Pipeline taxonomy availability mode. null preserves legacy NVD_TAXONOMY_OFFLINE behavior; read_only never downloads or rebuilds; missing prepares taxonomy only when required files are absent." + }, + "taxonomy_refresh": { + "type": "string", + "enum": [ + "missing", + "stale", + "force" + ], + "default": "missing", + "description": "Admin taxonomy preflight refresh policy. missing downloads/builds only when required files are absent; stale refreshes when older than taxonomy_max_age_days; force always refreshes." + }, + "taxonomy_max_age_days": { + "type": "integer", + "minimum": 1, + "default": 90, + "description": "Freshness threshold in days for taxonomy stale refreshes and warnings." + } + }, + "additionalProperties": false +} diff --git a/uv.lock b/uv.lock index df4e97bb..47bfb397 100644 --- a/uv.lock +++ b/uv.lock @@ -1585,7 +1585,7 @@ wheels = [ [[package]] name = "nvd" -version = "3.4.1" +version = "3.5.0" source = { editable = "." } dependencies = [ { name = "altair" }, From a148979eb4e160eca185a4cb7c07bf50f853ee5b Mon Sep 17 00:00:00 2001 From: Nick Minor Date: Thu, 30 Jul 2026 19:32:20 -0500 Subject: [PATCH 06/10] fix: close out samples without enriched reads Carry validated input and retained counts on Deacon outputs, route zero-retention samples directly to completion, and keep downstream reporting aware of samples that cannot produce contigs. Cover the empty route through workflow and MultiQC boundaries. --- lib/py_nvd/multiqc_domains.py | 4 +- modules/deacon.nf | 19 +++++- pixi.lock | 58 ++++++++++++++++++ pyproject.toml | 1 + subworkflows/preprocess_reads.nf | 53 +++++++---------- tests/test_no_contig_routing.py | 92 +++++++++++++++++++++++++++++ tests/test_nvd_multiqc_reporting.py | 25 +++++--- workflows/nvd_main.nf | 1 + 8 files changed, 208 insertions(+), 45 deletions(-) diff --git a/lib/py_nvd/multiqc_domains.py b/lib/py_nvd/multiqc_domains.py index 8517badc..3909795a 100644 --- a/lib/py_nvd/multiqc_domains.py +++ b/lib/py_nvd/multiqc_domains.py @@ -183,9 +183,7 @@ def deacon_row(package: ReportPackage) -> DeaconRow: return DeaconRow( sample_id=receipt.sample_id, query_class=receipt.query_class, - status="complete_empty" - if stats.seqs_in == 0 and stats.bp_in == 0 - else "observed", + status="complete_empty" if stats.seqs_out == 0 else "observed", reads_in=stats.seqs_in, reads_retained=stats.seqs_out, reads_removed=stats.seqs_removed, diff --git a/modules/deacon.nf b/modules/deacon.nf index ecb111ec..ba75b0f3 100644 --- a/modules/deacon.nf +++ b/modules/deacon.nf @@ -167,7 +167,7 @@ process DEACON_ENRICH_TARGET_READS { tuple val(meta), path(read_files, stageAs: "reads??????/*"), path(deacon_idx), val(target_enrichment_enabled) output: - tuple val(meta.id), val(meta.platform), val(meta.deacon_read_structure), path("${meta.id}.target_enriched.fastq.gz"), emit: reads + tuple val(meta.id), val(meta.platform), val(meta.deacon_read_structure), path("${meta.id}.target_enriched.fastq.gz"), env('NVD_INPUT_READ_COUNT'), env('NVD_ENRICHED_READ_COUNT'), emit: reads tuple val(meta.id), path("${meta.id}.deacon_filter.json"), emit: stats script: @@ -187,6 +187,9 @@ process DEACON_ENRICH_TARGET_READS { --output ${meta.id}.target_enriched.fastq.gz \ ${deacon_idx} \ ${r1_files[0]} + + NVD_INPUT_READ_COUNT=\$(jq -er '.seqs_in | select(type == "number" and . >= 0 and . == floor)' ${meta.id}.deacon_filter.json) + NVD_ENRICHED_READ_COUNT=\$(jq -er '.seqs_out | select(type == "number" and . >= 0 and . == floor)' ${meta.id}.deacon_filter.json) """ else if (meta.read_mode == "paired" && r1_files.size() == 1 && r2_files.size() == 1) """ @@ -199,6 +202,9 @@ process DEACON_ENRICH_TARGET_READS { --output ${meta.id}.target_enriched.fastq.gz \ ${deacon_idx} \ ${r1_files[0]} ${r2_files[0]} + + NVD_INPUT_READ_COUNT=\$(jq -er '.seqs_in | select(type == "number" and . >= 0 and . == floor)' ${meta.id}.deacon_filter.json) + NVD_ENRICHED_READ_COUNT=\$(jq -er '.seqs_out | select(type == "number" and . >= 0 and . == floor)' ${meta.id}.deacon_filter.json) """ else if (meta.read_mode == "single") """ @@ -215,6 +221,9 @@ process DEACON_ENRICH_TARGET_READS { --rel-threshold ${params.virus_rel_threshold} \ --summary ${meta.id}.deacon_filter.json \ --output ${meta.id}.target_enriched.fastq.gz + + NVD_INPUT_READ_COUNT=\$(jq -er '.seqs_in | select(type == "number" and . >= 0 and . == floor)' ${meta.id}.deacon_filter.json) + NVD_ENRICHED_READ_COUNT=\$(jq -er '.seqs_out | select(type == "number" and . >= 0 and . == floor)' ${meta.id}.deacon_filter.json) """ else """ @@ -233,6 +242,9 @@ process DEACON_ENRICH_TARGET_READS { --rel-threshold ${params.virus_rel_threshold} \ --summary ${meta.id}.deacon_filter.json \ --output ${meta.id}.target_enriched.fastq.gz + + NVD_INPUT_READ_COUNT=\$(jq -er '.seqs_in | select(type == "number" and . >= 0 and . == floor)' ${meta.id}.deacon_filter.json) + NVD_ENRICHED_READ_COUNT=\$(jq -er '.seqs_out | select(type == "number" and . >= 0 and . == floor)' ${meta.id}.deacon_filter.json) """ } @@ -250,7 +262,7 @@ process DEACON_ENRICH_SRA_READS { tuple val(id), val(platform), val(run_accession), path(deacon_idx), val(target_enrichment_enabled) output: - tuple val(id), val(platform), path("${id}.sra_read_structure.txt"), path("${id}.target_enriched.fastq.gz"), emit: reads + tuple val(id), val(platform), path("${id}.sra_read_structure.txt"), path("${id}.target_enriched.fastq.gz"), env('NVD_INPUT_READ_COUNT'), env('NVD_ENRICHED_READ_COUNT'), emit: reads tuple val(id), path("${id}.deacon_filter.json"), emit: stats script: @@ -334,6 +346,9 @@ process DEACON_ENRICH_SRA_READS { --output '${id}.target_enriched.fastq.gz' \ ${deacon_idx} \ "\${deacon_inputs[@]}" + + NVD_INPUT_READ_COUNT=\$(jq -er '.seqs_in | select(type == "number" and . >= 0 and . == floor)' '${id}.deacon_filter.json') + NVD_ENRICHED_READ_COUNT=\$(jq -er '.seqs_out | select(type == "number" and . >= 0 and . == floor)' '${id}.deacon_filter.json') """ } diff --git a/pixi.lock b/pixi.lock index 49c3226c..825c7787 100644 --- a/pixi.lock +++ b/pixi.lock @@ -814,6 +814,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/icu-78.3-h4350ee2_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/isa-l-2.32.1-ha1e9b39_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/isl-0.26-imath32_h2e86a7b_101.conda + - conda: https://conda.anaconda.org/conda-forge/osx-64/jq-1.8.2-ha3d0635_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/kaleido-core-0.2.1-h0d85af4_0.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-64/kiwisolver-1.5.0-py313h224b87c_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/krb5-1.22.2-h3ddfcb2_1.conda @@ -880,6 +881,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/mpi-1.0-openmpi.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-64/ncurses-6.6-hcc0dc9a_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/numpy-2.5.2-py313hb870fc3_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-64/oniguruma-6.9.10-ha3d0635_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/openjdk-23.0.2-h18c9476_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/openjpeg-2.5.4-h52bb76a_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/openmpi-4.1.6-h7406208_101.conda @@ -1101,6 +1103,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/icu-78.3-hc7cc350_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/isa-l-2.31.1-h5505292_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/isl-0.26-imath32_h347afa1_101.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/jq-1.8.2-h1a92334_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/kaleido-core-0.2.1-h27ca646_0.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/kiwisolver-1.5.0-py313h2af2deb_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/krb5-1.22.2-hfd3d5f3_1.conda @@ -1167,6 +1170,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mpi-1.0-openmpi.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ncurses-6.6-he64c551_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-2.5.2-py313hce9b930_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/oniguruma-6.9.10-h1a92334_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjdk-23.0.2-hfb9339a_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjpeg-2.5.4-hd9e9057_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openmpi-4.1.6-h526c993_101.conda @@ -2216,6 +2220,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/icu-78.3-h4350ee2_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/isa-l-2.32.1-ha1e9b39_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/isl-0.26-imath32_h2e86a7b_101.conda + - conda: https://conda.anaconda.org/conda-forge/osx-64/jq-1.8.2-ha3d0635_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/kaleido-core-0.2.1-h0d85af4_0.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-64/kiwisolver-1.5.0-py312hb1dc2e7_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/krb5-1.22.2-h3ddfcb2_1.conda @@ -2281,6 +2286,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/mpi-1.0-openmpi.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-64/ncurses-6.6-hcc0dc9a_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/numpy-2.5.2-py312h746d82c_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-64/oniguruma-6.9.10-ha3d0635_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/openjdk-23.0.2-h18c9476_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/openjpeg-2.5.4-h52bb76a_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/openmpi-4.1.6-h7406208_101.conda @@ -2602,6 +2608,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/icu-78.3-hc7cc350_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/isa-l-2.31.1-h5505292_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/isl-0.26-imath32_h347afa1_101.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/jq-1.8.2-h1a92334_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/kaleido-core-0.2.1-h27ca646_0.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/kiwisolver-1.5.0-py312h3093aea_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/krb5-1.22.2-hfd3d5f3_1.conda @@ -2667,6 +2674,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/mpi-1.0-openmpi.tar.bz2 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ncurses-6.6-he64c551_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-2.5.2-py312ha003a3f_0.conda + - conda: https://conda.anaconda.org/conda-forge/osx-arm64/oniguruma-6.9.10-h1a92334_1.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjdk-23.0.2-hfb9339a_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjpeg-2.5.4-hd9e9057_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openmpi-4.1.6-h526c993_101.conda @@ -12078,6 +12086,18 @@ packages: - isl 0.26.* size: 894410 timestamp: 1680649639107 +- conda: https://conda.anaconda.org/conda-forge/osx-64/jq-1.8.2-ha3d0635_1.conda + sha256: 891be95d8fd8110b6f8e1673d029761a540f099e9bc62e90a35b03ee17f2629e + md5: 99ea337e58c6216a04473c8f52815435 + depends: + - oniguruma 6.9.* + - __osx >=11.0 + - oniguruma >=6.9.10,<6.10.0a0 + license: MIT + purls: [] + run_exports: {} + size: 344199 + timestamp: 1786435076751 - conda: https://conda.anaconda.org/conda-forge/osx-64/kaleido-core-0.2.1-h0d85af4_0.tar.bz2 sha256: 86ad734be7c1f02810c6048506223c58bc9f55326c091704ebacd2c0e57e937c md5: 2d0e8226402dbd6ed0d9134be2aa2bcb @@ -13154,6 +13174,19 @@ packages: - numpy >=1.25,<3 size: 8236483 timestamp: 1786330709895 +- conda: https://conda.anaconda.org/conda-forge/osx-64/oniguruma-6.9.10-ha3d0635_1.conda + sha256: bd36b82b3553a458535068a91ce5e9632c7b4ddf6c122e6e58ccb31fc9a23ded + md5: 4abc9f10ed9f563b4e5c235e115ba9f8 + depends: + - __osx >=11.0 + license: BSD-2-Clause + license_family: BSD + purls: [] + run_exports: + weak: + - oniguruma >=6.9.10,<6.10.0a0 + size: 263274 + timestamp: 1786354383678 - conda: https://conda.anaconda.org/conda-forge/osx-64/openjdk-23.0.2-h18c9476_2.conda sha256: 00a2381baecf707b4cb6ed74c4d91ede955e8f6f859e7caba35aa11720af8dd6 md5: e8b746e8947c0dab5aeb7a26b8bce2ad @@ -14788,6 +14821,18 @@ packages: - isl 0.26.* size: 819937 timestamp: 1680649567633 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/jq-1.8.2-h1a92334_1.conda + sha256: 2d345464da308424f15322fb848b0f845291fdf2383f709866e94825d722f550 + md5: be1ad8bb4f91519d2a3eb2bcb6d9bc0b + depends: + - oniguruma 6.9.* + - __osx >=11.0 + - oniguruma >=6.9.10,<6.10.0a0 + license: MIT + purls: [] + run_exports: {} + size: 350634 + timestamp: 1786434970720 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/kaleido-core-0.2.1-h27ca646_0.tar.bz2 sha256: 4d25abd1b9c2934d1f5ef87b9b481b6ca6fd84f712b44ded2c584569a58a09b7 md5: 5ed406a3805fc24876d90575129a9f87 @@ -15869,6 +15914,19 @@ packages: - numpy >=1.25,<3 size: 7088694 timestamp: 1786330629373 +- conda: https://conda.anaconda.org/conda-forge/osx-arm64/oniguruma-6.9.10-h1a92334_1.conda + sha256: 92a6ca62564165da0c1d6c321555b8cd3e5a660512ee0466b89c46e3637a5bf7 + md5: ffe6286c38000935f186202d469aab0a + depends: + - __osx >=11.0 + license: BSD-2-Clause + license_family: BSD + purls: [] + run_exports: + weak: + - oniguruma >=6.9.10,<6.10.0a0 + size: 256169 + timestamp: 1786354290452 - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjdk-23.0.2-hfb9339a_2.conda sha256: f16143181808731b74e03b191a1d859cfb6dc7ed066f9d567808c438e9770f85 md5: 6e330bfe359f78e2313eaecf8ca1529c diff --git a/pyproject.toml b/pyproject.toml index 5149193f..2627dd8e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -62,6 +62,7 @@ vsearch = ">=2.31.0,<3" # System Dependencies # ----------------------------------------------------------------------------------- # compilers = ">=1.9.0,<2" +jq = ">=1.8.2,<2" nextflow = ">=25.10.2,<26" openssl = ">=3.5.0,<4" pkg-config = ">=0.29.2,<0.30" diff --git a/subworkflows/preprocess_reads.nf b/subworkflows/preprocess_reads.nf index 0184f8f4..8e7b5f47 100644 --- a/subworkflows/preprocess_reads.nf +++ b/subworkflows/preprocess_reads.nf @@ -71,9 +71,9 @@ workflow PREPROCESS_READS { // Step 2: Inlined preprocessing on target-enriched reads // ------------------------------------------------------------------------- ch_sra_target_reads = DEACON_ENRICH_SRA_READS.out.reads - .map { sample_id, platform, read_structure_file, reads -> + .map { sample_id, platform, read_structure_file, reads, input_read_count, enriched_read_count -> def read_structure = read_structure_file.toFile().text.trim() - tuple(sample_id, platform, read_structure, reads) + tuple(sample_id, platform, read_structure, reads, input_read_count, enriched_read_count) } ch_target_reads = DEACON_ENRICH_TARGET_READS.out.reads @@ -82,39 +82,30 @@ workflow PREPROCESS_READS { ch_target_enrichment_stats = DEACON_ENRICH_TARGET_READS.out.stats .mix(DEACON_ENRICH_SRA_READS.out.stats) - // Extract input and retained counts from Deacon's summary. A malformed - // summary is a process failure, never an empty-sample completion. - ch_enrichment_counts = ch_target_enrichment_stats - .map { sample_id, json_file -> - def summary = new groovy.json.JsonSlurper().parse(json_file.toFile()) - def input_count = summary.seqs_in - def retained_count = summary.seqs_out - assert input_count instanceof Number && input_count >= 0 && input_count == input_count.toLong() - assert retained_count instanceof Number && retained_count >= 0 && retained_count == retained_count.toLong() - tuple(sample_id, input_count.toLong().toString(), retained_count.toLong()) + // Deacon reports counts across all input and retained reads. Carry those + // scalar values on the channel so workflow routing never parses JSON. + ch_read_counts = ch_target_reads + .map { sample_id, _platform, _read_structure, _reads, input_read_count, _enriched_read_count -> + tuple(sample_id, input_read_count) } - ch_read_counts = ch_enrichment_counts - .map { sample_id, input_count, _retained_count -> tuple(sample_id, input_count) } - - ch_target_reads_by_retention = ch_target_reads - .join(ch_enrichment_counts, by: 0) - .branch { _sample_id, _platform, _read_structure, _reads, _input_count, retained_count -> - retained: retained_count > 0 - complete_empty: true - } + ch_target_reads_by_retention = ch_target_reads.branch { _id, _platform, _read_structure, _reads, _input_read_count, enriched_read_count -> + retained: enriched_read_count != "0" + complete_empty: true + } - ch_target_reads = ch_target_reads_by_retention.retained - .map { sample_id, platform, read_structure, reads, _input_count, _retained_count -> - tuple(sample_id, platform, read_structure, reads) - } + ch_complete_empty_samples = ch_target_reads_by_retention.complete_empty.map { + sample_id, platform, _read_structure, _reads, _input_read_count, _enriched_read_count -> + log.debug "nvd.contig_route sample_id=${sample_id} platform=${platform} outcome=no_contigs stage=target_enrichment" + tuple(sample_id, platform) + } - ch_complete_empty_samples = ch_target_reads_by_retention.complete_empty - .map { sample_id, platform, _read_structure, _reads, _input_count, _retained_count -> - tuple(sample_id, platform) - } + ch_retained_target_reads = ch_target_reads_by_retention.retained.map { + sample_id, platform, read_structure, reads, _input_read_count, _enriched_read_count -> + tuple(sample_id, platform, read_structure, reads) + } - ch_target_reads_by_layout = ch_target_reads.branch { _id, platform, read_structure, _reads -> + ch_target_reads_by_layout = ch_retained_target_reads.branch { _id, platform, read_structure, _reads -> mergeable: platform == "illumina" && read_structure == "interleaved" other: true } @@ -176,7 +167,7 @@ workflow PREPROCESS_READS { ? ch_merged_read_batches .mix(ch_unmerged_read_batches) .mix(ch_nonmergeable_read_batches) - : ch_target_reads.map { sample_id, platform, read_structure, reads -> + : ch_retained_target_reads.map { sample_id, platform, read_structure, reads -> tuple( [ id: sample_id, diff --git a/tests/test_no_contig_routing.py b/tests/test_no_contig_routing.py index c68811dd..641cb086 100644 --- a/tests/test_no_contig_routing.py +++ b/tests/test_no_contig_routing.py @@ -16,6 +16,7 @@ PREPARE_BLAST_QUERIES = ROOT / "subworkflows" / "prepare_blast_queries" SHORT_READ_ASSEMBLY = ROOT / "subworkflows" / "short_read_denovo_assembly" LONG_READ_ENSEMBLE = ROOT / "subworkflows" / "long_read_denovo_ensembly" +PREPROCESS_READS = ROOT / "subworkflows" / "preprocess_reads" FASTX_MODULE = ROOT / "modules" / "fastx" @@ -310,6 +311,97 @@ def write_deacon_fake(bin_dir: Path) -> None: ) +def test_empty_target_enrichment_closes_sample_before_preprocessing( + tmp_path: Path, +) -> None: + bin_dir = tmp_path / "bin" + bin_dir.mkdir() + write_deacon_fake(bin_dir) + + lib = tmp_path / "lib" + lib.mkdir() + for filename in ("NvdReporting.groovy", "NvdUtils.groovy"): + shutil.copy2(ROOT / "lib" / filename, lib / filename) + + reads = write_fastq(tmp_path / "reads.fastq.gz", "source-read") + target_index = tmp_path / "target.idx" + target_index.touch() + + workflow = tmp_path / "main.nf" + workflow.write_text( + f"""\ +nextflow.enable.dsl = 2 + +include {{ PREPROCESS_READS }} from '{PREPROCESS_READS}' + +params.skip_fastqc = true +params.no_enrichment = false +params.virus_index = '{target_index}' +params.virus_index_url = null +params.virus_reference_fasta = null +params.virus_abs_threshold = 1 +params.virus_rel_threshold = 0.0 +params.merge_pairs = false +params.dedup = false +params.dedup_seq = false +params.trim_adapters = false +params.host_index = null +params.host_index_url = null +params.host_contaminants_fasta = null +params.filter_reads = false +params.filter_low_complexity_reads = false +params.min_read_quality_illumina = 20 +params.min_read_quality_nanopore = 12 +params.min_read_length = 50 +params.min_consecutive_bases = 200 +params.max_concurrent_downloads = 1 + +workflow {{ + PREPROCESS_READS( + Channel.of(tuple( + [ + id: 'empty_enrichment', + platform: 'illumina', + read_mode: 'single', + r1_count: 1, + deacon_read_structure: 'single', + ], + [file('{reads}')], + )), + Channel.empty(), + ) + + PREPROCESS_READS.out.complete_empty_samples.view {{ sample_id, platform -> + "COMPLETE_EMPTY: ${{sample_id}}:${{platform}}" + }} + PREPROCESS_READS.out.read_counts.view {{ sample_id, count -> + "INPUT_READS: ${{sample_id}}:${{count}}" + }} + ch_unexpected_read_outputs = PREPROCESS_READS.out.reads + .map {{ sample_id, _platform, _read_structure, _reads -> sample_id }} + .mix(PREPROCESS_READS.out.read_batches.map {{ sample_id, _platform, _read_structure, _query_class, _reads -> sample_id }}) + .mix(PREPROCESS_READS.out.profiled_read_batches.map {{ meta, _reads, _profile, _histogram -> meta.id }}) + .mix(PREPROCESS_READS.out.profiled_batches_by_sample.map {{ meta, _batches -> meta.id }}) + .mix(PREPROCESS_READS.out.paired_reads_for_mapback.map {{ sample_id, _platform, _overlap_reads, _single_reads -> sample_id }}) + .mix(PREPROCESS_READS.out.single_reads_for_mapback.map {{ sample_id, _platform, _read_structure, _reads -> sample_id }}) + + ch_unexpected_read_outputs.view {{ sample_id -> + "UNEXPECTED_READ_OUTPUT: ${{sample_id}}" + }} +}} +""", + encoding="utf-8", + ) + + completed = run_nextflow(workflow, bin_dir=bin_dir) + diagnostics = f"stdout:\n{completed.stdout}\nstderr:\n{completed.stderr}" + + assert completed.returncode == 0, diagnostics + assert "COMPLETE_EMPTY: empty_enrichment:illumina" in completed.stdout + assert "INPUT_READS: empty_enrichment:1" in completed.stdout + assert "UNEXPECTED_READ_OUTPUT:" not in completed.stdout + + def test_no_contig_samples_route_full_reads_without_mapback(tmp_path: Path) -> None: """Read queries are a default capability; no --experimental flag is needed.""" bin_dir = tmp_path / "bin" diff --git a/tests/test_nvd_multiqc_reporting.py b/tests/test_nvd_multiqc_reporting.py index 3b6301ba..ae0db183 100644 --- a/tests/test_nvd_multiqc_reporting.py +++ b/tests/test_nvd_multiqc_reporting.py @@ -9,6 +9,7 @@ import subprocess from pathlib import Path +import pytest import yaml ROOT = Path(__file__).resolve().parents[1] @@ -386,8 +387,14 @@ def test_fastqc_unit_expansion_for_supported_bundle_shapes(tmp_path: Path) -> No assert "SKIP: no_fastqc_units" in completed.stdout +@pytest.mark.parametrize( + ("seqs_out", "expected_status"), + [(8, "observed"), (0, "complete_empty")], +) def test_target_enrichment_package_crosses_explicit_process_boundary( tmp_path: Path, + seqs_out: int, + expected_status: str, ) -> None: bin_dir = tmp_path / "bin" bin_dir.mkdir() @@ -402,15 +409,15 @@ def test_target_enrichment_package_crosses_explicit_process_boundary( { "version": "deacon-test", "seqs_in": 10, - "seqs_out": 8, - "seqs_removed": 2, - "seqs_out_proportion": 0.8, - "seqs_removed_proportion": 0.2, + "seqs_out": seqs_out, + "seqs_removed": 10 - seqs_out, + "seqs_out_proportion": seqs_out / 10, + "seqs_removed_proportion": (10 - seqs_out) / 10, "bp_in": 1000, - "bp_out": 800, - "bp_removed": 200, - "bp_out_proportion": 0.8, - "bp_removed_proportion": 0.2, + "bp_out": seqs_out * 100, + "bp_removed": (10 - seqs_out) * 100, + "bp_out_proportion": seqs_out / 10, + "bp_removed_proportion": (10 - seqs_out) / 10, }, ), encoding="utf-8", @@ -446,7 +453,7 @@ def test_target_enrichment_package_crosses_explicit_process_boundary( sections = sorted(tmp_path.glob("work/**/nvd_target_enrichment_mqc.yaml")) assert sections, diagnostics section_text = sections[-1].read_text(encoding="utf-8") - assert "status: observed" in section_text + assert f"status: {expected_status}" in section_text assert "reads_in: 10" in section_text assert "\n sample_A:\n" in section_text assert "row_0001" not in section_text diff --git a/workflows/nvd_main.nf b/workflows/nvd_main.nf index 17201320..663169d3 100644 --- a/workflows/nvd_main.nf +++ b/workflows/nvd_main.nf @@ -100,6 +100,7 @@ workflow NVD_MAIN { .mix(LONG_READ_DENOVO_ENSEMBLY.out.no_contigs) .mix(PROCESS_CONTIGS.out.no_contigs) .mix(ch_assembly_disabled) + .mix(PREPROCESS_READS.out.complete_empty_samples) ch_run_context = COMPUTE_RUN_CONTEXT.out.run_context ch_taxonomy_dir = ENSURE_TAXONOMY.out.taxonomy_dir From 4a5466abf4d2cf2e1f461170f033356af45beccf Mon Sep 17 00:00:00 2001 From: Nick Minor Date: Tue, 11 Aug 2026 10:19:15 -0500 Subject: [PATCH 07/10] feat: validate paired Deacon inputs Upgrade Deacon to 0.16 and use its explicit interleaved-input mode for paired streams and files. Add an opt-in NVD check-pairs parameter that asks Deacon to reject mismatched CASAVA or /1 and /2 record names without changing the default compatibility policy. Propagate the validation setting through direct, streamed, and SRA-backed paired filtering while leaving single-read, contig, and post-merge per-record filtering unchanged. --- bin/stream_fastqs_to_deacon.py | 7 + bin/test_stream_fastqs_to_deacon.py | 87 ++++++- lib/py_nvd/cli/commands/preset.py | 6 + lib/py_nvd/cli/commands/run.py | 10 + lib/py_nvd/cli/test_cli_smoke.py | 21 ++ lib/py_nvd/models.py | 7 + lib/py_nvd/params.py | 1 + lib/py_nvd/test_models.py | 24 ++ modules/deacon.nf | 24 +- nextflow.config | 2 + pixi.lock | 81 ++++--- pyproject.toml | 2 +- schemas/nvd-params.v3.5.0.schema.json | 5 + tests/test_deacon_pairing.py | 65 +++++- tests/test_target_enrichment.py | 316 ++++++++++++++++++++++++++ 15 files changed, 604 insertions(+), 54 deletions(-) diff --git a/bin/stream_fastqs_to_deacon.py b/bin/stream_fastqs_to_deacon.py index 467d598e..7595a297 100755 --- a/bin/stream_fastqs_to_deacon.py +++ b/bin/stream_fastqs_to_deacon.py @@ -42,6 +42,7 @@ class DeaconStreamConfig: rel_threshold: float deplete: bool deacon_bin: str + check_pairs: bool = False @dataclass(frozen=True) @@ -66,6 +67,7 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace: parser.add_argument("--abs-threshold", type=int, default=1) parser.add_argument("--rel-threshold", type=float, default=0.0) parser.add_argument("--deplete", action="store_true") + parser.add_argument("--check-pairs", action="store_true") parser.add_argument("--deacon-bin", default="deacon") return parser.parse_args(argv) @@ -94,6 +96,7 @@ def config_from_args(args: argparse.Namespace) -> DeaconStreamConfig: rel_threshold=args.rel_threshold, deplete=args.deplete, deacon_bin=args.deacon_bin, + check_pairs=args.check_pairs, ) @@ -144,6 +147,8 @@ def validate_config( if single == paired: message = "Provide either --reads-list for single-end input or both --r1-list and --r2-list for paired input" raise StreamError(message) + if config.check_pairs and not paired: + raise StreamError("check-pairs requires paired input") if single: validate_files("reads", config.reads) validate_sample_compression(config.reads) @@ -204,9 +209,11 @@ def stream_files_to_fifo( def deacon_command(config: DeaconStreamConfig, inputs: tuple[Path, ...]) -> list[str]: + pair_args = ["--check-pairs"] if config.check_pairs and config.r1 else [] command = [ config.deacon_bin, "filter", + *pair_args, "--threads", str(config.threads), "--abs-threshold", diff --git a/bin/test_stream_fastqs_to_deacon.py b/bin/test_stream_fastqs_to_deacon.py index 3c2fcbc5..f9d42ad0 100644 --- a/bin/test_stream_fastqs_to_deacon.py +++ b/bin/test_stream_fastqs_to_deacon.py @@ -18,6 +18,7 @@ DeaconStreamConfig, StreamError, config_from_args, + deacon_command, parse_args, run_deacon_stream, ) @@ -53,7 +54,7 @@ def command_io(command: list[str]) -> tuple[Path, Path, Path, tuple[Path, ...]]: summary = Path(command[command.index("--summary") + 1]) index = command.index("filter") + 1 while command[index].startswith("--"): - index += 2 + index += 1 if command[index] in {"--check-pairs", "--deplete"} else 2 return ( output, summary, @@ -194,6 +195,7 @@ def config( reads: tuple[Path, ...] = (), r1: tuple[Path, ...] = (), r2: tuple[Path, ...] = (), + check_pairs: bool = False, ) -> DeaconStreamConfig: index = tmp_path / "index.dcn" index.write_text("fake index", encoding="utf-8") @@ -210,6 +212,7 @@ def config( rel_threshold=0.0, deplete=False, deacon_bin=str(deacon_bin), + check_pairs=check_pairs, ) @@ -301,6 +304,34 @@ def test_paired_bundle_streams_both_mates_in_order(tmp_path: Path) -> None: assert output.index("@L1_A/2") < output.index("@L2_B/2") +def test_deacon_command_forwards_opt_in_pair_validation(tmp_path: Path) -> None: + """Only an enabled paired helper command asks Deacon to check names.""" + r1 = write_fastq(tmp_path / "r1.fastq.gz", ["pair/1"]) + r2 = write_fastq(tmp_path / "r2.fastq.gz", ["pair/2"]) + paired = config(tmp_path, r1=(r1,), r2=(r2,), check_pairs=True) + unchecked = config(tmp_path, r1=(r1,), r2=(r2,)) + + assert "--check-pairs" in deacon_command( + paired, + (Path("r1.fastq.gz"), Path("r2.fastq.gz")), + ) + assert "--check-pairs" not in deacon_command( + unchecked, + (Path("r1.fastq.gz"), Path("r2.fastq.gz")), + ) + + +def test_pair_validation_rejects_single_input_mode(tmp_path: Path) -> None: + """The helper cannot silently ignore pair validation on single reads.""" + reads = write_fastq(tmp_path / "reads.fastq.gz", ["single"]) + + with pytest.raises(StreamError, match="check-pairs requires paired input"): + run_deacon_stream( + config(tmp_path, reads=(reads,), check_pairs=True), + runner=copying_deacon_runner(), + ) + + def test_paired_gzip_fifo_preserves_record_ordinality_across_lanes( tmp_path: Path, ) -> None: @@ -562,6 +593,35 @@ def test_cli_list_files_are_supported( assert capsys.readouterr().err == "" +def test_cli_accepts_pair_validation_for_paired_lists(tmp_path: Path) -> None: + """The helper exposes Deacon pair-name validation for its paired path.""" + r1 = (write_fastq(tmp_path / "r1.fastq.gz", ["pair/1"]),) + r2 = (write_fastq(tmp_path / "r2.fastq.gz", ["pair/2"]),) + r1_list = write_path_list(tmp_path / "r1.txt", r1) + r2_list = write_path_list(tmp_path / "r2.txt", r2) + argv = [ + "--sample-id", + "S1", + "--index", + str(tmp_path / "index.dcn"), + "--output", + str(tmp_path / "out.fastq.gz"), + "--summary", + str(tmp_path / "summary.json"), + "--r1-list", + str(r1_list), + "--r2-list", + str(r2_list), + "--check-pairs", + ] + + loaded = config_from_args(parse_args(argv)) + + assert loaded.r1 == r1 + assert loaded.r2 == r2 + assert loaded.check_pairs is True + + @pytest.fixture(scope="module") def deacon_bin() -> str: deacon = shutil.which("deacon") @@ -643,6 +703,7 @@ def real_deacon_config( # noqa: PLR0913 reads: tuple[Path, ...] = (), r1: tuple[Path, ...] = (), r2: tuple[Path, ...] = (), + check_pairs: bool = False, ) -> DeaconStreamConfig: return DeaconStreamConfig( sample_id="real", @@ -657,6 +718,7 @@ def real_deacon_config( # noqa: PLR0913 rel_threshold=0.0, deplete=False, deacon_bin=deacon, + check_pairs=check_pairs, ) @@ -712,6 +774,7 @@ def test_real_deacon_reads_paired_fifo_stream(tmp_path: Path, deacon_bin: str) - index, r1=(r1_lane1, r1_lane2), r2=(r2_lane1, r2_lane2), + check_pairs=True, ), ) @@ -723,6 +786,28 @@ def test_real_deacon_reads_paired_fifo_stream(tmp_path: Path, deacon_bin: str) - assert output.index("@L1_A/2") < output.index("@L2_B/1") +def test_real_deacon_pair_validation_rejects_mismatched_names( + tmp_path: Path, + deacon_bin: str, +) -> None: + """Opt-in validation surfaces positional pairs with different names.""" + index = build_deacon_index(tmp_path, deacon_bin) + r1 = write_fastq(tmp_path / "R1.fastq.gz", ["left/1"]) + r2 = write_fastq(tmp_path / "R2.fastq.gz", ["right/2"]) + + with pytest.raises(StreamError, match="deacon filter failed"): + run_deacon_stream( + real_deacon_config( + tmp_path, + deacon_bin, + index, + r1=(r1,), + r2=(r2,), + check_pairs=True, + ), + ) + + def test_real_deacon_streams_larger_paired_bundle_in_pair_order( tmp_path: Path, deacon_bin: str, diff --git a/lib/py_nvd/cli/commands/preset.py b/lib/py_nvd/cli/commands/preset.py index 5dba9f45..b05eafb4 100644 --- a/lib/py_nvd/cli/commands/preset.py +++ b/lib/py_nvd/cli/commands/preset.py @@ -181,6 +181,11 @@ def preset_register( "--merge-pairs/--no-merge-pairs", help="Merge overlapping paired-end reads before contig mapback", ), + check_pairs: bool | None = typer.Option( + None, + "--check-pairs/--no-check-pairs", + help="Ask Deacon to validate paired FASTQ record names", + ), host_index: Path | None = typer.Option( None, "--host-index", @@ -297,6 +302,7 @@ def preset_register( "dedup_pos": dedup_pos, "trim_adapters": trim_adapters, "merge_pairs": merge_pairs, + "check_pairs": check_pairs, "host_index": host_index, "host_index_url": host_index_url, "host_contaminants_fasta": host_contaminants_fasta, diff --git a/lib/py_nvd/cli/commands/run.py b/lib/py_nvd/cli/commands/run.py index 1d421f04..ae6bdc71 100644 --- a/lib/py_nvd/cli/commands/run.py +++ b/lib/py_nvd/cli/commands/run.py @@ -351,6 +351,15 @@ def run( help="Merge overlapping paired-end reads before contig mapback (default: on)", rich_help_panel=PANEL_PREPROCESSING, ), + check_pairs: bool | None = typer.Option( + None, + "--check-pairs/--no-check-pairs", + help=( + "Ask Deacon to validate paired FASTQ record names and fail on " + "mismatches; accepts CASAVA names or /1 and /2 suffixes" + ), + rich_help_panel=PANEL_PREPROCESSING, + ), host_index: Path | None = typer.Option( None, "--host-index", @@ -623,6 +632,7 @@ def run( "dedup_pos": dedup_pos, "trim_adapters": trim_adapters, "merge_pairs": merge_pairs, + "check_pairs": check_pairs, "host_index": host_index, "host_index_url": host_index_url, "host_contaminants_fasta": host_contaminants_fasta, diff --git a/lib/py_nvd/cli/test_cli_smoke.py b/lib/py_nvd/cli/test_cli_smoke.py index 3733e717..f55045ce 100644 --- a/lib/py_nvd/cli/test_cli_smoke.py +++ b/lib/py_nvd/cli/test_cli_smoke.py @@ -227,6 +227,27 @@ def test_run_help_describes_moderate_read_entropy_default() -> None: assert "default: 0.5" in result.output +def test_run_accepts_check_pairs_flag(tmp_path: Path) -> None: + """The pair-validation CLI flag maps to the underscore Nextflow param.""" + samplesheet = tmp_path / "samples.csv" + samplesheet.write_text("sample_id,srr,platform,fastq1,fastq2\n", encoding="utf-8") + + result = runner.invoke( + app, + [ + "run", + "--samplesheet", + str(samplesheet), + "--check-pairs", + "--dry-run", + ], + ) + + assert result.exit_code == 0, result.output + assert "--check_pairs true" in result.output + assert "--check-pairs" not in result.output + + def test_samplesheet_generate_sanitizes_illumina_ids(tmp_path: Path) -> None: fastq_dir = tmp_path / "fastqs" fastq_dir.mkdir() diff --git a/lib/py_nvd/models.py b/lib/py_nvd/models.py index 99c7e592..b9ff48d5 100644 --- a/lib/py_nvd/models.py +++ b/lib/py_nvd/models.py @@ -254,6 +254,13 @@ def reject_removed_params(cls, data: Any) -> Any: # noqa: ANN401 # pydantic ho ), json_schema_extra={"category": "Preprocessing"}, ) + check_pairs: bool = Field( + default=False, + description=( + "Ask Deacon to validate paired FASTQ record names and fail on mismatches" + ), + json_schema_extra={"category": "Preprocessing"}, + ) dedup: bool = Field( default=False, description="Deduplicate reads (umbrella: enables both dedup_seq and dedup_pos)", diff --git a/lib/py_nvd/params.py b/lib/py_nvd/params.py index c146739b..c7862056 100644 --- a/lib/py_nvd/params.py +++ b/lib/py_nvd/params.py @@ -279,6 +279,7 @@ def _yaml_analysis_section( "dedup_pos", "trim_adapters", "merge_pairs", + "check_pairs", "filter_reads", "filter_low_complexity_reads", "min_read_entropy", diff --git a/lib/py_nvd/test_models.py b/lib/py_nvd/test_models.py index fd20105f..a06abb47 100644 --- a/lib/py_nvd/test_models.py +++ b/lib/py_nvd/test_models.py @@ -1,5 +1,6 @@ """Tests for retained py_nvd parameter models.""" +import json from pathlib import Path import pytest @@ -13,6 +14,8 @@ ) from py_nvd.params import load_params_file +ROOT = Path(__file__).resolve().parents[2] + class TestNvdParamsInstantiation: """Tests for basic NvdParams instantiation.""" @@ -481,6 +484,15 @@ def test_skip_stage_params(self) -> None: assert "--skip-blast" not in cmd assert "--skip-fastqc" not in cmd + def test_check_pairs_param(self) -> None: + """Pair validation reaches Nextflow with underscore naming.""" + p = NvdParams(check_pairs=True) + cmd = p.to_nextflow_args(Path("/pipeline")) + + check_pairs_idx = cmd.index("--check_pairs") + assert cmd[check_pairs_idx + 1] == "true" + assert "--check-pairs" not in cmd + def test_default_cutoff_percent(self) -> None: """Default cutoff_percent matches nextflow.config.""" assert NvdParams().cutoff_percent == 0.001 @@ -540,6 +552,18 @@ def test_default_merge_pairs(self) -> None: """Pair merging is on by default; --no-merge-pairs is the opt-out.""" assert NvdParams().merge_pairs is True + def test_default_check_pairs(self) -> None: + """Model, published schema, and Nextflow keep validation opt-in.""" + assert NvdParams().check_pairs is False + schema = json.loads( + (ROOT / "schemas" / "nvd-params.v3.5.0.schema.json").read_text( + encoding="utf-8", + ), + ) + assert schema["properties"]["check_pairs"]["default"] is False + config = (ROOT / "nextflow.config").read_text(encoding="utf-8") + assert "check_pairs = false" in config + def test_default_low_complexity_read_filter(self) -> None: """Low-complexity read filtering is opt-in with a dormant threshold.""" params = NvdParams() diff --git a/modules/deacon.nf b/modules/deacon.nf index ba75b0f3..ec693b29 100644 --- a/modules/deacon.nf +++ b/modules/deacon.nf @@ -105,22 +105,24 @@ process DEACON_DEPLETE { tuple val(meta), path("${meta.id}.${meta.query_class}.deacon.json"), emit: stats script: - def input_stream = meta.read_structure == "interleaved" - ? "gzip -dc ${reads} | " + def interleaved_arg = meta.read_structure == "interleaved" ? "--interleaved" : "" + def check_pairs_arg = meta.read_structure == "interleaved" && params.check_pairs + ? "--check-pairs" : "" - def filter_inputs = meta.read_structure == "interleaved" ? "- -" : "${reads}" """ set -euo pipefail - ${input_stream}deacon filter \ + deacon filter \ --deplete \ + ${interleaved_arg} \ + ${check_pairs_arg} \ --threads ${task.cpus} \ --abs-threshold ${params.host_abs_threshold} \ --rel-threshold ${params.host_rel_threshold} \ --summary ${meta.id}.${meta.query_class}.deacon.json \ --output ${meta.id}.${meta.query_class}.depleted.fastq.gz \ ${index} \ - ${filter_inputs} + ${reads} """ } @@ -176,6 +178,7 @@ process DEACON_ENRICH_TARGET_READS { def r1_files = files.take(r1_count) def r2_files = files.drop(r1_count) def deplete_arg = target_enrichment_enabled ? "" : "--deplete" + def check_pairs_arg = params.check_pairs ? "--check-pairs" : "" if (meta.read_mode == "single" && r1_files.size() == 1) """ deacon filter \ @@ -195,6 +198,7 @@ process DEACON_ENRICH_TARGET_READS { """ deacon filter \ ${deplete_arg} \ + ${check_pairs_arg} \ --threads ${task.cpus} \ --abs-threshold ${params.virus_abs_threshold} \ --rel-threshold ${params.virus_rel_threshold} \ @@ -235,6 +239,7 @@ process DEACON_ENRICH_TARGET_READS { --sample-id ${meta.id} \ --index ${deacon_idx} \ ${deplete_arg} \ + ${check_pairs_arg} \ --r1-list r1.list \ --r2-list r2.list \ --threads ${task.cpus} \ @@ -270,6 +275,7 @@ process DEACON_ENRICH_SRA_READS { def sracha_threads = Math.max(1, cpus.intdiv(2)) def deacon_threads = Math.max(1, cpus - sracha_threads) def deplete_arg = target_enrichment_enabled ? "" : "--deplete" + def check_pairs_arg = params.check_pairs ? "--check-pairs" : "" """ set -euo pipefail @@ -311,14 +317,14 @@ process DEACON_ENRICH_SRA_READS { exit 1 fi + deacon_args=(filter) case "\${layout}/\${nreads}" in SINGLE/1) read_structure='single' - deacon_inputs=(-) ;; PAIRED/2) read_structure='interleaved' - deacon_inputs=(- -) + deacon_args+=(--interleaved ${check_pairs_arg}) ;; *) printf 'Unsupported sracha read layout for %s: layout=%s nreads=%s\n' \ @@ -337,7 +343,7 @@ process DEACON_ENRICH_SRA_READS { --no-progress \ --yes \ '${run_accession}' \ - | deacon filter \ + | deacon "\${deacon_args[@]}" \ ${deplete_arg} \ --threads ${deacon_threads} \ --abs-threshold ${params.virus_abs_threshold} \ @@ -345,7 +351,7 @@ process DEACON_ENRICH_SRA_READS { --summary '${id}.deacon_filter.json' \ --output '${id}.target_enriched.fastq.gz' \ ${deacon_idx} \ - "\${deacon_inputs[@]}" + - NVD_INPUT_READ_COUNT=\$(jq -er '.seqs_in | select(type == "number" and . >= 0 and . == floor)' '${id}.deacon_filter.json') NVD_ENRICHED_READ_COUNT=\$(jq -er '.seqs_out | select(type == "number" and . >= 0 and . == floor)' '${id}.deacon_filter.json') diff --git a/nextflow.config b/nextflow.config index 9219943b..b0b37d58 100644 --- a/nextflow.config +++ b/nextflow.config @@ -83,6 +83,8 @@ params { // On by default; disable with `--no-merge-pairs` (nvd CLI) or // `--merge_pairs false` (native Nextflow). merge_pairs = true + // whether Deacon should reject mismatched CASAVA or /1 and /2 pair names + check_pairs = false // Skip optional or expensive stages. null keeps direct // Nextflow flag usage ergonomic: `--skip_assembly`, `--skip_blast`, diff --git a/pixi.lock b/pixi.lock index 825c7787..43917c36 100644 --- a/pixi.lock +++ b/pixi.lock @@ -35,7 +35,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/linux-64/alignoth-1.10.1-h9e2517d_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/bbmap-39.91-h09cc210_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/blast-2.17.0-h66d330f_0.conda - - conda: https://conda.anaconda.org/bioconda/linux-64/deacon-0.15.0-hdd79491_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-64/deacon-0.16.0-h3edb6b3_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/entrez-direct-24.0-he881be0_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/linux-64/fastp-1.3.6-h43da1c4_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/flye-2.9.6-py312h734f728_1.conda @@ -348,7 +348,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/linux-aarch64/alignoth-1.10.1-hb267aa8_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/bbmap-39.91-h7683cbd_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/blast-2.16.0-h6a93c2d_5.tar.bz2 - - conda: https://conda.anaconda.org/bioconda/linux-aarch64/deacon-0.15.0-hacc9189_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-aarch64/deacon-0.16.0-ha0850b2_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/entrez-direct-22.4-h8865c2f_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/linux-aarch64/fastp-1.3.6-h5eda1b2_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/flye-2.9.6-py312h01567af_1.conda @@ -669,7 +669,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-64/alignoth-1.10.1-h822ee72_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/bbmap-39.91-h991cdf7_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/blast-2.17.0-h53185c9_0.conda - - conda: https://conda.anaconda.org/bioconda/osx-64/deacon-0.15.0-h5126ea9_0.conda + - conda: https://conda.anaconda.org/bioconda/osx-64/deacon-0.16.0-h80f120d_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/entrez-direct-24.0-h193322a_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/osx-64/fastp-1.3.6-h3119cac_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/flye-2.9.6-py313h95f5313_1.conda @@ -957,7 +957,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-arm64/alignoth-1.10.1-h75efd46_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/bbmap-39.91-he046917_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/blast-2.16.0-hb260f6e_5.tar.bz2 - - conda: https://conda.anaconda.org/bioconda/osx-arm64/deacon-0.15.0-hc0d6d67_0.conda + - conda: https://conda.anaconda.org/bioconda/osx-arm64/deacon-0.16.0-h314a369_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/entrez-direct-22.4-hd5f1084_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/osx-arm64/fastp-1.3.6-ha1d0559_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/flye-2.9.6-py313h05f0e27_1.conda @@ -1245,7 +1245,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/linux-64/alignoth-1.10.1-h9e2517d_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/bbmap-39.91-h09cc210_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/blast-2.17.0-h66d330f_0.conda - - conda: https://conda.anaconda.org/bioconda/linux-64/deacon-0.15.0-hdd79491_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-64/deacon-0.16.0-h3edb6b3_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/entrez-direct-24.0-he881be0_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/linux-64/fastp-1.3.6-h43da1c4_0.conda - conda: https://conda.anaconda.org/bioconda/linux-64/flye-2.9.6-py312h734f728_1.conda @@ -1656,7 +1656,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/linux-aarch64/alignoth-1.10.1-hb267aa8_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/bbmap-39.91-h7683cbd_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/blast-2.16.0-h6a93c2d_5.tar.bz2 - - conda: https://conda.anaconda.org/bioconda/linux-aarch64/deacon-0.15.0-hacc9189_0.conda + - conda: https://conda.anaconda.org/bioconda/linux-aarch64/deacon-0.16.0-ha0850b2_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/entrez-direct-22.4-h8865c2f_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/linux-aarch64/fastp-1.3.6-h5eda1b2_0.conda - conda: https://conda.anaconda.org/bioconda/linux-aarch64/flye-2.9.6-py312h01567af_1.conda @@ -2075,7 +2075,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-64/alignoth-1.10.1-h822ee72_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/bbmap-39.91-h991cdf7_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/blast-2.17.0-h53185c9_0.conda - - conda: https://conda.anaconda.org/bioconda/osx-64/deacon-0.15.0-h5126ea9_0.conda + - conda: https://conda.anaconda.org/bioconda/osx-64/deacon-0.16.0-h80f120d_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/entrez-direct-24.0-h193322a_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/osx-64/fastp-1.3.6-h3119cac_0.conda - conda: https://conda.anaconda.org/bioconda/osx-64/flye-2.9.6-py312hf7bdf34_1.conda @@ -2462,7 +2462,7 @@ environments: - conda: https://conda.anaconda.org/bioconda/osx-arm64/alignoth-1.10.1-h75efd46_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/bbmap-39.91-he046917_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/blast-2.16.0-hb260f6e_5.tar.bz2 - - conda: https://conda.anaconda.org/bioconda/osx-arm64/deacon-0.15.0-hc0d6d67_0.conda + - conda: https://conda.anaconda.org/bioconda/osx-arm64/deacon-0.16.0-h314a369_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/entrez-direct-22.4-hd5f1084_0.tar.bz2 - conda: https://conda.anaconda.org/bioconda/osx-arm64/fastp-1.3.6-ha1d0559_0.conda - conda: https://conda.anaconda.org/bioconda/osx-arm64/flye-2.9.6-py312hf40edc0_1.conda @@ -2897,20 +2897,23 @@ packages: - blast >=2.17.0,<2.18.0a0 size: 84832339 timestamp: 1754909742570 -- conda: https://conda.anaconda.org/bioconda/linux-64/deacon-0.15.0-hdd79491_0.conda - sha256: 69c31e411d270f1311d31fceacb522ddea1b2ba3a04cbfa63e5669e919ecdee6 - md5: 310026878e1a138dfade306cc1ab0950 +- conda: https://conda.anaconda.org/bioconda/linux-64/deacon-0.16.0-h3edb6b3_0.conda + sha256: 05fd4cd0e01d44d746304e6e605af637294d2226adc26ba12bddf30f56d9f1b9 + md5: 816292026eccfd60728456ce183eb146 depends: - - openssl >=3.5.5,<4.0a0 + - __glibc >=2.17,<3.0.a0 + - libgcc >=14 + - openssl >=3.5.7,<4.0a0 constrains: - __glibc >=2.17 license: MIT + license_family: MIT purls: [] run_exports: weak: - - deacon >=0.15.0,<1.0a0 - size: 2160787 - timestamp: 1773419136781 + - deacon >=0.16.0,<0.17.0a0 + size: 2451074 + timestamp: 1786274921399 - conda: https://conda.anaconda.org/bioconda/linux-64/entrez-direct-24.0-he881be0_0.tar.bz2 sha256: 71a8f349659c9c18efa544663de2db1a20b5e3d32f8e4e88cd33110a0caf4eb3 md5: 52a3fabee9201c2c6093c13b4eaf29b4 @@ -3261,20 +3264,22 @@ packages: - blast >=2.16.0,<2.17.0a0 size: 145178710 timestamp: 1743180765609 -- conda: https://conda.anaconda.org/bioconda/linux-aarch64/deacon-0.15.0-hacc9189_0.conda - sha256: 49b7e515bc70b2a5838bdfd75088f3d1e098868935be819939ffa1e1a7a89e3f - md5: 049bd11a7743471d3dbd5bd438f161be +- conda: https://conda.anaconda.org/bioconda/linux-aarch64/deacon-0.16.0-ha0850b2_0.conda + sha256: d9e12ec77eb3b00929ae546c348d1c1bbb6c751a6400ec246c91ada7056e295e + md5: 13cc68a2750d208c2be56e66d86f27c2 depends: - - openssl >=3.5.5,<4.0a0 + - libgcc >=14 + - openssl >=3.5.7,<4.0a0 constrains: - __glibc >=2.17 license: MIT + license_family: MIT purls: [] run_exports: weak: - - deacon >=0.15.0,<1.0a0 - size: 2197786 - timestamp: 1773419019481 + - deacon >=0.16.0,<0.17.0a0 + size: 2448963 + timestamp: 1786274674979 - conda: https://conda.anaconda.org/bioconda/linux-aarch64/entrez-direct-22.4-h8865c2f_0.tar.bz2 sha256: 2ed85538e2e3ae19129faf078afed08b44b32ba12ff3d6982176ccb7879adf75 md5: 289ba7576234e65af7824ef380db651a @@ -3790,20 +3795,22 @@ packages: - blast >=2.17.0,<2.18.0a0 size: 152602258 timestamp: 1754920755551 -- conda: https://conda.anaconda.org/bioconda/osx-64/deacon-0.15.0-h5126ea9_0.conda - sha256: 0ce07edc1f37913b2b83cd75bd9e5c1ddbf8b69b3d6eaf9be9d4b006cab0b9aa - md5: 6e9dcc0f8de108fff7b71fa08572e08e +- conda: https://conda.anaconda.org/bioconda/osx-64/deacon-0.16.0-h80f120d_0.conda + sha256: a5c5e9cffdf74742c1491f8c5f6179218cf3a795133edc118c97e005395e960d + md5: 43dfcdf69a5043fdad90bdf90dc99cfe depends: - - openssl >=3.5.5,<4.0a0 + - __osx >=10.13 + - openssl >=3.5.7,<4.0a0 constrains: - __osx >=10.13 license: MIT + license_family: MIT purls: [] run_exports: weak: - - deacon >=0.15.0,<1.0a0 - size: 1395517 - timestamp: 1773420060364 + - deacon >=0.16.0,<0.17.0a0 + size: 2002368 + timestamp: 1786276123016 - conda: https://conda.anaconda.org/bioconda/osx-64/entrez-direct-24.0-h193322a_0.tar.bz2 sha256: a9257df9c0f41436190e6161fca969308b082b4a6b51e50f27600dfd035637c7 md5: 1a3b2f781a1a9f3a14f26dd9d7660998 @@ -4153,20 +4160,22 @@ packages: - blast >=2.16.0,<2.17.0a0 size: 157885023 timestamp: 1743182400505 -- conda: https://conda.anaconda.org/bioconda/osx-arm64/deacon-0.15.0-hc0d6d67_0.conda - sha256: 4ea72a07956e653e7d17ccd4addf5b02ede55d99abf371288dfd58869176ae65 - md5: 6ec8b3f1a537bcd5043bed62f9352a43 +- conda: https://conda.anaconda.org/bioconda/osx-arm64/deacon-0.16.0-h314a369_0.conda + sha256: a3b1d2a8eea8ff5fc9e572e5065ba1d1a9b0fed420561894445667ebb39e521f + md5: beb02a73e6694cf04bda74d97c648cda depends: - - openssl >=3.5.5,<4.0a0 + - __osx >=11.0 + - openssl >=3.5.7,<4.0a0 constrains: - __osx >=11.0 license: MIT + license_family: MIT purls: [] run_exports: weak: - - deacon >=0.15.0,<1.0a0 - size: 1297184 - timestamp: 1773418777553 + - deacon >=0.16.0,<0.17.0a0 + size: 1780226 + timestamp: 1786274455646 - conda: https://conda.anaconda.org/bioconda/osx-arm64/entrez-direct-22.4-hd5f1084_0.tar.bz2 sha256: 46efea9f61526340ae8db2c1e25666fc46a259f5a52e8d69d5ca62fedf46774d md5: 56c3860d52479d20d86dd78f24947e8e diff --git a/pyproject.toml b/pyproject.toml index 2627dd8e..986e4ed1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,7 +40,7 @@ platforms = ["linux-64", "linux-aarch64", "osx-64", "osx-arm64"] alignoth = ">=1.6.0,<2" bbmap = ">=39.19,<40" blast = ">=2.16.0,<3" -deacon = ">=0.15.0,<0.16" +deacon = ">=0.16.0,<0.17" fastp = ">=1.0.1,<2" fastqc = "==0.12.1" flye = ">=2.9.6,<3" diff --git a/schemas/nvd-params.v3.5.0.schema.json b/schemas/nvd-params.v3.5.0.schema.json index 127bea39..d9118221 100644 --- a/schemas/nvd-params.v3.5.0.schema.json +++ b/schemas/nvd-params.v3.5.0.schema.json @@ -190,6 +190,11 @@ "default": true, "description": "Merge overlapping paired-end reads before contig mapback. On by default; disable with --no-merge-pairs." }, + "check_pairs": { + "type": "boolean", + "default": false, + "description": "Ask Deacon to validate paired FASTQ record names and fail on mismatches. Accepted names use Illumina CASAVA mate fields or /1 and /2 suffixes; validation does not repair pairs." + }, "dedup": { "type": "boolean", "default": false, diff --git a/tests/test_deacon_pairing.py b/tests/test_deacon_pairing.py index e3d948c2..5fb83b12 100644 --- a/tests/test_deacon_pairing.py +++ b/tests/test_deacon_pairing.py @@ -32,6 +32,7 @@ def write_interleaved_fastq( path: Path, *, contaminant_mate: int, + canonical_names: bool = False, ) -> None: """Write one SRA-style interleaved pair with one contaminant mate.""" sequences = { @@ -41,7 +42,10 @@ def write_interleaved_fastq( with gzip.open(path, "wt", encoding="utf-8") as handle: for mate in (1, 2): sequence = sequences[mate] - handle.write(f"@spot.{mate}\n{sequence}\n+\n{'I' * len(sequence)}\n") + separator = "/" if canonical_names else "." + handle.write( + f"@spot{separator}{mate}\n{sequence}\n+\n{'I' * len(sequence)}\n", + ) def build_deacon_index(tmp_path: Path) -> Path: @@ -77,11 +81,17 @@ def run_depletion( *, read_structure: str, contaminant_mate: int = 1, -) -> str: - """Run the public DEACON_DEPLETE process and return its FASTQ text.""" + check_pairs: bool = False, + canonical_names: bool = False, +) -> tuple[str, str]: + """Run DEACON_DEPLETE and return its FASTQ text and rendered task command.""" reads = tmp_path / "reads.fastq.gz" if read_structure == "interleaved": - write_interleaved_fastq(reads, contaminant_mate=contaminant_mate) + write_interleaved_fastq( + reads, + contaminant_mate=contaminant_mate, + canonical_names=canonical_names, + ) else: with gzip.open(reads, "wt", encoding="utf-8") as handle: handle.write( @@ -96,6 +106,7 @@ def run_depletion( params.host_abs_threshold = 1 params.host_rel_threshold = 0.0 +params.check_pairs = {str(check_pairs).lower()} include {{ DEACON_DEPLETE }} from '{DEACON_MODULE}' @@ -130,7 +141,10 @@ def run_depletion( ) assert len(outputs) == 1, outputs with gzip.open(outputs[0], "rt", encoding="utf-8") as handle: - return handle.read() + output = handle.read() + commands = list((tmp_path / "work").glob("**/.command.sh")) + assert len(commands) == 1, commands + return output, commands[0].read_text(encoding="utf-8") @pytest.mark.parametrize("contaminant_mate", [MATE_ONE, MATE_TWO]) @@ -139,7 +153,7 @@ def test_depletion_removes_complete_interleaved_pair_when_either_mate_matches( contaminant_mate: int, ) -> None: """A contaminant match on either mate removes the complete pair.""" - output = run_depletion( + output, _command = run_depletion( tmp_path, read_structure="interleaved", contaminant_mate=contaminant_mate, @@ -150,6 +164,43 @@ def test_depletion_removes_complete_interleaved_pair_when_either_mate_matches( def test_depletion_keeps_direct_file_filtering_for_single_reads(tmp_path: Path) -> None: """Single reads retain Deacon's native direct gzip input path.""" - output = run_depletion(tmp_path, read_structure="single") + output, command = run_depletion(tmp_path, read_structure="single") assert "@single\n" in output + assert "reads.fastq.gz" in command + assert "--interleaved" not in command + assert "--check-pairs" not in command + assert "gzip -dc" not in command + + +def test_interleaved_depletion_uses_explicit_deacon_pair_mode(tmp_path: Path) -> None: + """Interleaved files use Deacon's native compressed-input interface.""" + _output, command = run_depletion(tmp_path, read_structure="interleaved") + + assert "--interleaved" in command + assert "gzip -dc" not in command + assert " - -" not in command + + +def test_interleaved_depletion_forwards_opt_in_pair_validation( + tmp_path: Path, +) -> None: + """NVD only asks Deacon to validate names when explicitly enabled.""" + unchecked_dir = tmp_path / "unchecked" + checked_dir = tmp_path / "checked" + unchecked_dir.mkdir() + checked_dir.mkdir() + + _unchecked_output, unchecked_command = run_depletion( + unchecked_dir, + read_structure="interleaved", + ) + _checked_output, checked_command = run_depletion( + checked_dir, + read_structure="interleaved", + check_pairs=True, + canonical_names=True, + ) + + assert "--check-pairs" not in unchecked_command + assert "--check-pairs" in checked_command diff --git a/tests/test_target_enrichment.py b/tests/test_target_enrichment.py index 3bb005ea..2624d09a 100644 --- a/tests/test_target_enrichment.py +++ b/tests/test_target_enrichment.py @@ -32,6 +32,11 @@ def write_fake_deacon(bin_dir: Path, *, seqs_out: int = 1) -> None: arguments = sys.argv[1:] summary = arguments[arguments.index("--summary") + 1] output = arguments[arguments.index("--output") + 1] +stdin_payload = sys.stdin.buffer.read() +with open("deacon.args", "w") as stream: + stream.write("\\n".join(arguments)) +with open("deacon.stdin", "wb") as stream: + stream.write(stdin_payload) with gzip.open(output, "wt") as stream: if {seqs_out}: stream.write("@retained\\nACGT\\n+\\nIIII\\n") @@ -43,6 +48,57 @@ def write_fake_deacon(bin_dir: Path, *, seqs_out: int = 1) -> None: executable.chmod(0o755) +def write_fake_sracha(bin_dir: Path, *, paired: bool = True) -> None: + """Install an SRA-stream stand-in with deterministic layout metadata.""" + layout = "PAIRED" if paired else "SINGLE" + nreads = 2 if paired else 1 + payload = ( + "@pair/1\nACGT\n+\nIIII\n@pair/2\nTGCA\n+\nIIII\n" + if paired + else "@single\nACGT\n+\nIIII\n" + ) + executable = bin_dir / "sracha" + executable.write_text( + f"""#!/usr/bin/env python3 +import sys + +if sys.argv[1] == "info": + print("accession\\tarchive_type\\tlayout\\tnreads\\tspots\\tsize_bytes\\tplatform\\tmd5") + print("SRR_TEST\\tSRA\\t{layout}\\t{nreads}\\t1\\t100\\tILLUMINA\\tdeadbeef") +elif sys.argv[1] == "get": + sys.stdout.write({payload!r}) +else: + raise SystemExit(2) +""", + encoding="utf-8", + ) + executable.chmod(0o755) + + +def write_fake_stream_helper(bin_dir: Path) -> None: + """Install a deterministic stand-in for the multi-file FIFO helper.""" + executable = bin_dir / "stream_fastqs_to_deacon.py" + executable.write_text( + """#!/usr/bin/env python3 +import gzip +import json +import sys + +arguments = sys.argv[1:] +output = arguments[arguments.index("--output") + 1] +summary = arguments[arguments.index("--summary") + 1] +with open("stream-helper.args", "w") as stream: + stream.write("\\n".join(arguments)) +with gzip.open(output, "wt") as stream: + stream.write("@retained/1\\nACGT\\n+\\nIIII\\n@retained/2\\nTGCA\\n+\\nIIII\\n") +with open(summary, "w") as stream: + json.dump({"seqs_in": 4, "seqs_out": 2, "seqs_removed": 2}, stream) +""", + encoding="utf-8", + ) + executable.chmod(0o755) + + def test_target_enrichment_outputs_are_published(tmp_path: Path) -> None: """Configured target enrichment publishes retained reads and its summary.""" bin_dir = tmp_path / "bin" @@ -129,6 +185,266 @@ def test_target_enrichment_outputs_are_published(tmp_path: Path) -> None: assert (target_results / "summaries" / "sample_A.deacon_filter.json").is_file() +@pytest.mark.parametrize("check_pairs", [False, True]) +def test_paired_target_enrichment_forwards_pair_validation( + tmp_path: Path, + check_pairs: bool, +) -> None: + """The direct two-file path forwards opt-in Deacon name validation.""" + bin_dir = tmp_path / "bin" + bin_dir.mkdir() + write_fake_deacon(bin_dir) + + r1 = tmp_path / "R1.fastq.gz" + r2 = tmp_path / "R2.fastq.gz" + for path, name in ((r1, "pair/1"), (r2, "pair/2")): + with gzip.open(path, "wt", encoding="utf-8") as handle: + handle.write(f"@{name}\nACGT\n+\nIIII\n") + target_index = tmp_path / "target.idx" + target_index.touch() + + workflow = tmp_path / "main.nf" + workflow.write_text( + f"""\ +nextflow.enable.dsl = 2 + +params.check_pairs = {str(check_pairs).lower()} +params.virus_abs_threshold = 1 +params.virus_rel_threshold = 0.0 + +include {{ DEACON_ENRICH_TARGET_READS }} from '{DEACON_MODULE}' + +workflow {{ + DEACON_ENRICH_TARGET_READS(Channel.of(tuple( + [ + id: 'sample_A', + platform: 'illumina', + read_mode: 'paired', + r1_count: 1, + deacon_read_structure: 'interleaved', + ], + [file('{r1}'), file('{r2}')], + file('{target_index}'), + true, + ))) +}} +""", + encoding="utf-8", + ) + + environment = os.environ.copy() + environment["PATH"] = f"{bin_dir}:{environment['PATH']}" + environment["NXF_ANSI_LOG"] = "false" + completed = subprocess.run( # noqa: S603 + ["nextflow", "-C", "/dev/null", "run", str(workflow)], # noqa: S607 + cwd=tmp_path, + env=environment, + text=True, + capture_output=True, + check=False, + timeout=45, + ) + + diagnostics = f"stdout:\n{completed.stdout}\nstderr:\n{completed.stderr}" + assert completed.returncode == 0, diagnostics + argument_files = list((tmp_path / "work").glob("**/deacon.args")) + assert len(argument_files) == 1, argument_files + arguments = argument_files[0].read_text(encoding="utf-8").splitlines() + assert ("--check-pairs" in arguments) is check_pairs + assert "--interleaved" not in arguments + + +@pytest.mark.parametrize("check_pairs", [False, True]) +def test_paired_sra_enrichment_uses_explicit_checked_interleaving( + tmp_path: Path, + check_pairs: bool, +) -> None: + """The paired SRA stream reaches Deacon as one checked interleaved stdin.""" + bin_dir = tmp_path / "bin" + bin_dir.mkdir() + write_fake_deacon(bin_dir) + write_fake_sracha(bin_dir) + target_index = tmp_path / "target.idx" + target_index.touch() + + workflow = tmp_path / "main.nf" + workflow.write_text( + f"""\ +nextflow.enable.dsl = 2 + +params.check_pairs = {str(check_pairs).lower()} +params.max_concurrent_downloads = 1 +params.virus_abs_threshold = 1 +params.virus_rel_threshold = 0.0 + +include {{ DEACON_ENRICH_SRA_READS }} from '{DEACON_MODULE}' + +workflow {{ + DEACON_ENRICH_SRA_READS(Channel.of(tuple( + 'sample_A', + 'illumina', + 'SRR_TEST', + file('{target_index}'), + true, + ))) +}} +""", + encoding="utf-8", + ) + + environment = os.environ.copy() + environment["PATH"] = f"{bin_dir}:{environment['PATH']}" + environment["NXF_ANSI_LOG"] = "false" + completed = subprocess.run( # noqa: S603 + ["nextflow", "-C", "/dev/null", "run", str(workflow)], # noqa: S607 + cwd=tmp_path, + env=environment, + text=True, + capture_output=True, + check=False, + timeout=45, + ) + + diagnostics = f"stdout:\n{completed.stdout}\nstderr:\n{completed.stderr}" + assert completed.returncode == 0, diagnostics + argument_files = list((tmp_path / "work").glob("**/deacon.args")) + assert len(argument_files) == 1, argument_files + arguments = argument_files[0].read_text(encoding="utf-8").splitlines() + assert "--interleaved" in arguments + assert ("--check-pairs" in arguments) is check_pairs + assert arguments.count("-") == 1 + stdin_files = list((tmp_path / "work").glob("**/deacon.stdin")) + assert len(stdin_files) == 1, stdin_files + stdin_payload = stdin_files[0].read_text(encoding="utf-8") + assert "@pair/1\n" in stdin_payload + assert "@pair/2\n" in stdin_payload + + +def test_single_sra_enrichment_omits_pair_arguments(tmp_path: Path) -> None: + """A single-read SRA stream expands the empty pair-argument array safely.""" + bin_dir = tmp_path / "bin" + bin_dir.mkdir() + write_fake_deacon(bin_dir) + write_fake_sracha(bin_dir, paired=False) + target_index = tmp_path / "target.idx" + target_index.touch() + + workflow = tmp_path / "main.nf" + workflow.write_text( + f"""\ +nextflow.enable.dsl = 2 + +params.check_pairs = true +params.max_concurrent_downloads = 1 +params.virus_abs_threshold = 1 +params.virus_rel_threshold = 0.0 + +include {{ DEACON_ENRICH_SRA_READS }} from '{DEACON_MODULE}' + +workflow {{ + DEACON_ENRICH_SRA_READS(Channel.of(tuple( + 'sample_A', 'illumina', 'SRR_TEST', file('{target_index}'), true, + ))) +}} +""", + encoding="utf-8", + ) + + environment = os.environ.copy() + environment["PATH"] = f"{bin_dir}:{environment['PATH']}" + environment["NXF_ANSI_LOG"] = "false" + completed = subprocess.run( # noqa: S603 + ["nextflow", "-C", "/dev/null", "run", str(workflow)], # noqa: S607 + cwd=tmp_path, + env=environment, + text=True, + capture_output=True, + check=False, + timeout=45, + ) + + diagnostics = f"stdout:\n{completed.stdout}\nstderr:\n{completed.stderr}" + assert completed.returncode == 0, diagnostics + [argument_file] = list((tmp_path / "work").glob("**/deacon.args")) + arguments = argument_file.read_text(encoding="utf-8").splitlines() + assert "--interleaved" not in arguments + assert "--check-pairs" not in arguments + assert arguments.count("-") == 1 + [stdin_file] = list((tmp_path / "work").glob("**/deacon.stdin")) + assert "@single\n" in stdin_file.read_text(encoding="utf-8") + + +@pytest.mark.parametrize("check_pairs", [False, True]) +def test_multifile_paired_enrichment_forwards_validation_to_stream_helper( + tmp_path: Path, + check_pairs: bool, +) -> None: + """The multi-lane paired branch carries the opt-in policy to its helper.""" + bin_dir = tmp_path / "bin" + bin_dir.mkdir() + write_fake_stream_helper(bin_dir) + reads = [] + for lane in (1, 2): + for mate in (1, 2): + path = tmp_path / f"L{lane}_R{mate}.fastq.gz" + with gzip.open(path, "wt", encoding="utf-8") as handle: + handle.write(f"@lane{lane}/{mate}\nACGT\n+\nIIII\n") + reads.append(path) + ordered_reads = [reads[0], reads[2], reads[1], reads[3]] + target_index = tmp_path / "target.idx" + target_index.touch() + + workflow = tmp_path / "main.nf" + staged_reads = ", ".join(f"file('{path}')" for path in ordered_reads) + workflow.write_text( + f"""\ +nextflow.enable.dsl = 2 + +params.check_pairs = {str(check_pairs).lower()} +params.virus_abs_threshold = 1 +params.virus_rel_threshold = 0.0 + +include {{ DEACON_ENRICH_TARGET_READS }} from '{DEACON_MODULE}' + +workflow {{ + DEACON_ENRICH_TARGET_READS(Channel.of(tuple( + [ + id: 'sample_A', + platform: 'illumina', + read_mode: 'paired', + r1_count: 2, + deacon_read_structure: 'interleaved', + ], + [{staged_reads}], + file('{target_index}'), + true, + ))) +}} +""", + encoding="utf-8", + ) + + environment = os.environ.copy() + environment["PATH"] = f"{bin_dir}:{environment['PATH']}" + environment["NXF_ANSI_LOG"] = "false" + completed = subprocess.run( # noqa: S603 + ["nextflow", "-C", "/dev/null", "run", str(workflow)], # noqa: S607 + cwd=tmp_path, + env=environment, + text=True, + capture_output=True, + check=False, + timeout=45, + ) + + diagnostics = f"stdout:\n{completed.stdout}\nstderr:\n{completed.stderr}" + assert completed.returncode == 0, diagnostics + argument_files = list((tmp_path / "work").glob("**/stream-helper.args")) + assert len(argument_files) == 1, argument_files + arguments = argument_files[0].read_text(encoding="utf-8").splitlines() + assert ("--check-pairs" in arguments) is check_pairs + + def test_zero_target_reads_stop_before_preprocessing(tmp_path: Path) -> None: """A sample with no retained target reads emits no preprocessing batch.""" bin_dir = tmp_path / "bin" From 41bd4bf5444d5ef2afca5b068458e7dab980a834 Mon Sep 17 00:00:00 2001 From: Nick Minor Date: Sat, 8 Aug 2026 08:52:45 -0500 Subject: [PATCH 08/10] feat: reapply target enrichment after pair merging Apply the configured Deacon target index independently to both BBMerge output classes when target enrichment and pair merging are enabled. This removes nonmatching former mates retained by the pair-atomic initial pass while preserving genuinely interleaved, merge-disabled, and originally single-read paths. Both post-merge classes are refiltered: overlap_merged_pair, and the single_read class BBMerge emits for pairs that did not overlap. Reads that never entered MERGE_PAIRS, meaning anything that is not Illumina interleaved, keep the initial pass untouched, because pair-atomic retention was never a question for them. merge_pairs defaults to true as of v3.4.0, so merging and then refiltering is ordinary behaviour for every run that supplies a target index, not an opt-in experiment. The gate is target enrichment plus merging, so --no-merge-pairs or --no-enrichment turns it off. Emit one Deacon summary per post-merge query class and publish those summaries beneath the existing target-enrichment results at summaries/postmerge. MultiQC reporting is unchanged; the summaries are retained for later reporting decisions. Exercise the semantic premise directly against the locked Deacon CLI: a single FASTQ containing one matching and one nonmatching former mate retains only the match and reports two sequences in and one out. Channel routing remains intentionally reserved for end-to-end validation. --- conf/results.config | 9 ++++++ modules/deacon.nf | 29 ++++++++++++++++++ subworkflows/preprocess_reads.nf | 20 +++++++++++-- tests/test_deacon_pairing.py | 50 ++++++++++++++++++++++++++++++++ 4 files changed, 106 insertions(+), 2 deletions(-) diff --git a/conf/results.config b/conf/results.config index 00b3306b..10902f3a 100644 --- a/conf/results.config +++ b/conf/results.config @@ -126,6 +126,15 @@ process { ] } + withName: 'DEACON_REENRICH_POSTMERGE_READS' { + publishDir = [ + path: { params.target_enrichment_summaries + "/postmerge" }, + mode: 'copy', + pattern: '*.json', + enabled: !params.no_enrichment && (params.virus_index || params.virus_index_url || params.virus_reference_fasta) + ] + } + // ========================================================================= // 02 — Preprocessed reads (final preprocessing output) // ========================================================================= diff --git a/modules/deacon.nf b/modules/deacon.nf index ec693b29..29fea0a0 100644 --- a/modules/deacon.nf +++ b/modules/deacon.nf @@ -253,6 +253,35 @@ process DEACON_ENRICH_TARGET_READS { """ } +process DEACON_REENRICH_POSTMERGE_READS { + /* Reapply target enrichment independently to records emitted by pair merging. */ + + tag "${meta.id}, ${meta.query_class}" + label "medium" + + errorStrategy { task.attempt < 3 ? 'retry' : 'ignore' } + maxRetries 2 + + input: + tuple val(meta), path(reads), path(deacon_idx) + + output: + tuple val(meta), path("${meta.id}.${meta.query_class}.per_read_target_enriched.fastq.gz"), emit: reads + tuple val(meta), path("${meta.id}.${meta.query_class}.per_read_deacon.json"), emit: stats + + script: + """ + deacon filter \ + --threads ${task.cpus} \ + --abs-threshold ${params.virus_abs_threshold} \ + --rel-threshold ${params.virus_rel_threshold} \ + --summary ${meta.id}.${meta.query_class}.per_read_deacon.json \ + --output ${meta.id}.${meta.query_class}.per_read_target_enriched.fastq.gz \ + ${deacon_idx} \ + ${reads} + """ +} + process DEACON_ENRICH_SRA_READS { /* Stream one resolved SRA run through deacon without materializing decoded FASTQ files. */ diff --git a/subworkflows/preprocess_reads.nf b/subworkflows/preprocess_reads.nf index 8e7b5f47..4693a74b 100644 --- a/subworkflows/preprocess_reads.nf +++ b/subworkflows/preprocess_reads.nf @@ -5,6 +5,7 @@ include { DEACON_BUILD_INDEX_FROM_FASTA } from "../modules/dea include { DEACON_UNION_INDEXES } from "../modules/deacon" include { DEACON_ENRICH_TARGET_READS } from "../modules/deacon" include { DEACON_ENRICH_SRA_READS } from "../modules/deacon" +include { DEACON_REENRICH_POSTMERGE_READS } from "../modules/deacon" include { DEACON_DEPLETE } from "../modules/deacon" include { MERGE_PAIRS ; DEDUP_WITH_CLUMPIFY ; TRIM_ADAPTERS ; FILTER_READS } from "../modules/bbmap" include { PROFILE_FASTX as PROFILE_READS ; PLOT_FASTX_LENGTH_PROFILE as PLOT_READ_LENGTH_PROFILES ; PLOT_FASTX_QUALITY_PROFILE as PLOT_READ_QUALITY_PROFILES } from "../modules/fastx" @@ -162,10 +163,25 @@ workflow PREPROCESS_READS { ) } + ch_postmerge_read_batches = ch_merged_read_batches + .mix(ch_unmerged_read_batches) + + // The initial Deacon pass evaluates paired inputs pair-atomically. Once + // BBMerge has converted both output classes to single-read batches, reapply + // the same target index so each resulting FASTQ record must match by itself. + def should_reenrich_postmerge = target_enrichment_enabled && params.merge_pairs + if (should_reenrich_postmerge) { + DEACON_REENRICH_POSTMERGE_READS( + ch_postmerge_read_batches.combine(ch_target_index) + ) + ch_postmerge_batches_for_qc = DEACON_REENRICH_POSTMERGE_READS.out.reads + } else { + ch_postmerge_batches_for_qc = ch_postmerge_read_batches + } + // With pair merging disabled, every sample likewise has one single_read batch. ch_read_batches = params.merge_pairs - ? ch_merged_read_batches - .mix(ch_unmerged_read_batches) + ? ch_postmerge_batches_for_qc .mix(ch_nonmergeable_read_batches) : ch_retained_target_reads.map { sample_id, platform, read_structure, reads -> tuple( diff --git a/tests/test_deacon_pairing.py b/tests/test_deacon_pairing.py index 5fb83b12..5b900fee 100644 --- a/tests/test_deacon_pairing.py +++ b/tests/test_deacon_pairing.py @@ -3,6 +3,7 @@ from __future__ import annotations import gzip +import json import os import shutil import subprocess @@ -204,3 +205,52 @@ def test_interleaved_depletion_forwards_opt_in_pair_validation( assert "--check-pairs" not in unchecked_command assert "--check-pairs" in checked_command + + +def test_single_input_enrichment_filters_records_independently( + tmp_path: Path, +) -> None: + """A single-input target pass can discard one former mate by itself.""" + reads = tmp_path / "postmerge.fastq.gz" + with gzip.open(reads, "wt", encoding="utf-8") as handle: + handle.write(f"@target\n{CONTAMINANT}\n+\n{'I' * len(CONTAMINANT)}\n") + handle.write( + f"@nontarget\n{NONCONTAMINANT}\n+\n{'I' * len(NONCONTAMINANT)}\n", + ) + + index = build_deacon_index(tmp_path) + output = tmp_path / "enriched.fastq.gz" + summary = tmp_path / "summary.json" + completed = subprocess.run( # noqa: S603 + [ + "deacon", + "filter", + "--threads", + "1", + "--abs-threshold", + "1", + "--rel-threshold", + "0.0", + "--summary", + str(summary), + "--output", + str(output), + str(index), + str(reads), + ], # noqa: S607 + check=False, + capture_output=True, + text=True, + timeout=45, + ) + + diagnostics = f"stdout:\n{completed.stdout}\nstderr:\n{completed.stderr}" + assert completed.returncode == 0, diagnostics + with gzip.open(output, "rt", encoding="utf-8") as handle: + enriched_reads = handle.read() + assert "@target\n" in enriched_reads + assert "@nontarget\n" not in enriched_reads + + summary_data = json.loads(summary.read_text(encoding="utf-8")) + assert summary_data["seqs_in"] == 2 + assert summary_data["seqs_out"] == 1 From 54b37cf5e130d26ed47ad448c3291717ca6598ee Mon Sep 17 00:00:00 2001 From: wkgardner Date: Wed, 9 Sep 2026 14:03:38 -0500 Subject: [PATCH 09/10] fix: restore the TestNvdParamsDefaults class boundary "feat!: remove sourmash rapid screening" deleted test_sourmash_reference_params along with the six reference and threshold params it exercised, and took the following `class TestNvdParamsDefaults:` header with it. The 24 default-value tests below that point silently became TestNvdParamsToNextflowArgs methods. Nothing stopped running, so nothing failed and nothing warned. pytest collects the same 110 tests in lib/py_nvd/test_models.py before and after this change; they are only reparented back to the class whose docstring describes them. Reinstate the header immediately before test_default_cutoff_percent, which is where it sat before the removal. --- lib/py_nvd/test_models.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lib/py_nvd/test_models.py b/lib/py_nvd/test_models.py index a06abb47..60b88c9d 100644 --- a/lib/py_nvd/test_models.py +++ b/lib/py_nvd/test_models.py @@ -493,6 +493,10 @@ def test_check_pairs_param(self) -> None: assert cmd[check_pairs_idx + 1] == "true" assert "--check-pairs" not in cmd + +class TestNvdParamsDefaults: + """Tests for default values matching nextflow.config.""" + def test_default_cutoff_percent(self) -> None: """Default cutoff_percent matches nextflow.config.""" assert NvdParams().cutoff_percent == 0.001 From 411f7527540fecac2a4f423b83d956888c498373 Mon Sep 17 00:00:00 2001 From: wkgardner Date: Wed, 9 Sep 2026 14:04:12 -0500 Subject: [PATCH 10/10] docs: describe the v3.5.0 enrichment params as target enrichment "feat: retain BLAST hits across taxa" reframed the virus_* parameters as generic target-enrichment controls in models.py, deliberately leaving the 3.x parameter names alone. schemas/nvd-params.v3.5.0.schema.json did not exist at that point, so it inherited the v3.4.0 wording and still described a prebuilt "vertebrate-infecting virus deacon index" driving "virus read enrichment". That is the text IDEs show on hover, so the published contract described a narrower feature than the pipeline implements. Rewords eight descriptions to match the model. The schema's longer phrasing is kept deliberately: these descriptions are hover documentation and say more than the terse CLI help strings, so this is not a blanket sync of the two. schemas/nvd-params.v3.4.0.schema.json is left alone. It is the published contract for a released version whose enrichment really was described that way. --- schemas/nvd-params.v3.5.0.schema.json | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/schemas/nvd-params.v3.5.0.schema.json b/schemas/nvd-params.v3.5.0.schema.json index d9118221..28d20c7c 100644 --- a/schemas/nvd-params.v3.5.0.schema.json +++ b/schemas/nvd-params.v3.5.0.schema.json @@ -103,7 +103,7 @@ }, "virus_index_version": { "type": "string", - "description": "Virus enrichment index version" + "description": "Target-enrichment index version" }, "nvd_files": { "type": "string", @@ -125,7 +125,7 @@ "null" ], "default": null, - "description": "Path to a prebuilt vertebrate-infecting virus deacon index (.idx file). Providing this enables virus read enrichment." + "description": "Path to a prebuilt Deacon target-enrichment index (.idx file). Providing this enables target enrichment." }, "virus_index_url": { "type": [ @@ -133,7 +133,7 @@ "null" ], "default": null, - "description": "URL to download a prebuilt vertebrate-infecting virus deacon index. Providing this enables virus read enrichment when virus_index is not set." + "description": "URL to download a prebuilt Deacon target-enrichment index. Providing this enables target enrichment when virus_index is not set." }, "virus_reference_fasta": { "type": [ @@ -141,7 +141,7 @@ "null" ], "default": null, - "description": "Custom vertebrate-infecting virus FASTA to build and union with other virus indexes. Providing this enables virus read enrichment." + "description": "Custom target FASTA to build and union with other target-enrichment indexes. Providing this enables target enrichment." }, "no_enrichment": { "type": "boolean", @@ -152,26 +152,26 @@ "type": "integer", "default": 31, "minimum": 1, - "description": "K-mer size for building a custom virus enrichment index" + "description": "K-mer size for building a custom target-enrichment index" }, "virus_window_size": { "type": "integer", "default": 1, "minimum": 1, - "description": "Minimizer window size for building a custom virus enrichment index" + "description": "Minimizer window size for building a custom target-enrichment index" }, "virus_abs_threshold": { "type": "integer", "default": 1, "minimum": 1, - "description": "Minimum absolute minimizer hits for virus read enrichment" + "description": "Minimum absolute minimizer hits for target enrichment" }, "virus_rel_threshold": { "type": "number", "default": 0.0, "minimum": 0, "maximum": 1, - "description": "Minimum relative proportion of minimizers for virus read enrichment (0.0-1.0)" + "description": "Minimum relative proportion of minimizers for target enrichment (0.0-1.0)" }, "sourmash_ksize": { "type": "integer",