diff --git a/README.md b/README.md index 28a8c680..2b676ee1 100755 --- a/README.md +++ b/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 `.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. diff --git a/neat/cli/commands/model_qual_score.py b/neat/cli/commands/model_qual_score.py index aea15edf..422aa15d 100644 --- a/neat/cli/commands/model_qual_score.py +++ b/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, ) diff --git a/neat/cli/commands/model_sequencing_error.py b/neat/cli/commands/model_sequencing_error.py index 4c3420b0..c9527219 100644 --- a/neat/cli/commands/model_sequencing_error.py +++ b/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, diff --git a/neat/model_quality_score/runner.py b/neat/model_quality_score/runner.py index 2c355c7b..0e9ca77a 100644 --- a/neat/model_quality_score/runner.py +++ b/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 diff --git a/neat/models/error_models.py b/neat/models/error_models.py index c9ed0793..c3e44411 100644 --- a/neat/models/error_models.py +++ b/neat/models/error_models.py @@ -49,7 +49,8 @@ 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 @@ -57,6 +58,10 @@ def __init__( 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): diff --git a/neat/quality_score_modeling/markov_utils.py b/neat/quality_score_modeling/markov_utils.py index d75c7a47..75f9088e 100644 --- a/neat/quality_score_modeling/markov_utils.py +++ b/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 diff --git a/neat/quality_score_modeling/presets.py b/neat/quality_score_modeling/presets.py new file mode 100644 index 00000000..b6b1f40f --- /dev/null +++ b/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], +} \ No newline at end of file diff --git a/tests/test_models/test_markov_utils.py b/tests/test_models/test_markov_utils.py index 333ccf36..c45389d9 100644 --- a/tests/test_models/test_markov_utils.py +++ b/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 # --------------------------------------------------------------------------- diff --git a/tests/test_models/test_qual_score_models.py b/tests/test_models/test_qual_score_models.py index 4e4eef9e..192e91ea 100644 --- a/tests/test_models/test_qual_score_models.py +++ b/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 diff --git a/tests/test_models/test_quality_presets.py b/tests/test_models/test_quality_presets.py new file mode 100644 index 00000000..d234a878 --- /dev/null +++ b/tests/test_models/test_quality_presets.py @@ -0,0 +1,71 @@ +""" +Tests for neat/quality_score_modeling/presets.py and the --quality-preset +CLI wiring in neat/cli/commands/model_qual_score.py. +""" +import argparse + +import pytest + +from neat.quality_score_modeling.presets import QUALITY_PRESETS + + +# =========================================================================== +# Preset data integrity +# =========================================================================== + +def test_all_presets_are_sorted(): + for name, bins in QUALITY_PRESETS.items(): + assert bins == sorted(bins), f"preset {name!r} bins are not sorted" + + +def test_all_preset_bins_are_positive(): + for name, bins in QUALITY_PRESETS.items(): + assert all(b > 0 for b in bins), f"preset {name!r} has non-positive bin" + + +def test_novaseq_has_four_bins(): + assert len(QUALITY_PRESETS["novaseq"]) == 4 + + +def test_nextseq2000_has_four_bins(): + assert len(QUALITY_PRESETS["nextseq2000"]) == 4 + + +def test_nextseq500_has_five_bins(): + assert len(QUALITY_PRESETS["nextseq500"]) == 5 + + +def test_novaseq_bins_match_illumina_spec(): + assert QUALITY_PRESETS["novaseq"] == [2, 12, 23, 37] + + +def test_nextseq2000_bins_match_illumina_spec(): + assert QUALITY_PRESETS["nextseq2000"] == [2, 12, 26, 37] + + +# =========================================================================== +# CLI --quality-preset wiring +# =========================================================================== + +def _parse_qual_score_args(argv): + """Run the model-qual-score argument parser against argv.""" + from neat.cli.commands.model_qual_score import Command + parser = argparse.ArgumentParser() + Command(parser) + return parser.parse_args(argv) + + +def test_quality_preset_accepted_by_cli(): + args = _parse_qual_score_args(["-i", "f.fq", "-o", "/tmp", "--quality-preset", "novaseq"]) + assert args.quality_preset == "novaseq" + + +def test_unknown_preset_rejected_by_cli(): + with pytest.raises(SystemExit): + _parse_qual_score_args(["-i", "f.fq", "-o", "/tmp", "--quality-preset", "bad_preset"]) + + +def test_all_preset_names_accepted_by_cli(): + for name in QUALITY_PRESETS: + args = _parse_qual_score_args(["-i", "f.fq", "-o", "/tmp", "--quality-preset", name]) + assert args.quality_preset == name \ No newline at end of file