perf(long): pipe minimap2 straight into the plasmid FASTQ extraction (removes the last SAM in long mode)

src/plassembler/__init__.py

@@ -55,6 +55,7 @@
 from plassembler.utils.sam_to_fastq import (
     extract_long_fastqs_fast,
     extract_long_fastqs_slow_keep_fastqs,
+    map_and_extract_long_fastqs,
 )
 from plassembler.utils.test_incompatibility import incompatbility
 from plassembler.utils.util import get_version, print_citation
@@ -1562,19 +1563,15 @@ def long(
         else:
             coverage = 50
 
-        logger.info("Mapping long reads.")
+        logger.info("Mapping long reads and extracting plasmid Fastqs.")
         input_long_reads: Path = Path(outdir) / "chopper_long_reads.fastq.gz"
         fasta: Path = Path(outdir) / "flye_renamed.fasta"
-        samfile: Path = Path(outdir) / "long_read.sam"
-        minimap_long_reads(
-            input_long_reads, fasta, samfile, threads, pacbio_model, logdir
-        )
-
-        # for long, custom function is quick enough
-        logger.info("Processing Sam/Bam Files and extracting Fastqs.")
-        samfile: Path = Path(outdir) / "long_read.sam"
         plasmidfastqs: Path = Path(outdir) / "plasmid_long.fastq"
-        extract_long_fastqs_fast(samfile, plasmidfastqs, threads)
+        # long-only mode is the one case where the sam has a single consumer, so
+        # the mapping feeds the extraction directly instead of via a file
+        map_and_extract_long_fastqs(
+            input_long_reads, fasta, plasmidfastqs, threads, pacbio_model, logdir
+        )
 
         # map the validated --pacbio_model to canu's read-type flag, error rate,
         # and whether read correction should be skipped (HiFi reads)

src/plassembler/utils/sam_to_fastq.py

@@ -5,6 +5,9 @@
 
 import pysam
 
+from plassembler.utils.external_tools import ExternalTool
+from plassembler.utils.mapping import minimap2_model_for
+
 
 def extract_long_fastqs_slow_keep_fastqs(out_dir, samname, plasmidname):
     #################################################
@@ -132,8 +135,62 @@ def extract_long_fastqs_slow_keep_fastqs(out_dir, samname, plasmidname):
 """
 
 
+# keep reads whose primary alignment (flag 0 or 16) is to a contig whose name
+# contains "plas", plus every unmapped read (flag 4), and emit them as fastq
+PLASMID_READ_AWK = (
+    '{if((($3 ~ /plas/)&& ($2 == "0"|| $2 == "16"))||($2 == "4"))'
+    ' print "@"$1"\\n"$10"\\n+"$1"\\n"$11}'
+)
+
+
 def extract_long_fastqs_fast(sam_name, plasmidfile, threads):
-    cmd = f'samtools  view -@ {threads} {sam_name} | awk \'{{if((($3 ~ /plas/)&& ($2 == "0"|| $2 == "16"))||($2 == "4")) print "@"$1"\\n"$10"\\n+"$1"\\n"$11}}\' > {plasmidfile}'
+    cmd = f"samtools  view -@ {threads} {sam_name} | awk '{PLASMID_READ_AWK}' > {plasmidfile}"
     # shell=True is required for the samtools | awk pipeline; check=True surfaces
     # a failing samtools instead of silently leaving an empty plasmid FASTQ.
     sp.run(cmd, shell=True, check=True)
+
+
+def map_and_extract_long_fastqs(
+    input_long_reads, fasta, plasmidfile, threads, pacbio_model, logdir
+):
+    """Map long reads and pull the plasmid/unmapped ones straight out of the stream.
+
+    In long-only mode `long_read.sam` has exactly one consumer - this extraction
+    - so writing minimap2's full uncompressed SAM to disk only to have samtools
+    read it straight back is pure I/O. Piping the three stages together removes
+    it entirely: on a real ONT isolate that is ~0.6 GiB written and re-read.
+
+    :param input_long_reads: reads to map
+    :param fasta: reference (renamed flye assembly)
+    :param plasmidfile: output fastq of plasmid and unmapped reads
+    :param threads: threads
+    :param pacbio_model: pacbio_model
+    :param logdir: logdir
+    :return:
+    """
+    minimap2_model = minimap2_model_for(pacbio_model)
+    minimap2 = ExternalTool(
+        tool="minimap2",
+        input="",
+        output="",
+        params=f" -ax {minimap2_model} -t {threads} {fasta} {input_long_reads}",
+        logdir=logdir,
+        outfile="",
+    )
+    samtools_view = ExternalTool(
+        tool="samtools",
+        input="",
+        output="",
+        params=f" view -@ {threads}",
+        logdir=logdir,
+        outfile="",
+    )
+    awk = ExternalTool(
+        tool="awk",
+        input="",
+        output="",
+        params=f" '{PLASMID_READ_AWK}'",
+        logdir=logdir,
+        outfile="",
+    )
+    ExternalTool.run_piped((minimap2, samtools_view, awk), outfile=plasmidfile)

tests/test_long_extraction.py

@@ -0,0 +1,89 @@
+"""Tests for mapping long reads straight into the plasmid FASTQ."""
+
+import hashlib
+import shutil
+from pathlib import Path
+
+import pytest
+
+from src.plassembler.utils.mapping import minimap_long_reads
+from src.plassembler.utils.sam_to_fastq import (
+    PLASMID_READ_AWK,
+    extract_long_fastqs_fast,
+    map_and_extract_long_fastqs,
+)
+
+TEST_DATA = Path("tests/test_data")
+TOOLS_PRESENT = all(shutil.which(tool) for tool in ("minimap2", "samtools", "awk"))
+
+
+def test_awk_program_is_the_documented_filter():
+    """Primary plasmid alignments (flag 0/16) plus every unmapped read (flag 4)."""
+    assert "$3 ~ /plas/" in PLASMID_READ_AWK
+    assert '$2 == "0"' in PLASMID_READ_AWK
+    assert '$2 == "16"' in PLASMID_READ_AWK
+    assert '$2 == "4"' in PLASMID_READ_AWK
+
+
+@pytest.fixture
+def reference(tmp_path):
+    """A chromosome and a plasmid, named the way identify_chromosome_process_*
+    names them - the awk filter keys off "plas" in the contig name."""
+    fasta = tmp_path / "flye_renamed.fasta"
+    fasta.write_text(
+        ">chromosome\n" + "ACGTTGCA" * 400 + "\n>plasmid_1\n" + "GGCCTTAA" * 200 + "\n"
+    )
+    return fasta
+
+
+@pytest.mark.requires_tool
+@pytest.mark.skipif(not TOOLS_PRESENT, reason="minimap2/samtools/awk not installed")
+def test_piped_extraction_matches_the_two_step_form(tmp_path, reference):
+    """The piped pipeline must produce exactly the FASTQ that writing a SAM and
+    running samtools|awk over it produced."""
+    reads = str(TEST_DATA / "test_long.fastq.gz")
+    logdir = tmp_path / "logs"
+
+    two_step = tmp_path / "two_step.fastq"
+    sam = tmp_path / "long_read.sam"
+    minimap_long_reads(reads, reference, sam, "2", "nothing", logdir)
+    extract_long_fastqs_fast(sam, two_step, "2")
+
+    piped = tmp_path / "piped.fastq"
+    map_and_extract_long_fastqs(reads, reference, piped, "2", "nothing", logdir)
+
+    assert (
+        hashlib.sha256(piped.read_bytes()).hexdigest()
+        == hashlib.sha256(two_step.read_bytes()).hexdigest()
+    )
+
+
+@pytest.mark.requires_tool
+@pytest.mark.skipif(not TOOLS_PRESENT, reason="minimap2/samtools/awk not installed")
+def test_piped_extraction_writes_no_intermediate_sam(tmp_path, reference):
+    """The point of the change: nothing lands on disk between the stages."""
+    reads = str(TEST_DATA / "test_long.fastq.gz")
+    out = tmp_path / "plasmid_long.fastq"
+    map_and_extract_long_fastqs(
+        reads, reference, out, "2", "nothing", tmp_path / "logs"
+    )
+
+    assert out.exists()
+    assert list(tmp_path.glob("*.sam")) == []
+
+
+@pytest.mark.requires_tool
+@pytest.mark.skipif(not TOOLS_PRESENT, reason="minimap2/samtools/awk not installed")
+def test_piped_extraction_output_is_valid_fastq(tmp_path, reference):
+    reads = str(TEST_DATA / "test_long.fastq.gz")
+    out = tmp_path / "plasmid_long.fastq"
+    map_and_extract_long_fastqs(
+        reads, reference, out, "2", "nothing", tmp_path / "logs"
+    )
+
+    lines = out.read_text().splitlines()
+    assert len(lines) % 4 == 0
+    for i in range(0, len(lines), 4):
+        assert lines[i].startswith("@")
+        assert lines[i + 2].startswith("+")
+        assert len(lines[i + 1]) == len(lines[i + 3])