diff --git a/src/plassembler/__init__.py b/src/plassembler/__init__.py index 5bc0476..bc70af0 100644 --- a/src/plassembler/__init__.py +++ b/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) diff --git a/src/plassembler/utils/sam_to_fastq.py b/src/plassembler/utils/sam_to_fastq.py index 00ec25b..1049b81 100644 --- a/src/plassembler/utils/sam_to_fastq.py +++ b/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) diff --git a/tests/test_long_extraction.py b/tests/test_long_extraction.py new file mode 100644 index 0000000..3d8d306 --- /dev/null +++ b/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])