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", + )