perf(concat): block-copy concatenation instead of SeqRecord round-trip (89x faster, -2.6 GB)

HISTORY.md

@@ -1,5 +1,23 @@
 # History
 
+1.8.5 (2026-09-07)
+------------------
+
+* Concatenates reads and contigs by block copy rather than round-tripping every record through BioPython. `concatenate_single_fastq` parsed both inputs into a list of `SeqRecord`s and wrote them out again, which costs roughly 5-10x the file size in RAM; on 250 MiB of short reads that was 27s and 1.8 GB of peak RSS, and is now 0.7s and 46 MB. Only the non-chromosomal short reads pass through here, so the saving is a transient GB or two rather than the whole run - but it landed immediately before Unicycler, which needs the memory itself. `Plass.get_depth_long` no longer builds its own list of records either. Thanks @[sanjaynagi-eit](https://github.com/sanjaynagi-eit) ([#90](https://github.com/gbouras13/plassembler/pull/90))
+* A wrong-format input is now caught by checking the first byte against the format's record marker (`@` or `>`) instead of by parsing every record, keeping that check at O(1) instead of O(file). One consequence: a truncated *uncompressed* FASTQ is no longer rejected, where full parsing used to raise. Truncated `.gz` input still fails, because the decompressor raises before the copy finishes
+* Concatenation now writes to a `.tmp` sibling and renames it into place. A wrong-format second input, a truncated gzip or a full disk can no longer leave a half-written FASTQ where the run expects a complete one, and a failed re-run leaves any previous output intact
+* An empty `chromosome.fasta` reaching `Assembly.combine_input_fastas` is now a fatal error with a message explaining it. `list(SeqIO.parse(...))[0]` used to raise `IndexError` there; streaming the records has no equivalent, and without an explicit check every depth and copy number would have been computed against a `combined.fasta` containing no chromosome
+
+1.8.4 (2026-08-18)
+------------------
+
+* Fixes the 75bp head and tail cropping of long reads, which had been silently ignored since `chopper` v0.11.0. From v0.11.0, `chopper` only applies `--headcrop`/`--tailcrop` when `--trim-approach fixed-crop` is also specified, so `plassembler` now passes this
+* Bumps the minimum `chopper` version to v0.11.0
+* As a result, filtered long reads are 150bp shorter as originally intended, and reads falling below `--min_length` after cropping are now removed. Expect small changes to long read depths and plasmid copy number estimates compared to v1.8.3
+* Makes `chopper` failures fatal. Previously only the last process in the read filtering pipeline was checked, so a `chopper` that exited non-zero (for example, an old `chopper` rejecting `--trim-approach`) left a valid but empty `chopper_long_reads.fastq.gz` and the assembly continued with zero reads. `plassembler` now checks every stage, reports `chopper`'s own error message rather than only the path to the logfile, rejects empty filtered output, and exits. The per-stage process handling this builds on came from @[sanjaynagi-eit](https://github.com/sanjaynagi-eit) ([#88](https://github.com/gbouras13/plassembler/pull/88)), which also uncovered that `tests/test_data/end_to_end/input_half.fastq.gz` had always ended mid-record, so `chopper` had been failing on it unnoticed
+* `plassembler` now warns if the installed `chopper` is older than v0.11.0, as version pins only bind when an environment is first created
+* Compresses filtered long reads with `bgzip -@` rather than `gzip`, which is roughly 10x faster on a multi-core machine and produces slightly smaller output. BGZF is a valid gzip stream, so nothing downstream changes, and `gzip` is still used where `bgzip` is unavailable. Note that `bgzip` reaches plassembler only via `samtools`' dependency on `htslib`, which does not currently resolve on Apple Silicon - M-series users get the `gzip` fallback and no speedup. Thanks @[sanjaynagi-eit](https://github.com/sanjaynagi-eit) ([#88](https://github.com/gbouras13/plassembler/pull/88))
+
 1.8.3 (2026-07-05)
 ------------------
 

README.md

@@ -232,7 +232,7 @@ You will then need to install the external dependencies separately, which can be
 * [Unicycler](https://github.com/rrwick/Unicycler) >=0.4.8
 * [Minimap2](https://github.com/lh3/minimap2) >=2.11
 * [fastp](https://github.com/OpenGene/fastp) >=0.24.2
-* [chopper](https://github.com/wdecoster/chopper) >=0.5.0
+* [chopper](https://github.com/wdecoster/chopper) >=0.11.0
 * [mash](https://github.com/marbl/Mash) >=2.2
 * [Raven](https://github.com/lbcb-sci/raven) >=1.8
 * [Samtools](https://github.com/samtools/samtools) >=0.15.0

build/environment.yaml

@@ -8,7 +8,7 @@ dependencies:
   - unicycler >=0.4.8
   - minimap2 >=2.11
   - fastp >=0.24.2
-  - chopper >=0.5.0
+  - chopper >=0.11.0
   - mash >=2.2
   - raven-assembler >=1.8
   - samtools >=0.15.0

docs/install.md

@@ -32,7 +32,7 @@ You will then need to install the external dependencies separately, which can be
 * [Unicycler](https://github.com/rrwick/Unicycler) >=0.4.8
 * [Minimap2](https://github.com/lh3/minimap2) >=2.11
 * [fastp](https://github.com/OpenGene/fastp) >=0.24.2
-* [chopper](https://github.com/wdecoster/chopper) >=0.5.0
+* [chopper](https://github.com/wdecoster/chopper) >=0.11.0
 * [mash](https://github.com/marbl/Mash) >=2.2
 * [Raven](https://github.com/lbcb-sci/raven) >=1.8
 * [Samtools](https://github.com/samtools/samtools) >=0.15.0

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "plassembler"
-version = "1.8.3"
+version = "1.8.5"
 description = "Quickly and accurately assemble plasmids in hybrid sequenced bacterial isolates"
 authors = [
     { name = "George Bouras", email = "george.bouras@adelaide.edu.au" }
@@ -75,7 +75,7 @@ flye = ">=2.9"
 unicycler = ">=0.4.8"
 minimap2 = ">=2.11"
 fastp = ">=0.24.2"
-chopper = ">=0.5.0"
+chopper = ">=0.11.0"
 mash = ">=2.2"
 raven-assembler = ">=1.8"
 samtools = ">=0.15.0"

src/plassembler/__init__.py

@@ -491,7 +491,7 @@ def run(
                 longreads,
                 Path(f"{outdir}/chopper_long_reads.fastq"),
             )
-            gzip_file(Path(f"{outdir}/chopper_long_reads.fastq"))
+            gzip_file(Path(f"{outdir}/chopper_long_reads.fastq"), threads)
             remove_file(Path(f"{outdir}/chopper_long_reads.fastq"))
 
     # Raven for long only or '--use_raven'
@@ -1103,7 +1103,7 @@ def assembled(
                         longreads,
                         Path(f"{outdir}/chopper_long_reads.fastq"),
                     )
-                    gzip_file(Path(f"{outdir}/chopper_long_reads.fastq"))
+                    gzip_file(Path(f"{outdir}/chopper_long_reads.fastq"), threads)
                     remove_file(Path(f"{outdir}/chopper_long_reads.fastq"))
 
         if short_flag is True:
@@ -1454,7 +1454,7 @@ def long(
                 longreads,
                 Path(f"{outdir}/chopper_long_reads.fastq"),
             )
-            gzip_file(Path(f"{outdir}/chopper_long_reads.fastq"))
+            gzip_file(Path(f"{outdir}/chopper_long_reads.fastq"), threads)
             remove_file(Path(f"{outdir}/chopper_long_reads.fastq"))
 
     # flye - skip directory an option here

src/plassembler/utils/concat.py

@@ -1,9 +1,12 @@
 import gzip
+import os
 from pathlib import Path
 
-from Bio import SeqIO
 from loguru import logger
 
+# copied between files in 1 MiB blocks
+COPY_BLOCK = 1024 * 1024
+
 
 def concatenate_short_fastqs(out_dir):
     """moves and copies files
@@ -32,49 +35,94 @@ def concatenate_short_fastqs(out_dir):
         logger.error("Error with concatenate_fastqs\n")
 
 
+def _append_file(source: Path, out_handle, record_marker: bytes):
+    """Append source to an open binary handle, decompressing it if gzipped.
+
+    Copies in fixed blocks, so memory does not depend on file size, and ensures
+    the appended block ends with a newline: without that, a source whose final
+    line is unterminated would run into the next file's first header.
+
+    A non-empty file must start with record_marker (b"@" for fastq, b">" for
+    fasta). Parsing every record through BioPython used to catch a wrong-format
+    input; this keeps that check at O(1) instead of O(file).
+
+    :param source: file to append
+    :param out_handle: destination, opened in binary mode
+    :param record_marker: first byte every record of this format starts with
+    :raises ValueError: if source is non-empty and does not start with the marker
+    """
+    opener = gzip.open if Path(source).suffix == ".gz" else open
+    last_byte = b""
+    first = True
+    with opener(source, "rb") as in_handle:
+        while True:
+            block = in_handle.read(COPY_BLOCK)
+            if not block:
+                break
+            if first:
+                if block[:1] != record_marker:
+                    raise ValueError(
+                        f"{source} does not look like a "
+                        f"{record_marker.decode()}-delimited file: it starts with "
+                        f"{block[:1]!r}"
+                    )
+                first = False
+            out_handle.write(block)
+            last_byte = block[-1:]
+    # an empty file contributes nothing, which is normal (e.g. no unmapped reads)
+    if last_byte and last_byte != b"\n":
+        out_handle.write(b"\n")
+
+
+def _concatenate(sources, out_path, record_marker: bytes):
+    """Block-copy sources into out_path, leaving no output behind on failure.
+
+    The copy goes to a sibling .tmp and is renamed into place only once every
+    source has been read, so a wrong-format second input, a truncated gzip or a
+    full disk cannot leave a half-written FASTQ/FASTA where the run expects a
+    complete one. The SeqIO version got this for free by parsing everything
+    before it opened the output.
+
+    :param sources: files to copy, in order
+    :param out_path: destination path
+    :param record_marker: first byte every record of this format starts with
+    """
+    out_path = Path(out_path)
+    tmp_path = out_path.with_name(f"{out_path.name}.tmp")
+    try:
+        with open(tmp_path, "wb") as out_handle:
+            for source in sources:
+                _append_file(source, out_handle, record_marker)
+    except BaseException:
+        tmp_path.unlink(missing_ok=True)
+        raise
+    # same directory, so this is an atomic replace rather than a copy
+    os.replace(tmp_path, out_path)
+
+
 def concatenate_single_fastq(fastq_in1: Path, fastq_in2: Path, fastq_out: Path):
     """concatenates 2 fastq files
+
+    Concatenating reads needs no parsing. Round-tripping them through
+    SeqIO.parse into a list of SeqRecords cost roughly 5-10x the file size in
+    RAM. Only the non-chromosomal short reads pass through here, so that is a
+    transient spike of a GB or two on a typical isolate rather than the whole
+    run - but it landed immediately before Unicycler, which needs the memory
+    itself. A block copy is constant-memory and far faster.
+
     :param fastq_in1:  fastq_in1 input fastq 1
     :param fastq_in2: fastq_in1 input fastq 2
     :param fastq_out: fastq_out output fastq 2
-    :param logger: logger
     :return:
     """
-
-    records = []
-
-    # Read and append records from the first FASTQ file
-    if fastq_in1.suffix == ".gz":
-        with gzip.open(fastq_in1, "rt") as handle:
-            records.extend(SeqIO.parse(handle, "fastq"))
-    else:
-        with open(fastq_in1, "r") as handle:
-            records.extend(SeqIO.parse(handle, "fastq"))
-
-    # Read and append records from the second FASTQ file
-    if fastq_in2.suffix == ".gz":
-        with gzip.open(fastq_in2, "rt") as handle:
-            records.extend(SeqIO.parse(handle, "fastq"))
-    else:
-        with open(fastq_in2, "r") as handle:
-            records.extend(SeqIO.parse(handle, "fastq"))
-
-    # Write the concatenated records to the output FASTQ file
-    with open(fastq_out, "w") as handle:
-        SeqIO.write(records, handle, "fastq")
+    _concatenate([fastq_in1, fastq_in2], fastq_out, b"@")
 
 
 def concatenate_single_fasta(file1: Path, file2: Path, output_file: Path):
-    sequences = []
-
-    # Read sequences from the first file
-    with open(file1, "r") as f1:
-        sequences.extend(SeqIO.parse(f1, "fasta"))
-
-    # Read sequences from the second file
-    with open(file2, "r") as f2:
-        sequences.extend(SeqIO.parse(f2, "fasta"))
-
-    # Write concatenated sequences to the output file
-    with open(output_file, "w") as output:
-        SeqIO.write(sequences, output, "fasta")
+    """concatenates 2 fasta files
+    :param file1: input fasta 1
+    :param file2: input fasta 2
+    :param output_file: output fasta
+    :return:
+    """
+    _concatenate([file1, file2], output_file, b">")

src/plassembler/utils/input_commands.py

@@ -188,6 +188,23 @@ def validate_flye_assembly_info(flye_assembly, flye_info):
     return skip_assembly
 
 
+# plassembler passes `--trim-approach fixed-crop`, which chopper only gained in
+# v0.11.0; older versions reject the flag outright
+MIN_CHOPPER_VERSION = (0, 11, 0)
+
+
+def parse_chopper_version(version_output: str):
+    """Extract (major, minor, patch) from `chopper --version` output.
+
+    :param version_output: stdout of ``chopper --version``, e.g. "chopper 0.11.0".
+    :return: the version as a tuple of ints, or None if it cannot be parsed.
+    """
+    match = re.search(r"(\d+)\.(\d+)\.(\d+)", version_output)
+    if match is None:
+        return None
+    return tuple(int(part) for part in match.groups())
+
+
 def parse_unicycler_version(version_output: str):
     """Extract (major, minor, patch) from `unicycler --version` output.
 
@@ -330,6 +347,8 @@ def check_dependencies():
         logger.error("fastp not found.")
 
     # chopper
+    # bound up front: the gate below still runs if the version could not be read
+    chopper_version = ""
     try:
         process = sp.Popen(["chopper", "--version"], stdout=sp.PIPE, stderr=sp.PIPE)
         chopper_out, _ = process.communicate()
@@ -340,6 +359,21 @@ def check_dependencies():
     except Exception:
         logger.error("chopper not found.")
 
+    # a version pin only binds at install time, so an environment built before the
+    # bump can still hold a chopper that rejects --trim-approach. Warn rather than
+    # exit, so anyone deliberately on an older chopper can still run
+    parsed_chopper_version = parse_chopper_version(chopper_version)
+    min_chopper_version = ".".join(str(part) for part in MIN_CHOPPER_VERSION)
+    if parsed_chopper_version is None:
+        message = f"Could not determine the chopper version from '{chopper_version}'. Plassembler needs chopper >=v{min_chopper_version}."
+        logger.warning(message)
+    elif parsed_chopper_version < MIN_CHOPPER_VERSION:
+        message = f"chopper v{chopper_version} is older than v{min_chopper_version} and will reject the --trim-approach flag Plassembler passes, so long read filtering will fail. Please update chopper, see instructions at https://github.com/gbouras13/plassembler."
+        logger.warning(message)
+    else:
+        message = "chopper version is ok."
+        logger.info(message)
+
     # mash
     try:
         process = sp.Popen(["mash", "version"], stdout=sp.PIPE, stderr=sp.PIPE)