From d76ab2c4c0baec706f277b90aff3dfec22566f0d Mon Sep 17 00:00:00 2001 From: Sanjay Nagi Date: Thu, 13 Aug 2026 21:19:38 +0000 Subject: [PATCH] perf(mapping): pipe minimap2 straight into samtools sort Every depth calculation wrote minimap2's full uncompressed SAM to disk and then had samtools sort read it straight back. The SAM had exactly one consumer; it existed only as a pipe buffer. For one real ONT isolate that is 0.64 GiB written and re-read per mapping, and `plassembler run` does this twice (long and short). Nextflow traces of a 90-isolate plate show plassembler doing ~2.2-2.7 GB of writes per sample, a large part of which is this. Add ExternalTool.run_piped, which runs a chain of tools connected by pipes without a shell: it closes the parent's copy of each upstream read end, waits on every stage downstream-first, and raises CalledProcessError for the earliest failing stage. Then add minimap_long_reads_to_sorted_bam / minimap_short_reads_to_sorted_bam on top of it and use them at the six map-then-sort sites in plass_class. The resulting bam is identical: alignment records compare byte-for-byte and `samtools depth` output has the same md5. Only the @PG header line differs, because it records the command line that produced the file. On a fast local NVMe with a warm page cache the wall-clock difference is small (20.4s -> 20.0s); the win is the 1.27 GiB of avoided disk traffic per mapping, which matters when many isolates run concurrently and the page cache is contended. minimap_long_reads / minimap_short_reads are kept: the hybrid path's long_read.sam has two consumers and is handled separately. If a stage fails to start partway through building the chain - a missing binary, a fork failure - the stages already started are killed and reaped rather than left blocked on a pipe nobody will read. --- src/plassembler/utils/external_tools.py | 64 +++++++++++++++ src/plassembler/utils/mapping.py | 105 ++++++++++++++++++++---- src/plassembler/utils/plass_class.py | 47 ++++------- tests/test_piped_tools.py | 83 +++++++++++++++++++ 4 files changed, 251 insertions(+), 48 deletions(-) create mode 100644 tests/test_piped_tools.py diff --git a/src/plassembler/utils/external_tools.py b/src/plassembler/utils/external_tools.py index 05d413f..8c785dd 100644 --- a/src/plassembler/utils/external_tools.py +++ b/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 diff --git a/src/plassembler/utils/mapping.py b/src/plassembler/utils/mapping.py index ad800af..744c37a 100644 --- a/src/plassembler/utils/mapping.py +++ b/src/plassembler/utils/mapping.py @@ -1,3 +1,5 @@ +from pathlib import Path + from plassembler.utils.external_tools import ExternalTool ################################# @@ -5,6 +7,87 @@ ################################# +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, ) diff --git a/src/plassembler/utils/plass_class.py b/src/plassembler/utils/plass_class.py index b8e7318..97795af 100644 --- a/src/plassembler/utils/plass_class.py +++ b/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 diff --git a/tests/test_piped_tools.py b/tests/test_piped_tools.py new file mode 100644 index 0000000..7c5bfab --- /dev/null +++ b/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", + )