Feature/binned quality scoring

README.md

@@ -586,6 +586,42 @@ Similarly, use `-i2` to produce a model for paired-ended data. `-q` denotes the
 
 Finally, `-o` is the output directory for the model file and `-p` is the prefix for the output model, such that the file will be written as `<prefix>.p.gz` inside the output folder.
 
+#### Binned quality scoring for modern Illumina instruments
+
+Modern Illumina instruments (NovaSeq 6000, NovaSeq X, NextSeq 2000) compress
+Phred quality scores into a small discrete set of bins rather than emitting a
+continuous range. To train a model that faithfully reproduces this behaviour,
+use either the named `--quality-preset` flag or an explicit bin list via `-Q`.
+
+**Named presets** (recommended):
+
+```bash
+# NovaSeq 6000 / NovaSeq X — Q2, Q12, Q23, Q37
+neat model-qual-score -i reads.fastq.gz --quality-preset novaseq \
+    -o /path/to/models -p novaseq_model
+
+# NextSeq 2000 — Q2, Q12, Q26, Q37
+neat model-qual-score -i reads.fastq.gz --quality-preset nextseq2000 \
+    -o /path/to/models -p nextseq2000_model
+
+# NextSeq 500 / MiniSeq — Q2, Q12, Q23, Q27, Q37
+neat model-qual-score -i reads.fastq.gz --quality-preset nextseq500 \
+    -o /path/to/models -p nextseq500_model
+```
+
+`--quality-preset` implies `--markov`; you do not need to pass both.
+
+**Explicit bin list** — if your instrument uses non-standard bins, pass them
+directly with `-Q`:
+
+```bash
+neat model-qual-score -i reads.fastq.gz -Q 2 12 23 37 --markov \
+    -o /path/to/models -p custom_binned_model
+```
+
+When bins are specified, both the Markov and traditional quality models will
+constrain simulation output to those exact Phred values.
+
 ### `neat model-gc-bias`
 
 Computes GC-bias model from a BAM file and reference genome. It calculates the relative weight of fragments based on their GC content.

neat/cli/commands/model_qual_score.py

@@ -1,6 +1,7 @@
 import argparse
 
 from ...model_quality_score import model_qual_score_runner
+from ...quality_score_modeling.presets import QUALITY_PRESETS
 from .base import BaseCommand
 from .options import output_group
 
@@ -38,7 +39,20 @@ def add_arguments(self, parser: argparse.ArgumentParser):
             type=int,
             nargs="+",
             default=[42],
-            help="Max quality or explicit list of quality scores [42].",
+            help="Maximum quality score (single int) or explicit list of bin values "
+                 "(e.g. -Q 2 12 23 37). A list enables Markov binning: all observed "
+                 "scores are down-binned to the nearest value and simulation output is "
+                 "constrained to those values. Overridden by --quality-preset. [42]",
+        )
+
+        parser.add_argument(
+            "--quality-preset",
+            dest="quality_preset",
+            choices=list(QUALITY_PRESETS),
+            default=None,
+            metavar="PRESET",
+            help="Named bin preset for common Illumina instruments. Implies --markov. "
+                 f"Choices: {', '.join(QUALITY_PRESETS)}. Overrides -Q when set.",
         )
 
         parser.add_argument(
@@ -69,11 +83,15 @@ def add_arguments(self, parser: argparse.ArgumentParser):
 
     def execute(self, arguments: argparse.Namespace):
 
-        if len(arguments.quality_scores) == 1:
-            qual_scores: int | list[int] = arguments.quality_scores[0]
-
+        if arguments.quality_preset:
+            qual_scores: int | list[int] = QUALITY_PRESETS[arguments.quality_preset]
+            use_markov = True
+        elif len(arguments.quality_scores) == 1:
+            qual_scores = arguments.quality_scores[0]
+            use_markov = arguments.use_markov
         else:
             qual_scores = arguments.quality_scores
+            use_markov = arguments.use_markov
 
         model_qual_score_runner(
             files=arguments.input_files,
@@ -83,5 +101,5 @@ def execute(self, arguments: argparse.Namespace):
             overwrite=arguments.overwrite,
             output_dir=arguments.output_dir,
             output_prefix=arguments.prefix,
-            use_markov=arguments.use_markov,
+            use_markov=use_markov,
         )

neat/cli/commands/model_sequencing_error.py

@@ -49,8 +49,9 @@ def add_arguments(self, parser: argparse.ArgumentParser):
                             nargs='+',
                             required=False,
                             default=42,
-                            help="Quality score max. The default 42. The lowest possible score is 1. To used binned"
-                                 "scoring, enter a space separated list of scores, e.g., -Q 2 15 23 37")
+                            help="Maximum quality score [42], or a space-separated list of bin values "
+                                 "(e.g., -Q 2 12 23 37) to constrain the model to those discrete levels. "
+                                 "For named instrument presets use `neat model-qual-score --quality-preset`.")
 
         parser.add_argument('-m',
                             type=int,

neat/model_quality_score/runner.py

@@ -133,6 +133,7 @@ def model_qual_score_runner(
             average_error=average_error,
             quality_scores=np.array(final_quality_scores),
             qual_score_probs=read_parameters[idx],
+            quality_bins=allowed_bins,
         )
 
         markov_model: Optional[MarkovQualityModel] = None

neat/models/error_models.py

@@ -49,14 +49,19 @@ def __init__(
             transition_matrix: np.ndarray = default_error_transition_matrix,
             quality_scores: np.ndarray = default_quality_scores,
             qual_score_probs: np.ndarray = default_qual_score_probs,
-            is_uniform: bool = False
+            is_uniform: bool = False,
+            quality_bins: list[int] | None = None,
     ):
 
         self.transition_matrix = transition_matrix
         self.quality_scores = quality_scores
         self.quality_score_probabilities = qual_score_probs
         self.is_uniform = is_uniform
         self.average_error = average_error
+        # When set, generated scores are snapped to the nearest bin value ≤ the
+        # drawn score (down-binning), matching the discrete quality levels produced
+        # by instruments like NovaSeq (Q2/Q12/Q23/Q37).
+        self.quality_bins: list[int] | None = sorted(quality_bins) if quality_bins else None
 
         # pre-compute the error rate for each quality score. This is the inverse of the phred score equation
         self.quality_score_error_rate: dict[int, float] = {x: 10. ** (-x / 10) for x in self.quality_scores}
@@ -101,7 +106,13 @@ def get_quality_scores(
         means = self.quality_score_probabilities[quality_index_map, 0]
         scales = self.quality_score_probabilities[quality_index_map, 1]
         scores = rng.normal(means, scales)
-        return np.clip(np.rint(scores).astype(int), 1, 42)
+        scores = np.clip(np.rint(scores).astype(int), 1, 42)
+        if self.quality_bins:
+            bins = np.array(self.quality_bins)
+            idx = np.searchsorted(bins, scores, side="right") - 1
+            idx = np.clip(idx, 0, len(bins) - 1)
+            scores = bins[idx]
+        return scores
 
 
 class SequencingErrorModel(SnvModel, DeletionModel, InsertionModel):

neat/quality_score_modeling/markov_utils.py

@@ -17,6 +17,7 @@
 _LOG = logging.getLogger(__name__)
 
 __all__ = [
+    "down_bin_quality",
     "read_quality_lists",
     "compute_initial_distribution",
     "compute_position_distributions",
@@ -25,7 +26,7 @@
 ]
 
 
-def _down_bin_quality(q: int, allowed: List[int]) -> int:
+def down_bin_quality(q: int, allowed: List[int]) -> int:
     """
     Map q to the greatest allowed value <= q (down-binning).
     If q is below the smallest allowed, map to allowed[0].
@@ -110,7 +111,7 @@ def read_quality_lists(
                     continue
 
                 if allowed_sorted is not None:
-                    qlist = [_down_bin_quality(q, allowed_sorted) for q in qlist]
+                    qlist = [down_bin_quality(q, allowed_sorted) for q in qlist]
 
                 qualities.append(qlist)
                 reads_read += 1

neat/quality_score_modeling/presets.py

@@ -0,0 +1,20 @@
+"""
+Named quality-bin presets for common Illumina instruments.
+
+Each preset maps to the discrete Phred score levels that the instrument emits.
+Pass the preset name to ``neat model-qual-score --quality-preset`` instead of
+specifying raw bin values with ``-Q``.
+
+References:
+  NovaSeq 6000 / NovaSeq X: Illumina RTA3 4-level binning (Q2, Q12, Q23, Q37)
+  NextSeq 2000: 4-level binning matching NovaSeq X (Q2, Q12, Q26, Q37)
+  NextSeq 500 / MiniSeq: 5-level binning (Q2, Q12, Q23, Q27, Q37)
+"""
+
+__all__ = ["QUALITY_PRESETS"]
+
+QUALITY_PRESETS: dict[str, list[int]] = {
+    "novaseq":     [2, 12, 23, 37],
+    "nextseq2000": [2, 12, 26, 37],
+    "nextseq500":  [2, 12, 23, 27, 37],
+}

tests/test_models/test_markov_utils.py

@@ -5,7 +5,7 @@
 import pytest
 
 from neat.quality_score_modeling.markov_utils import (
-    _down_bin_quality,
+    down_bin_quality,
     read_quality_lists,
     compute_initial_distribution,
     compute_position_distributions,
@@ -27,33 +27,33 @@ def _write_fastq(path, reads):
 
 
 # ---------------------------------------------------------------------------
-# _down_bin_quality
+# down_bin_quality
 # ---------------------------------------------------------------------------
 
 def test_down_bin_exact_match():
-    assert _down_bin_quality(30, [10, 20, 30, 40]) == 30
+    assert down_bin_quality(30, [10, 20, 30, 40]) == 30
 
 
 def test_down_bin_between_bins_maps_down():
-    assert _down_bin_quality(25, [10, 20, 30, 40]) == 20
+    assert down_bin_quality(25, [10, 20, 30, 40]) == 20
 
 
 def test_down_bin_below_min_maps_to_first_bin():
-    assert _down_bin_quality(5, [10, 20, 30]) == 10
+    assert down_bin_quality(5, [10, 20, 30]) == 10
 
 
 def test_down_bin_above_max_maps_to_last_bin():
-    assert _down_bin_quality(99, [10, 20, 30]) == 30
+    assert down_bin_quality(99, [10, 20, 30]) == 30
 
 
 def test_down_bin_empty_allowed_returns_q_unchanged():
-    assert _down_bin_quality(25, []) == 25
+    assert down_bin_quality(25, []) == 25
 
 
 def test_down_bin_single_bin():
-    assert _down_bin_quality(0, [20]) == 20
-    assert _down_bin_quality(20, [20]) == 20
-    assert _down_bin_quality(40, [20]) == 20
+    assert down_bin_quality(0, [20]) == 20
+    assert down_bin_quality(20, [20]) == 20
+    assert down_bin_quality(40, [20]) == 20
 
 
 # ---------------------------------------------------------------------------

tests/test_models/test_qual_score_models.py

@@ -1,11 +1,14 @@
 """
-Unit tests for MarkovQualityModel
+Unit tests for MarkovQualityModel and TraditionalQualityModel (binned scoring).
 """
 
+import gzip
+import pickle
 import pytest
 import numpy as np
 from numpy.random import default_rng
 
+from neat.models.error_models import TraditionalQualityModel
 from neat.models.markov_quality_model import MarkovQualityModel
 
 # ---------------------------------------------------------------------------
@@ -164,3 +167,84 @@ def test_position_index_for_length_single_position_read():
     """length=1 always returns index 0 regardless of pos."""
     qm = _simple_model(read_length=151)
     assert qm._position_index_for_length(0, 1) == 0
+
+
+# ===========================================================================
+# TraditionalQualityModel — binned scoring
+# ===========================================================================
+
+def _trad_model(quality_bins=None):
+    """Minimal TraditionalQualityModel with a flat μ=30, σ=5 distribution."""
+    read_length = 151
+    qual_score_probs = np.full((read_length, 2), [30.0, 5.0])
+    quality_scores = np.arange(0, 43)
+    return TraditionalQualityModel(
+        quality_scores=quality_scores,
+        qual_score_probs=qual_score_probs,
+        quality_bins=quality_bins,
+    )
+
+
+def test_traditional_model_unbinned_produces_varied_scores():
+    """Without bins, scores span a range (not locked to a small set)."""
+    rng = default_rng(42)
+    model = _trad_model()
+    scores = model.get_quality_scores(151, 151, rng)
+    assert len(scores) == 151
+    assert len(set(scores.tolist())) > 4
+
+
+def test_traditional_model_binned_output_constrained():
+    """All output scores must be members of the supplied bin set."""
+    bins = [2, 12, 23, 37]
+    rng = default_rng(42)
+    model = _trad_model(quality_bins=bins)
+    scores = model.get_quality_scores(151, 151, rng)
+    assert set(scores.tolist()).issubset(set(bins))
+
+
+def test_traditional_model_binned_novaseq_bins():
+    """NovaSeq preset bins produce only {2, 12, 23, 37} across many reads."""
+    bins = [2, 12, 23, 37]
+    rng = default_rng(7)
+    model = _trad_model(quality_bins=bins)
+    all_scores = set()
+    for _ in range(20):
+        all_scores.update(model.get_quality_scores(151, 151, rng).tolist())
+    assert all_scores.issubset(set(bins))
+
+
+def test_traditional_model_binned_below_min_maps_to_first_bin():
+    """A score below the lowest bin should map to the first (lowest) bin."""
+    # Force very low scores: μ=1, σ=0.1 → always clips to 1 → below Q2 → Q2
+    read_length = 10
+    qual_score_probs = np.full((read_length, 2), [1.0, 0.1])
+    model = TraditionalQualityModel(
+        quality_scores=np.arange(0, 43),
+        qual_score_probs=qual_score_probs,
+        quality_bins=[2, 12, 23, 37],
+    )
+    rng = default_rng(0)
+    scores = model.get_quality_scores(read_length, read_length, rng)
+    assert set(scores.tolist()) == {2}
+
+
+def test_traditional_model_binned_persists_through_pickle(tmp_path):
+    """quality_bins must survive a pickle/unpickle round-trip."""
+    bins = [2, 12, 23, 37]
+    model = _trad_model(quality_bins=bins)
+    path = tmp_path / "model.pkl"
+    with open(path, "wb") as f:
+        pickle.dump(model, f)
+    with open(path, "rb") as f:
+        loaded = pickle.load(f)
+    assert loaded.quality_bins == bins
+    rng = default_rng(1)
+    scores = loaded.get_quality_scores(151, 151, rng)
+    assert set(scores.tolist()).issubset(set(bins))
+
+
+def test_traditional_model_none_bins_is_unbinned():
+    """Passing quality_bins=None (explicit) is identical to the default."""
+    model = _trad_model(quality_bins=None)
+    assert model.quality_bins is None