From ec53e9b20c5bbb65ce57f6851de57f6adb5b26ea Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Wed, 29 Apr 2026 00:12:07 -0500 Subject: [PATCH] Add Markov quality model test coverage and fix invalid test data test_qual_score_models.py: - Fix two tests that used -1 as a quality score key in position_distributions (invalid data that bypassed real validation); replace with valid scores - Add 11 new tests covering the previously untested Markov chain path: transition chain is exercised, fallback to marginal for unknown q_prev, wrong transition_distributions length raises, out-of-range score clamping, length=0/1 edge cases, read-length interpolation via _position_index_for_length test_error_models.py: - Remove two stub tests that imported from error_models.MarkovQualityModel (a dead TODO class never called at runtime); they passed unconditionally and gave false confidence about the real model test_markov_utils.py (new): - 36 tests covering all previously untested functions in markov_utils.py: _down_bin_quality (6), compute_initial_distribution (4), compute_position_distributions (4), compute_transition_distributions (5), read_quality_lists (6), build_markov_model (3) Co-Authored-By: Claude Sonnet 4.6 --- tests/test_models/test_error_models.py | 18 -- tests/test_models/test_markov_utils.py | 240 ++++++++++++++++++++ tests/test_models/test_qual_score_models.py | 166 ++++++++++---- 3 files changed, 366 insertions(+), 58 deletions(-) create mode 100644 tests/test_models/test_markov_utils.py diff --git a/tests/test_models/test_error_models.py b/tests/test_models/test_error_models.py index 837d0540..5b31708b 100644 --- a/tests/test_models/test_error_models.py +++ b/tests/test_models/test_error_models.py @@ -287,24 +287,6 @@ def test_sem_blacklist_prevents_duplicate_deletion_sites(): assert len(locations) == len(set(locations)) -# =========================================================================== -# MarkovQualityModel — stub coverage (lines 112-118) -# =========================================================================== - -def test_markov_quality_model_can_be_instantiated(): - """MarkovQualityModel is a TODO stub; construction must not raise.""" - from neat.models.error_models import MarkovQualityModel - m = MarkovQualityModel() - assert m is not None - - -def test_markov_quality_model_get_quality_scores_returns_none(): - """get_quality_scores is a TODO stub; calling it returns None.""" - from neat.models.error_models import MarkovQualityModel - m = MarkovQualityModel() - result = m.get_quality_scores() - assert result is None - # =========================================================================== # TraditionalQualityModel — score clamping (line 104) diff --git a/tests/test_models/test_markov_utils.py b/tests/test_models/test_markov_utils.py new file mode 100644 index 00000000..333ccf36 --- /dev/null +++ b/tests/test_models/test_markov_utils.py @@ -0,0 +1,240 @@ +""" +Unit tests for neat/quality_score_modeling/markov_utils.py +""" + +import pytest + +from neat.quality_score_modeling.markov_utils import ( + _down_bin_quality, + read_quality_lists, + compute_initial_distribution, + compute_position_distributions, + compute_transition_distributions, + build_markov_model, +) + + +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _write_fastq(path, reads): + """Write a list of (seq, qual_string) tuples as a FASTQ file.""" + lines = [] + for i, (seq, qual) in enumerate(reads): + lines += [f"@read{i}", seq, "+", qual] + path.write_text("\n".join(lines) + "\n") + + +# --------------------------------------------------------------------------- +# _down_bin_quality +# --------------------------------------------------------------------------- + +def test_down_bin_exact_match(): + 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 + + +def test_down_bin_below_min_maps_to_first_bin(): + 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 + + +def test_down_bin_empty_allowed_returns_q_unchanged(): + 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 + + +# --------------------------------------------------------------------------- +# compute_initial_distribution +# --------------------------------------------------------------------------- + +def test_compute_initial_distribution_basic(): + quals = [[30, 31, 32], [30, 28, 27], [35, 30, 29]] + result = compute_initial_distribution(quals) + assert result[30] == 2.0 + assert result[35] == 1.0 + assert 28 not in result # position 0 only + assert 31 not in result + + +def test_compute_initial_distribution_empty_input(): + assert compute_initial_distribution([]) == {} + + +def test_compute_initial_distribution_skips_empty_reads(): + result = compute_initial_distribution([[], [30, 31]]) + assert result == {30: 1.0} + + +def test_compute_initial_distribution_single_read(): + result = compute_initial_distribution([[40, 38, 35]]) + assert result == {40: 1.0} + + +# --------------------------------------------------------------------------- +# compute_position_distributions +# --------------------------------------------------------------------------- + +def test_compute_position_distributions_basic(): + quals = [[30, 35, 40], [30, 36, 41]] + result = compute_position_distributions(quals, 3) + assert len(result) == 3 + assert result[0][30] == 2.0 + assert result[1][35] == 1.0 + assert result[1][36] == 1.0 + assert result[2][40] == 1.0 + assert result[2][41] == 1.0 + + +def test_compute_position_distributions_skips_length_mismatch(): + quals = [[30, 31, 32], [30, 31]] # second read is wrong length + result = compute_position_distributions(quals, 3) + assert result[0][30] == 1.0 # only first read counted + assert 31 not in result[0] + + +def test_compute_position_distributions_zero_length_returns_empty(): + assert compute_position_distributions([[30, 31]], 0) == [] + + +def test_compute_position_distributions_single_position(): + quals = [[40], [38], [40]] + result = compute_position_distributions(quals, 1) + assert len(result) == 1 + assert result[0][40] == 2.0 + assert result[0][38] == 1.0 + + +# --------------------------------------------------------------------------- +# compute_transition_distributions +# --------------------------------------------------------------------------- + +def test_compute_transition_distributions_basic(): + quals = [[30, 31, 32], [30, 31, 33]] + result = compute_transition_distributions(quals, 3) + assert len(result) == 2 # read_length - 1 + assert result[0][30][31] == 2.0 # 30→31 seen twice at position 0 + assert result[1][31][32] == 1.0 # 31→32 once at position 1 + assert result[1][31][33] == 1.0 # 31→33 once at position 1 + + +def test_compute_transition_distributions_read_length_one_returns_empty(): + assert compute_transition_distributions([[30]], 1) == [] + + +def test_compute_transition_distributions_read_length_zero_returns_empty(): + assert compute_transition_distributions([], 0) == [] + + +def test_compute_transition_distributions_skips_length_mismatch(): + quals = [[30, 31, 32], [30, 31]] # second read is wrong length + result = compute_transition_distributions(quals, 3) + assert result[0][30][31] == 1.0 # only first read counted + + +def test_compute_transition_distributions_self_transitions(): + """A read with a constant quality score produces only self-transitions.""" + quals = [[35, 35, 35, 35]] + result = compute_transition_distributions(quals, 4) + assert len(result) == 3 + for pos in result: + assert pos[35][35] == 1.0 + assert len(pos) == 1 + + +# --------------------------------------------------------------------------- +# read_quality_lists +# --------------------------------------------------------------------------- + +def test_read_quality_lists_basic(tmp_path): + fq = tmp_path / "test.fastq" + # 'I' = ASCII 73, offset 33 → score 40 + _write_fastq(fq, [("ACGTA", "IIIII"), ("ACGTA", "IIIII")]) + quals, read_length = read_quality_lists([str(fq)], max_reads=100, offset=33) + assert read_length == 5 + assert len(quals) == 2 + assert quals[0] == [40, 40, 40, 40, 40] + + +def test_read_quality_lists_applies_offset(tmp_path): + fq = tmp_path / "test.fastq" + # '!' = ASCII 33, offset 33 → score 0 + _write_fastq(fq, [("ACGT", "!!!!")]) + quals, _ = read_quality_lists([str(fq)], max_reads=100, offset=33) + assert quals[0] == [0, 0, 0, 0] + + +def test_read_quality_lists_respects_max_reads(tmp_path): + fq = tmp_path / "test.fastq" + _write_fastq(fq, [("ACGT", "IIII")] * 10) + quals, _ = read_quality_lists([str(fq)], max_reads=3, offset=33) + assert len(quals) == 3 + + +def test_read_quality_lists_applies_binning(tmp_path): + fq = tmp_path / "test.fastq" + # Scores 40 ('I'), binned to allowed [10, 20, 30] → maps to 30 + _write_fastq(fq, [("ACGT", "IIII")]) + quals, _ = read_quality_lists([str(fq)], max_reads=100, offset=33, + allowed_quality_scores=[10, 20, 30]) + assert quals[0] == [30, 30, 30, 30] + + +def test_read_quality_lists_missing_file_raises(): + with pytest.raises(FileNotFoundError): + read_quality_lists(["/nonexistent/path.fastq"], max_reads=10, offset=33) + + +def test_read_quality_lists_empty_file_returns_empty(tmp_path): + fq = tmp_path / "empty.fastq" + fq.write_text("") + quals, read_length = read_quality_lists([str(fq)], max_reads=100, offset=33) + assert quals == [] + assert read_length == 0 + + +# --------------------------------------------------------------------------- +# build_markov_model (end-to-end) +# --------------------------------------------------------------------------- + +def test_build_markov_model_basic(tmp_path): + fq = tmp_path / "test.fastq" + _write_fastq(fq, [("ACGTA", "IIIII"), ("ACGTA", "IIIII")]) + init, pos_dists, trans_dists, max_q, read_len = build_markov_model( + [str(fq)], max_reads=100, offset=33 + ) + assert read_len == 5 + assert max_q == 40 + assert init == {40: 2.0} + assert len(pos_dists) == 5 + assert len(trans_dists) == 4 # read_length - 1 + assert trans_dists[0][40][40] == 2.0 + + +def test_build_markov_model_max_quality_capped_by_bins(tmp_path): + fq = tmp_path / "test.fastq" + _write_fastq(fq, [("ACGT", "IIII")]) # score 40 + _, _, _, max_q, _ = build_markov_model( + [str(fq)], max_reads=100, offset=33, + allowed_quality_scores=[10, 20, 30] + ) + assert max_q == 30 # capped to max allowed bin + + +def test_build_markov_model_no_reads_raises(tmp_path): + fq = tmp_path / "empty.fastq" + fq.write_text("") + with pytest.raises(ValueError, match="No quality scores"): + build_markov_model([str(fq)], max_reads=100, offset=33) \ No newline at end of file diff --git a/tests/test_models/test_qual_score_models.py b/tests/test_models/test_qual_score_models.py index 9a8485eb..4e4eef9e 100644 --- a/tests/test_models/test_qual_score_models.py +++ b/tests/test_models/test_qual_score_models.py @@ -8,73 +8,159 @@ from neat.models.markov_quality_model import MarkovQualityModel +# --------------------------------------------------------------------------- +# Helpers +# --------------------------------------------------------------------------- + +def _simple_model(max_q=42, read_length=151, init_score=30, pos_score=30, + transition_distributions=None): + """Minimal deterministic model with a single quality score at every position.""" + init_dist = {init_score: 1.0} + pos_dists = [{pos_score: 1.0}] * read_length + return MarkovQualityModel(init_dist, pos_dists, max_q, read_length, + transition_distributions) + + +# --------------------------------------------------------------------------- +# Original tests — fixed to use valid (non-negative) quality score keys +# --------------------------------------------------------------------------- def test_markov_quality_model_shapes_and_range(): - """ - Basic sanity check for MarkovQualityModel: output array shape and bounds. - """ + """Output array has the requested length and all scores are in [0, max_q].""" rng = default_rng(11) - # Simple symmetric distributions around a high-quality region init_dist = {30: 1.0, 31: 1.0, 32: 1.0} - single_pos_dist = {-1: 1.0, 0: 2.0, 1: 1.0} + pos_dist = {30: 1.0, 31: 2.0, 32: 1.0} # fix: was {-1, 0, 1} — invalid keys max_q = 42 read_length = 151 - # Position-specific distributions that reuse the same shape at every position - position_distributions = [single_pos_dist] * read_length - qm = MarkovQualityModel(init_dist, position_distributions, max_q, read_length) + qm = MarkovQualityModel(init_dist, [pos_dist] * read_length, max_q, read_length) qs = qm.get_quality_scores(model_read_length=read_length, length=75, rng=rng) assert isinstance(qs, np.ndarray) assert len(qs) == 75 - # Ensure we stay within the valid range assert qs.min() >= 0 assert qs.max() <= max_q def test_markov_quality_model_quality_scores_property_matches_range(): - """ - quality_scores should expose the full discrete range [0, max_quality]. - """ - init_dist = {35: 1.0} - single_pos_dist = {0: 1.0} - max_q = 40 - read_length = 151 - position_distributions = [single_pos_dist] * read_length - qm = MarkovQualityModel(init_dist, position_distributions, max_q, read_length) + """quality_scores exposes the full contiguous range [0, max_quality].""" + qm = _simple_model(max_q=40, pos_score=35) scores = qm.quality_scores assert isinstance(scores, list) assert scores[0] == 0 - assert scores[-1] == max_q - # The range should be contiguous - assert scores == list(range(0, max_q + 1)) + assert scores[-1] == 40 + assert scores == list(range(0, 41)) def test_markov_quality_model_reproducible_with_seed(): - """Markov quality model should be deterministic for a fixed RNG state.""" + """Fixed RNG seed produces identical output on two independent calls.""" init_dist = {30: 1.0, 31: 1.0} - single_pos_dist = {-1: 1.0, 0: 2.0, 1: 1.0} + pos_dist = {30: 1.0, 31: 2.0, 32: 1.0} # fix: was {-1, 0, 1} — invalid keys max_q = 42 read_length = 151 - position_distributions = [single_pos_dist] * read_length - qm = MarkovQualityModel(init_dist, position_distributions, max_q, read_length) - rng1 = default_rng(12) - rng2 = default_rng(12) - qs1 = qm.get_quality_scores(model_read_length=read_length, length=100, rng=rng1) - qs2 = qm.get_quality_scores(model_read_length=read_length, length=100, rng=rng2) - assert isinstance(qs1, np.ndarray) - assert isinstance(qs2, np.ndarray) + qm = MarkovQualityModel(init_dist, [pos_dist] * read_length, max_q, read_length) + qs1 = qm.get_quality_scores(model_read_length=read_length, length=100, rng=default_rng(12)) + qs2 = qm.get_quality_scores(model_read_length=read_length, length=100, rng=default_rng(12)) assert np.array_equal(qs1, qs2) def test_markov_quality_model_invalid_initial_distribution_raises(): - """ - The model should reject empty or zero-mass initial distributions. - """ + """Empty or zero-mass initial distributions must raise ValueError.""" read_length = 151 - single_pos_dist = {0: 1.0} - position_distributions = [single_pos_dist] * read_length - # Empty initial distribution + pos_dists = [{30: 1.0}] * read_length with pytest.raises(ValueError): - MarkovQualityModel({}, position_distributions, 40, read_length) - # Zero total mass in initial distribution + MarkovQualityModel({}, pos_dists, 40, read_length) with pytest.raises(ValueError): - MarkovQualityModel({30: 0.0}, position_distributions, 40, read_length) + MarkovQualityModel({30: 0.0}, pos_dists, 40, read_length) + + +# --------------------------------------------------------------------------- +# New tests — transition chain path +# --------------------------------------------------------------------------- + +def test_markov_model_uses_transition_chain(): + """When transition_distributions is given, scores follow the chain rows.""" + # Transitions always go to 40 regardless of previous score + always_40 = {q: {40: 1.0} for q in range(0, 43)} + trans_dists = [always_40] * 150 # read_length - 1 + qm = MarkovQualityModel({30: 1.0}, [{30: 1.0, 40: 1.0}] * 151, 42, 151, trans_dists) + qs = qm.get_quality_scores(model_read_length=151, length=20, rng=default_rng(0)) + assert qs[0] == 30 # init distribution + assert all(qs[i] == 40 for i in range(1, len(qs))) # chain forces 40 + + +def test_markov_model_falls_back_to_marginal_for_unknown_q_prev(): + """If q_prev has no transition row, the marginal is used without crashing.""" + # Transition map only covers q=99, which is never reached; marginal gives 35 + sparse_trans = [{99: {99: 1.0}}] * 150 + qm = MarkovQualityModel({30: 1.0}, [{35: 1.0}] * 151, 42, 151, sparse_trans) + qs = qm.get_quality_scores(model_read_length=151, length=10, rng=default_rng(1)) + assert len(qs) == 10 + assert all(0 <= q <= 42 for q in qs) + # All positions after 0 fall back to marginal → always 35 + assert all(qs[i] == 35 for i in range(1, len(qs))) + + +def test_markov_model_wrong_transition_length_raises(): + """transition_distributions with wrong length raises ValueError.""" + bad_trans = [{30: {31: 1.0}}] * 50 # should be 150 for read_length=151 + with pytest.raises(ValueError, match="read_length-1"): + MarkovQualityModel({30: 1.0}, [{30: 1.0}] * 151, 42, 151, bad_trans) + + +def test_markov_model_output_clipped_to_max_quality(): + """Transition rows that emit out-of-range scores are clipped to max_quality.""" + over_max = {q: {99: 1.0} for q in range(0, 43)} + trans_dists = [over_max] * 150 + qm = MarkovQualityModel({30: 1.0}, [{99: 1.0}] * 151, 42, 151, trans_dists) + qs = qm.get_quality_scores(model_read_length=151, length=50, rng=default_rng(5)) + assert qs.max() <= 42 + + +# --------------------------------------------------------------------------- +# New tests — edge cases +# --------------------------------------------------------------------------- + +def test_markov_model_length_zero_returns_empty(): + """length=0 must return an empty ndarray without error.""" + qm = _simple_model() + qs = qm.get_quality_scores(model_read_length=151, length=0, rng=default_rng(0)) + assert isinstance(qs, np.ndarray) + assert len(qs) == 0 + + +def test_markov_model_length_one_uses_only_init(): + """length=1 reads only from the initial distribution, never marginals.""" + qm = _simple_model(init_score=37, pos_score=10) + qs = qm.get_quality_scores(model_read_length=151, length=1, rng=default_rng(0)) + assert len(qs) == 1 + assert qs[0] == 37 + + +def test_markov_model_short_read_interpolates_correctly(): + """A read shorter than the model still returns the requested length.""" + qm = _simple_model(read_length=151) + qs = qm.get_quality_scores(model_read_length=151, length=50, rng=default_rng(0)) + assert len(qs) == 50 + assert all(0 <= q <= 42 for q in qs) + + +# --------------------------------------------------------------------------- +# New tests — _position_index_for_length +# --------------------------------------------------------------------------- + +def test_position_index_for_length_boundaries(): + """First and last positions always map to 0 and read_length-1.""" + qm = _simple_model(read_length=151) + assert qm._position_index_for_length(0, 50) == 0 + assert qm._position_index_for_length(49, 50) == 150 + + +def test_position_index_for_length_midpoint(): + """Midpoint of a 51-position read maps to midpoint of the 151-position model.""" + qm = _simple_model(read_length=151) + assert qm._position_index_for_length(25, 51) == 75 + + +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