perf(mapping): pipe minimap2 straight into samtools sort (removes ~1.3 GiB of disk traffic per mapping)

src/plassembler/utils/external_tools.py

@@ -2,6 +2,7 @@
 import shlex
 import subprocess
 import sys
+from contextlib import ExitStack
 from pathlib import Path
 from typing import List, Optional, Tuple
 
@@ -72,6 +73,69 @@ def run_to_stdout(
     def _run_core(command: List[str], stdout_fh, stderr_fh) -> None:
         subprocess.check_call(command, stdout=stdout_fh, stderr=stderr_fh)
 
+    @staticmethod
+    def run_piped(
+        tools: Tuple["ExternalTool", ...], outfile: Optional[Path] = None
+    ) -> None:
+        """Run several tools as one pipeline, without a shell.
+
+        Used to avoid materialising intermediates that exist only so the next
+        command can read them back - most importantly the full uncompressed SAM
+        that each mapping step used to write to disk before samtools re-read it.
+
+        :param tools: stages, in order; each stage's stdout feeds the next
+        :param outfile: file for the last stage's stdout. If None the last stage
+            writes its own output (e.g. `samtools sort -o`) and its stdout goes
+            to that tool's .out log, matching ExternalTool.run.
+        """
+        joined = " | ".join(tool.command_as_str for tool in tools)
+        logger.info(f"Started running {joined} ...")
+
+        procs: List[subprocess.Popen] = []
+        with ExitStack() as stack:
+            if outfile is None:
+                last_stdout = stack.enter_context(open(tools[-1].out_log, "w"))
+            else:
+                last_stdout = stack.enter_context(open(outfile, "wb"))
+
+            upstream_stdout = None
+            try:
+                for index, tool in enumerate(tools):
+                    stderr_fh = stack.enter_context(open(tool.err_log, "w"))
+                    print(f"Command line: {tool.command_as_str}", file=stderr_fh)
+                    is_last = index == len(tools) - 1
+                    proc = subprocess.Popen(
+                        tool.command,
+                        stdin=upstream_stdout,
+                        stdout=last_stdout if is_last else subprocess.PIPE,
+                        stderr=stderr_fh,
+                    )
+                    # the parent must drop its copy of the upstream read end,
+                    # otherwise that stage never sees EOF
+                    if upstream_stdout is not None:
+                        upstream_stdout.close()
+                    upstream_stdout = proc.stdout
+                    procs.append(proc)
+            except OSError:
+                # a stage failed to start (missing binary, fork failure). Do not
+                # leave the earlier stages of a half-built chain running.
+                for proc in procs:
+                    proc.kill()
+                    proc.wait()
+                raise
+
+            # wait downstream-first: an upstream stage blocked writing into a
+            # dead pipe only gets its SIGPIPE once the reader is gone
+            for proc in reversed(procs):
+                proc.wait()
+
+        # report the earliest failing stage, which is the informative one
+        for tool, proc in zip(tools, procs):
+            if proc.returncode != 0:
+                raise subprocess.CalledProcessError(proc.returncode, tool.command)
+
+        logger.info(f"Done running {joined}")
+
     @staticmethod
     def run_tools(
         tools_to_run: Tuple["ExternalTool", ...], ctx: Optional[click.Context] = None

src/plassembler/utils/mapping.py

@@ -1,10 +1,93 @@
+from pathlib import Path
+
 from plassembler.utils.external_tools import ExternalTool
 
 #################################
 # original mapping
 #################################
 
 
+def minimap2_model_for(pacbio_model):
+    """maps plassembler's --pacbio_model onto a minimap2 preset
+    :param pacbio_model: pacbio_model
+    :return: minimap2 preset name
+    """
+    if pacbio_model in ("--pacbio-raw", "--pacbio-corr"):
+        return "map-pb"
+    if pacbio_model == "--pacbio-hifi":
+        return "map-hifi"
+    # ONT, and the "nothing" default
+    return "map-ont"
+
+
+def _minimap2_tool(params, logdir, outfile=""):
+    return ExternalTool(
+        tool="minimap2",
+        input="",
+        output="",
+        params=params,
+        logdir=logdir,
+        outfile=outfile,
+    )
+
+
+def _samtools_sort_tool(sorted_bam, threads, logdir):
+    return ExternalTool(
+        tool="samtools",
+        input="",
+        output="",
+        params=f" sort -@ {threads} -o {sorted_bam}",
+        logdir=logdir,
+        outfile="",
+    )
+
+
+def minimap_long_reads_to_sorted_bam(
+    input_long_reads, fasta, sorted_bam: Path, threads, pacbio_model, logdir
+):
+    """maps long reads with minimap2 straight into a sorted bam
+
+    minimap2's SAM was previously written to disk in full and then re-read by
+    samtools sort. For an ONT isolate that is ~0.6 GiB written and read back per
+    mapping, twice per run, purely as a pipe buffer.
+
+    :param input_long_reads: reads to map
+    :param fasta: reference
+    :param sorted_bam: output sorted bam
+    :param threads: threads
+    :param pacbio_model: pacbio_model
+    :param logdir: logdir
+    :return:
+    """
+    minimap2_model = minimap2_model_for(pacbio_model)
+    ExternalTool.run_piped(
+        (
+            _minimap2_tool(
+                f" -ax {minimap2_model} -t {threads} {fasta} {input_long_reads}", logdir
+            ),
+            _samtools_sort_tool(sorted_bam, threads, logdir),
+        )
+    )
+
+
+def minimap_short_reads_to_sorted_bam(r1, r2, fasta, sorted_bam: Path, threads, logdir):
+    """maps short reads with minimap2 straight into a sorted bam
+    :param r1: R1 reads
+    :param r2: R2 reads
+    :param fasta: reference
+    :param sorted_bam: output sorted bam
+    :param threads: threads
+    :param logdir: logdir
+    :return:
+    """
+    ExternalTool.run_piped(
+        (
+            _minimap2_tool(f" -ax sr -t {threads} {fasta} {r1} {r2}", logdir),
+            _samtools_sort_tool(sorted_bam, threads, logdir),
+        )
+    )
+
+
 def minimap_long_reads(input_long_reads, fasta, sam, threads, pacbio_model, logdir):
     """maps long reads using minimap2
     :param threads: threads
@@ -14,25 +97,11 @@ def minimap_long_reads(input_long_reads, fasta, sam, threads, pacbio_model, logd
     :return:
     """
 
-    # ONT
-    minimap2_model = "map-ont"
-
-    # Pacbio
-    if pacbio_model == "nothing":
-        minimap2_model = "map-ont"
-    elif pacbio_model == "--pacbio-raw":
-        minimap2_model = "map-pb"
-    elif pacbio_model == "--pacbio-corr":
-        minimap2_model = "map-pb"
-    elif pacbio_model == "--pacbio-hifi":
-        minimap2_model = "map-hifi"
+    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,
+    minimap2 = _minimap2_tool(
+        f" -ax {minimap2_model} -t {threads} {fasta} {input_long_reads}",
+        logdir,
         outfile=sam,
     )
 

src/plassembler/utils/plass_class.py

@@ -6,7 +6,6 @@
 from Bio.SeqRecord import SeqRecord
 from loguru import logger
 
-from plassembler.utils.bam import sam_to_sorted_bam
 from plassembler.utils.depth import (
     collate_depths,
     combine_depth_dfs,
@@ -16,7 +15,10 @@
     get_contig_lengths,
     get_depths_from_bam,
 )
-from plassembler.utils.mapping import minimap_long_reads, minimap_short_reads
+from plassembler.utils.mapping import (
+    minimap_long_reads_to_sorted_bam,
+    minimap_short_reads_to_sorted_bam,
+)
 from plassembler.utils.run_mash import get_contig_count, is_file_empty
 
 
@@ -406,27 +408,20 @@ def get_depth(self, logdir, pacbio_model, threads):
 
         input_long_reads: Path = Path(outdir) / "chopper_long_reads.fastq.gz"
         fasta: Path = Path(outdir) / "combined.fasta"
-        sam_file: Path = Path(outdir) / "combined_long.sam"
         sorted_bam: Path = Path(outdir) / "combined_sorted_long.bam"
 
-        # map
-        minimap_long_reads(
-            input_long_reads, fasta, sam_file, threads, pacbio_model, logdir
+        # map and sort in one pipeline - no intermediate SAM on disk
+        minimap_long_reads_to_sorted_bam(
+            input_long_reads, fasta, sorted_bam, threads, pacbio_model, logdir
         )
-        # sort
-        sam_to_sorted_bam(sam_file, sorted_bam, threads, logdir)
 
         # short reads
         r1: Path = Path(outdir) / "trimmed_R1.fastq"
         r2: Path = Path(outdir) / "trimmed_R2.fastq"
         fasta: Path = Path(outdir) / "combined.fasta"
-        sam_file: Path = Path(outdir) / "combined_short.sam"
         sorted_bam: Path = Path(outdir) / "combined_sorted_short.bam"
 
-        # map
-        minimap_short_reads(r1, r2, fasta, sam_file, threads, logdir)
-        # sort
-        sam_to_sorted_bam(sam_file, sorted_bam, threads, logdir)
+        minimap_short_reads_to_sorted_bam(r1, r2, fasta, sorted_bam, threads, logdir)
 
         # get contig lengths
 
@@ -463,7 +458,6 @@ def get_depth_long(self, logdir, pacbio_model, threads, plas_fasta):
         input_long_reads: Path = Path(outdir) / "chopper_long_reads.fastq.gz"
         chromosome: Path = Path(outdir) / "chromosome.fasta"
         combined_fasta: Path = Path(outdir) / "long_combined.fasta"
-        sam_file: Path = Path(outdir) / "combined_long.sam"
         sorted_bam: Path = Path(outdir) / "combined_sorted_long.bam"
 
         # # write to combined fasta
@@ -482,12 +476,10 @@ def get_depth_long(self, logdir, pacbio_model, threads, plas_fasta):
         # Write the combined sequences to the output file
         SeqIO.write(combined_sequences, combined_fasta, "fasta")
 
-        # map
-        minimap_long_reads(
-            input_long_reads, combined_fasta, sam_file, threads, pacbio_model, logdir
+        # map and sort in one pipeline - no intermediate SAM on disk
+        minimap_long_reads_to_sorted_bam(
+            input_long_reads, combined_fasta, sorted_bam, threads, pacbio_model, logdir
         )
-        # sort
-        sam_to_sorted_bam(sam_file, sorted_bam, threads, logdir)
 
         # get contig lengths
         contig_lengths = get_contig_lengths(combined_fasta)
@@ -970,29 +962,24 @@ def get_depth(self, logdir, threads, pacbio_model):
 
         input_long_reads: Path = Path(outdir) / "chopper_long_reads.fastq.gz"
         fasta: Path = Path(outdir) / "combined.fasta"
-        sam_file: Path = Path(outdir) / "combined_long.sam"
         sorted_bam: Path = Path(outdir) / "combined_sorted_long.bam"
 
-        # map
+        # map and sort in one pipeline - no intermediate SAM on disk
         if self.long_flag is True:
-            minimap_long_reads(
-                input_long_reads, fasta, sam_file, threads, pacbio_model, logdir
+            minimap_long_reads_to_sorted_bam(
+                input_long_reads, fasta, sorted_bam, threads, pacbio_model, logdir
             )
-            # sort
-            sam_to_sorted_bam(sam_file, sorted_bam, threads, logdir)
 
         # short reads
         r1: Path = Path(outdir) / "trimmed_R1.fastq"
         r2: Path = Path(outdir) / "trimmed_R2.fastq"
         fasta: Path = Path(outdir) / "combined.fasta"
-        sam_file: Path = Path(outdir) / "combined_short.sam"
         sorted_bam: Path = Path(outdir) / "combined_sorted_short.bam"
 
-        # map
         if self.short_flag is True:
-            minimap_short_reads(r1, r2, fasta, sam_file, threads, logdir)
-            # sort
-            sam_to_sorted_bam(sam_file, sorted_bam, threads, logdir)
+            minimap_short_reads_to_sorted_bam(
+                r1, r2, fasta, sorted_bam, threads, logdir
+            )
 
         # get contig lengths
 

tests/test_piped_tools.py

@@ -0,0 +1,83 @@
+"""Tests for ExternalTool.run_piped, used to map straight into a sorted bam."""
+
+import subprocess
+from pathlib import Path
+
+import pytest
+
+from src.plassembler.utils.external_tools import ExternalTool
+
+
+def tool(cmd, params, logdir, outfile=""):
+    return ExternalTool(
+        tool=cmd, input="", output="", params=params, logdir=logdir, outfile=outfile
+    )
+
+
+def test_run_piped_writes_last_stage_stdout(tmp_path):
+    """Stages are connected by pipes and the last one's stdout lands in outfile."""
+    logdir = tmp_path / "logs"
+    out = tmp_path / "out.txt"
+    ExternalTool.run_piped(
+        (
+            tool("printf", r" 'b\na\nc\n'", logdir),
+            tool("sort", "", logdir),
+            tool("tr", " a-z A-Z", logdir),
+        ),
+        outfile=out,
+    )
+    assert out.read_text() == "A\nB\nC\n"
+
+
+def test_run_piped_last_stage_writes_its_own_file(tmp_path):
+    """With outfile=None the last stage writes its own output, as `samtools sort
+    -o` does; its stdout goes to that tool's .out log."""
+    logdir = tmp_path / "logs"
+    target = tmp_path / "written_by_tee.txt"
+    ExternalTool.run_piped(
+        (tool("printf", " 'hello\\n'", logdir), tool("tee", f" {target}", logdir))
+    )
+    assert target.read_text() == "hello\n"
+
+
+def test_run_piped_raises_on_failing_stage(tmp_path):
+    """A non-zero exit anywhere in the chain must surface, not be swallowed the
+    way an unwaited-on process would be."""
+    logdir = tmp_path / "logs"
+    with pytest.raises(subprocess.CalledProcessError):
+        ExternalTool.run_piped(
+            (
+                tool("printf", " 'x\\n'", logdir),
+                tool("cat", " /nonexistent/path/xyz", logdir),
+            ),
+            outfile=tmp_path / "out.txt",
+        )
+
+
+def test_run_piped_reports_the_earliest_failing_stage(tmp_path):
+    """When an upstream stage fails, its return code is the informative one."""
+    logdir = tmp_path / "logs"
+    with pytest.raises(subprocess.CalledProcessError) as excinfo:
+        ExternalTool.run_piped(
+            (
+                tool("cat", " /nonexistent/path/xyz", logdir),
+                tool("cat", "", logdir),
+            ),
+            outfile=tmp_path / "out.txt",
+        )
+    assert Path(excinfo.value.cmd[0]).name == "cat"
+    assert excinfo.value.returncode != 0
+
+
+def test_run_piped_does_not_leave_a_half_built_chain_running(tmp_path):
+    """If a later stage's binary is missing, the stages already started must be
+    killed rather than left waiting on a pipe nobody will read."""
+    logdir = tmp_path / "logs"
+    with pytest.raises(OSError):
+        ExternalTool.run_piped(
+            (
+                tool("cat", "", logdir),
+                tool("definitely_not_a_real_binary_xyz", "", logdir),
+            ),
+            outfile=tmp_path / "out.txt",
+        )