#88 plus review fixes

HISTORY.md

@@ -1,5 +1,15 @@
 # History
 
+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.4"
 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/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)

src/plassembler/utils/qc.py

@@ -1,13 +1,89 @@
 import gzip
 import shutil
 import subprocess as sp
+import sys
+from contextlib import ExitStack
 from pathlib import Path
 
 from loguru import logger
 
 from plassembler.utils.external_tools import ExternalTool
 
 
+def gzip_compressor_cmd(threads):
+    """Command that reads plain data on stdin and writes a gzip stream to stdout.
+
+    Prefers bgzip, which ships with htslib/samtools (already a hard plassembler
+    dependency) and compresses in parallel. BGZF is a valid gzip stream, so every
+    downstream reader - flye, minimap2, chopper, gunzip, python's gzip module -
+    handles the result unchanged.
+
+    Serial gzip dominated the chopper stage: on a 300 MB ONT fastq the chain took
+    47.7s, of which only 5.4s was gunzip+chopper and 42s was gzip.
+
+    :param threads: thread count (str or int) to hand to bgzip
+    :return: argv list for the compressor
+    """
+    if shutil.which("bgzip"):
+        return ["bgzip", "-@", str(threads), "-c"]
+    # bgzip should always be present, but never fail QC over a missing binary
+    return ["gzip"]
+
+
+def _chopper_stderr_tail(err_log, max_lines=20):
+    """Last few lines of a chopper logfile, ready to embed in an error message.
+
+    Whatever actually went wrong - an unrecognised flag, a bad value - is written
+    by chopper to its logfile and nowhere else, and users rarely open it. Reading
+    it back puts the real diagnosis where it will be seen.
+
+    :param err_log: path of the chopper .err logfile
+    :param max_lines: how many trailing lines to keep
+    :return: the trailing lines, or "" if the log is missing or empty
+    """
+    try:
+        log_text = Path(err_log).read_text(errors="replace")
+    except OSError:
+        return ""
+    return "\n".join(log_text.strip().splitlines()[-max_lines:])
+
+
+def _fastq_has_reads(filtered_long_reads):
+    """Whether a gzipped fastq actually holds any reads.
+
+    A chopper that dies still leaves a well formed - but empty - gzip member
+    behind, written by the compressor at the end of the chain, so neither the
+    file existing nor it being valid gzip says anything about the run.
+
+    :param filtered_long_reads: path of the gzipped fastq to inspect
+    :return: True if at least one byte of reads can be read back
+    """
+    try:
+        with gzip.open(filtered_long_reads, "rb") as fh:
+            return bool(fh.read(1))
+    except (OSError, EOFError):
+        return False
+
+
+def _fail_chopper(reason, err_log):
+    """Reports a fatal chopper problem and stops the run.
+
+    Everything goes into a single `logger.error`: under the CLI that sink exits
+    the process, so a second call would never be reached.
+
+    :param reason: what went wrong, e.g. "chopper (return code 2)"
+    :param err_log: path of the chopper .err logfile
+    """
+    message = f"Error with chopper: {reason}. Please check {err_log}"
+    stderr_tail = _chopper_stderr_tail(err_log)
+    if stderr_tail:
+        message = f"{message}\nchopper stderr:\n{stderr_tail}"
+    logger.error(message)
+    # the CLI exits inside logger.error above; this covers qc being driven as a
+    # library, where no ERROR sink is installed
+    sys.exit(1)
+
+
 def chopper(
     input_long_reads, outdir, min_length, min_quality, gzip_flag, threads, logdir
 ):
@@ -34,34 +110,81 @@ def chopper(
         threads,
         "-l",
         min_length,
+        # chopper >=0.11.0 only applies --headcrop/--tailcrop under this approach;
+        # without it they are silently ignored
+        "--trim-approach",
+        "fixed-crop",
         "--headcrop",
         "75",
         "--tailcrop",
         "75",
     ]
-    # `with` guarantees the log and output handles are closed even on error;
-    # locals are named *_proc to avoid shadowing the `gzip` module / `chopper`
-    # function name
-    with open(f"{logfile_prefix}.err", "w") as err_log, open(
-        filtered_long_reads, "wb"
-    ) as f:
+    compressor_cmd = gzip_compressor_cmd(threads)
+
+    # ExitStack guarantees the log, output and pipe handles are closed even if a
+    # Popen raises partway through building the chain
+    with ExitStack() as stack:
+        err_log = stack.enter_context(open(f"{logfile_prefix}.err", "w"))
+        out_fh = stack.enter_context(open(filtered_long_reads, "wb"))
+
+        stages = []
         try:
             if gzip_flag is True:
                 source_proc = sp.Popen(
-                    ["gunzip", "-c", input_long_reads], stdout=sp.PIPE
+                    ["gunzip", "-c", input_long_reads], stdout=sp.PIPE, stderr=err_log
                 )
+                stages.append(("gunzip", source_proc))
+                chopper_stdin = source_proc.stdout
             else:
-                source_proc = sp.Popen(["cat", input_long_reads], stdout=sp.PIPE)
+                # plain fastq needs no decompressor: hand the file straight to
+                # chopper rather than spawning a `cat` to copy it through a pipe
+                chopper_stdin = stack.enter_context(open(input_long_reads, "rb"))
+
             chopper_proc = sp.Popen(
-                chopper_cmd,
-                stdin=source_proc.stdout,
-                stdout=sp.PIPE,
-                stderr=err_log,
+                chopper_cmd, stdin=chopper_stdin, stdout=sp.PIPE, stderr=err_log
             )
-            gzip_proc = sp.Popen(["gzip"], stdin=chopper_proc.stdout, stdout=f)
-            gzip_proc.communicate()
-        except Exception:
-            logger.error("Error with chopper")
+            stages.append(("chopper", chopper_proc))
+            # the parent must drop its copy of each upstream read end, otherwise
+            # the downstream stage never sees EOF
+            if gzip_flag is True:
+                source_proc.stdout.close()
+
+            compress_proc = sp.Popen(
+                compressor_cmd, stdin=chopper_proc.stdout, stdout=out_fh, stderr=err_log
+            )
+            stages.append((compressor_cmd[0], compress_proc))
+            chopper_proc.stdout.close()
+        except OSError as e:
+            for _, proc in stages:
+                proc.kill()
+                # reap it too, or the error path leaks the zombies the rest of
+                # this function exists to avoid
+                proc.wait()
+            # _fail_chopper exits, so the stages killed above are never reached
+            # by the wait loop below and reported a second time
+            _fail_chopper(str(e), f"{logfile_prefix}.err")
+
+        # every stage must be waited on. Previously only the last one was, so a
+        # failing chopper was silently ignored and left a zombie behind
+        failures = []
+        for name, proc in reversed(stages):
+            if proc.wait() != 0:
+                failures.append(f"{name} (return code {proc.returncode})")
+
+    if failures:
+        _fail_chopper(", ".join(reversed(failures)), f"{logfile_prefix}.err")
+
+    # a dead stage is not the only way to end up with nothing: whatever the
+    # cause, an empty file here would flow on into Flye/Raven as a zero-read
+    # assembly, so refuse to hand one over
+    if not _fastq_has_reads(filtered_long_reads):
+        _fail_chopper(
+            f"no reads survived filtering. Check that {input_long_reads} holds "
+            f"reads longer than {min_length}bp once 150bp of cropping is applied, "
+            f"with quality above Q{min_quality}",
+            f"{logfile_prefix}.err",
+        )
+
     logger.info("Finished running chopper")
 
 
@@ -104,10 +227,27 @@ def copy_sr_fastq_file(infile: Path, outfile: Path):
         logger.error("Error with copy_sr_fastq_file")
 
 
-def gzip_file(input_path):
+def gzip_file(input_path, threads=1):
+    """gzips a file, in parallel where bgzip is available
+
+    Used by --skip_qc to compress the copied long reads. python's gzip module is
+    both single threaded and slower than the gzip binary, which is a poor fit for
+    a multi-GB ONT fastq; fall back to it only if spawning the compressor fails.
+
+    :param input_path: file to compress
+    :param threads: threads to give the compressor
+    :return: path of the compressed file
+    """
     input_path = Path(input_path)
     output_path = input_path.with_suffix(input_path.suffix + ".gz")
 
+    try:
+        with open(input_path, "rb") as f_in, open(output_path, "wb") as f_out:
+            sp.run(gzip_compressor_cmd(threads), stdin=f_in, stdout=f_out, check=True)
+        return output_path
+    except (OSError, sp.CalledProcessError) as e:
+        logger.warning(f"Falling back to python gzip for {input_path}: {e}")
+
     with open(input_path, "rb") as f_in:
         with gzip.open(output_path, "wb") as f_out:
             shutil.copyfileobj(f_in, f_out)

tests/test_plassembler.py

@@ -30,7 +30,9 @@
 
 # import functions
 from src.plassembler.utils.input_commands import (
+    MIN_CHOPPER_VERSION,
     check_dependencies,
+    parse_chopper_version,
     parse_unicycler_version,
     validate_fasta,
     validate_fastas_assembled_mode,
@@ -199,6 +201,23 @@ def test_parse_unicycler_version_missing_raises(self):
         with self.assertRaises(ValueError):
             parse_unicycler_version("bash: unicycler: command not found\n")
 
+    # chopper version parsing
+    def test_parse_chopper_version_clean(self):
+        self.assertEqual(parse_chopper_version("chopper 0.11.0\n"), (0, 11, 0))
+
+    def test_parse_chopper_version_unparseable(self):
+        # check_dependencies warns on None rather than crashing, so a missing or
+        # reworded version string must come back as None, not raise
+        self.assertIsNone(parse_chopper_version("chopper not found"))
+        self.assertIsNone(parse_chopper_version(""))
+
+    def test_chopper_minimum_version_ordering(self):
+        # plassembler passes --trim-approach, added in chopper v0.11.0
+        self.assertLess(parse_chopper_version("chopper 0.10.0"), MIN_CHOPPER_VERSION)
+        self.assertGreaterEqual(
+            parse_chopper_version("chopper 0.13.0"), MIN_CHOPPER_VERSION
+        )
+
     # checks all external dependencies are installed
     @pytest.mark.slow
     def test_deps(self):