From 44a470ebd12ad6166ef46a56f304656f70ae5fae Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Tue, 11 Nov 2025 23:09:01 -0600 Subject: [PATCH 1/7] Updating model, read names, and a few other details. Overall goal is to implement avg_seq_error in such a way that actually affects the output. --- neat/models/error_models.py | 23 ++++++++++++++----- neat/read_simulator/runner.py | 25 +++++++++++++++++++++ neat/read_simulator/single_runner.py | 13 +++++++++-- neat/read_simulator/utils/generate_reads.py | 14 ++++++++---- neat/read_simulator/utils/read.py | 7 +++++- 5 files changed, 70 insertions(+), 12 deletions(-) diff --git a/neat/models/error_models.py b/neat/models/error_models.py index 386cbf5c..457c8b60 100644 --- a/neat/models/error_models.py +++ b/neat/models/error_models.py @@ -10,6 +10,7 @@ from Bio.Seq import Seq from Bio import SeqRecord +from numpy import median from neat import variants @@ -165,8 +166,9 @@ def __init__( def get_sequencing_errors( self, padding: int, - reference_segment: SeqRecord, + reference_segment: Seq, quality_scores: np.ndarray, + num_errors, rng ): """ @@ -175,8 +177,9 @@ def get_sequencing_errors( :param padding: this is the amount of space we have in the read for deletions. :param reference_segment: The section of the reference from which the read is drawn :param quality_scores: Array of quality scores for the read - :return: Modified sequence and associated quality scores + :param num_errors: The estimated number of errors to add. :param rng: random number generator. + :return: Modified sequence and associated quality scores """ error_indexes = [] @@ -189,9 +192,19 @@ def get_sequencing_errors( if self.average_error == 0: return introduced_errors else: - for i in range(len(quality_scores)): - if rng.random() < quality_score_error_rate[quality_scores[i]]: - error_indexes.append(i) + i = len(quality_scores) + while len(error_indexes) <= num_errors and i > 0: + index = rng.choice(list(range(len(quality_scores)))) + if rng.random() < quality_score_error_rate[quality_scores[index]]: + error_indexes.append(index) + i -= 1 + # This should fill in any errors to make sure we aren't coming up short + median_score = median(quality_scores) + while len(error_indexes) < num_errors: + index = rng.choice(quality_scores) + score = quality_scores[index] + if score < median_score: + error_indexes.append(index) total_indel_length = 0 # To prevent deletion collisions diff --git a/neat/read_simulator/runner.py b/neat/read_simulator/runner.py index 424b6ed2..c47cba6c 100644 --- a/neat/read_simulator/runner.py +++ b/neat/read_simulator/runner.py @@ -1,11 +1,15 @@ """ Runner for generate_reads task """ +import gzip import logging +import os +import pickle import shutil import subprocess import time import multiprocessing as mp +from math import ceil from pathlib import Path @@ -16,6 +20,7 @@ from .utils import Options, OutputFileWriter, parse_beds, parse_input_vcf from ..common import validate_input_path, validate_output_path from .single_runner import read_simulator_single +from ..models import SequencingErrorModel from ..variants import ContigVariants from .utils.split_inputs import main as split_main from .utils.stitch_outputs import main as stitch_main @@ -74,6 +79,25 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str): reference_index = SeqIO.index(str(options.reference), "fasta") reference_keys_with_lens = {key: len(value) for key, value in reference_index.items()} + # We need sequencing errors to get the quality score attributes, even for the vcf + if options.error_model: + error_models = pickle.load(gzip.open(options.error_model)) + error_model = error_models["error_model1"] + else: + # Use all the default values + error_model = SequencingErrorModel() + + # Update error to user specified input + if options.avg_seq_error: + error_model.average_error = options.avg_seq_error + + # _LOG.debug('Sequencing error and quality score models loaded') + # We need to estimate how many total errors to add + total_reference_length = sum(reference_keys_with_lens.values()) + total_errors = ceil(error_model.average_error * total_reference_length) + normalized_counts = {k: v / total_reference_length for (k, v) in reference_keys_with_lens.items()} + errors_per_contig = {k: ceil(v * total_errors) for (k, v) in normalized_counts.items()} + count = 0 for contig in reference_keys_with_lens: count += reference_keys_with_lens[contig] @@ -175,6 +199,7 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str): target_regions_dict[contig], discard_regions_dict[contig], mutation_rate_dict[contig], + errors_per_contig[contig], ) _LOG.info(f"Completed simulating contig {contig}.") # TODO Remove if not needed diff --git a/neat/read_simulator/single_runner.py b/neat/read_simulator/single_runner.py index 57965b8d..965a239f 100644 --- a/neat/read_simulator/single_runner.py +++ b/neat/read_simulator/single_runner.py @@ -11,6 +11,8 @@ import logging from pathlib import Path +from math import ceil, floor + from .utils import OutputFileWriter, \ generate_variants, generate_reads, Options, recalibrate_mutation_regions from ..variants import ContigVariants @@ -25,20 +27,21 @@ def read_simulator_single( thread_idx: int, block_start: int, local_options: Options, - bam_header: list | None, + bam_header: dict | None, contig_name: str, contig_index: int, input_variants_local: ContigVariants, target_regions: list, discard_regions: list, mutation_regions: list, + errors_for_contig: int, ) -> tuple[int, str, ContigVariants, dict[str, Path], ]: """ inputs: :param thread_idx: index of current thread :param block_start: Where on the reference does this block start? For a full contig, this will be 0. :param local_options: options for current thread and reference chunk - :param bam_header: Pass in the outer bam header. + :param bam_header: Pass in the outer bam header, which is a dictionary of the contigs plus lengths. :param contig_name: The original list of contig names. :param contig_index: The index of the contig which this chunk comes from :param input_variants_local: The input variants for this block @@ -47,6 +50,7 @@ def read_simulator_single( :param target_regions: Target regions for the run :param discard_regions: discard regions for the run :param mutation_regions: mutation regions (unchecked) for the run + :param errors_for_contig: How many errors this contig will receive Ideally this should work for either a file chunk or contig. We'll assume here that we're getting either an entire contig or a file chunk, and that no new subdivisions are needed. @@ -108,11 +112,16 @@ def read_simulator_single( options=local_options, ) + # This gives the percentage of the contig this particular block is + contig_percentage = len(local_seq_record) * bam_header[contig_name] + # How many errors to add to this block + errors_to_add = ceil(contig_percentage * errors_for_contig) if local_options.produce_fastq or local_options.produce_bam: reads_to_write = generate_reads( thread_idx, local_seq_record, seq_error_model, + errors_to_add, qual_score_model, fraglen_model, local_variants, diff --git a/neat/read_simulator/utils/generate_reads.py b/neat/read_simulator/utils/generate_reads.py index 52bb619e..49c32d47 100644 --- a/neat/read_simulator/utils/generate_reads.py +++ b/neat/read_simulator/utils/generate_reads.py @@ -149,6 +149,7 @@ def generate_reads( thread_index: int, reference: SeqRecord, error_model: SequencingErrorModel, + errors_in_contig: int, qual_model: TraditionalQualityModel, fraglen_model: FragmentLengthModel, contig_variants: ContigVariants, @@ -166,6 +167,7 @@ def generate_reads( :param thread_index: Index of current thread :param reference: The reference segment that reads will be drawn from. :param error_model: The error model for this run, the forward strand + :param errors_in_contig: Total number of errors to add to contig :param qual_model: The quality score model for this run, forward strand :param fraglen_model: The fragment length model for this run :param contig_variants: An object containing all input and randomly generated variants to be included. @@ -258,6 +260,8 @@ def generate_reads( padding = options.read_len//5 segment = reference[read1[0]: read1[1] + padding].seq + errors_per_read = max((options.read_len // len(reference)) * errors_in_contig, 1) + # if we're at the end of the contig, this may not pick up the full padding actual_padding = len(segment) - options.read_len @@ -279,15 +283,16 @@ def generate_reads( fastq_handle = ofw.files_to_write[ofw.fq1] else: fastq_handle = None - read_1.finalize_read_and_write( + num_errors = read_1.finalize_read_and_write( error_model, qual_model, fastq_handle, options.quality_offset, options.produce_fastq, + errors_per_read, options.rng ) - + errors_in_contig -= num_errors # skip over read 2 for single ended reads. if options.paired_ended: # Padding, as above @@ -313,19 +318,20 @@ def generate_reads( ) read_2.mutations = find_applicable_mutations(read_2, contig_variants) - if options.produce_fastq: fastq_handle = ofw.files_to_write[ofw.fq2] else: fastq_handle = None - read_2.finalize_read_and_write( + num_errors = read_2.finalize_read_and_write( error_model, qual_model, fastq_handle, options.quality_offset, options.produce_fastq, + errors_per_read, options.rng ) + errors_in_contig -= num_errors reads_to_write.append((read_1, read_2)) else: reads_to_write.append((read_1, None)) diff --git a/neat/read_simulator/utils/read.py b/neat/read_simulator/utils/read.py index 6b81d79c..cd14ea57 100644 --- a/neat/read_simulator/utils/read.py +++ b/neat/read_simulator/utils/read.py @@ -308,6 +308,7 @@ def finalize_read_and_write( fastq_handle, quality_offset: int, produce_fastq: bool, + num_errors: int, rng: Generator, ): """ @@ -319,6 +320,7 @@ def finalize_read_and_write( :param quality_offset: the quality offset for this run :param produce_fastq: If true, this will write out the temp fastqs. If false, this will only write out the tsams to create the bam files. + :param num_errors: Estimated number of errors to add to this read. :param rng: the random number generator for this run """ @@ -336,11 +338,12 @@ def finalize_read_and_write( # set the read sequence to match the reference, then modify self.read_sequence = deepcopy(self.reference_segment) - # Get errors for the rea and update the quality score + # Get errors for the read and update the quality score self.errors, self.padding = err_model.get_sequencing_errors( self.padding, self.reference_segment, self.quality_array, + num_errors, rng ) @@ -362,6 +365,8 @@ def finalize_read_and_write( fastq_record = f'@{self.name}\n{str(self.read_sequence)}\n+\n{self.read_quality_string}\n' fastq_handle.write(fastq_record) + return len(self.errors) + def convert_masking(self, quality_model: TraditionalQualityModel): """ From 308cbd9559cc2d9df3e304f351167a6a57faaf2d Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Wed, 12 Nov 2025 08:26:42 -0600 Subject: [PATCH 2/7] Expanding error model somewhat --- neat/read_simulator/runner.py | 37 +++++++++++++++------ neat/read_simulator/single_runner.py | 12 ++----- neat/read_simulator/utils/generate_reads.py | 8 ++--- neat/read_simulator/utils/options.py | 6 +++- 4 files changed, 36 insertions(+), 27 deletions(-) diff --git a/neat/read_simulator/runner.py b/neat/read_simulator/runner.py index c47cba6c..a2834a0e 100644 --- a/neat/read_simulator/runner.py +++ b/neat/read_simulator/runner.py @@ -80,22 +80,28 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str): reference_keys_with_lens = {key: len(value) for key, value in reference_index.items()} # We need sequencing errors to get the quality score attributes, even for the vcf - if options.error_model: + if options.avg_seq_error: + average_error = options.avg_seq_error + elif options.error_model: error_models = pickle.load(gzip.open(options.error_model)) - error_model = error_models["error_model1"] + average_error = error_models["error_model1"].average_error + # We just need the error value + del error_models else: - # Use all the default values - error_model = SequencingErrorModel() - - # Update error to user specified input - if options.avg_seq_error: - error_model.average_error = options.avg_seq_error + # Use the default value + average_error = 0.009228843915252066 # _LOG.debug('Sequencing error and quality score models loaded') # We need to estimate how many total errors to add total_reference_length = sum(reference_keys_with_lens.values()) - total_errors = ceil(error_model.average_error * total_reference_length) + # This is an estimate due to some random effects + number_of_bases_in_analysis = total_reference_length * options.coverage + # Each base called has an "average_error" chance of being an error + total_errors = ceil(average_error * number_of_bases_in_analysis) + # Normalization gives the percent value for each item with a sum of all values being 1.0 normalized_counts = {k: v / total_reference_length for (k, v) in reference_keys_with_lens.items()} + # Multiply the normalized count by the total errors. This gives the number of errors that should be + # introduced into the contig errors_per_contig = {k: ceil(v * total_errors) for (k, v) in normalized_counts.items()} count = 0 @@ -165,6 +171,14 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str): for contig in splits_files_dict: contig_index = contig_dict[contig] for ((start, length), splits_file) in splits_files_dict[contig].items(): + block_percentage = length / reference_keys_with_lens[contig] + block_errors = errors_per_contig[contig] * block_percentage + estimated_number_of_reads = (length // options.read_len) * options.coverage + errors_per_read = round(block_errors / estimated_number_of_reads) + if errors_per_read < 1.0 and block_errors > 0: + # We know we need a few errors, but it's a small number total + if options.rng.random() < average_error: + errors_per_read += 1 current_output_dir = options.temp_dir_path / splits_file.stem current_output_dir.mkdir(parents=True, exist_ok=True) # Create local filenames based on fasta indexing scheme. @@ -199,7 +213,7 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str): target_regions_dict[contig], discard_regions_dict[contig], mutation_rate_dict[contig], - errors_per_contig[contig], + errors_per_read, ) _LOG.info(f"Completed simulating contig {contig}.") # TODO Remove if not needed @@ -221,7 +235,8 @@ def read_simulator_runner(config: str, output_dir: str, file_prefix: str): thread_input_variants, thread_target_regions, thread_discard_regions, - thread_mutation_regions + thread_mutation_regions, + errors_per_read, )) thread_idx += 1 diff --git a/neat/read_simulator/single_runner.py b/neat/read_simulator/single_runner.py index 965a239f..1db8d762 100644 --- a/neat/read_simulator/single_runner.py +++ b/neat/read_simulator/single_runner.py @@ -34,7 +34,7 @@ def read_simulator_single( target_regions: list, discard_regions: list, mutation_regions: list, - errors_for_contig: int, + errors_per_read: int, ) -> tuple[int, str, ContigVariants, dict[str, Path], ]: """ inputs: @@ -45,12 +45,10 @@ def read_simulator_single( :param contig_name: The original list of contig names. :param contig_index: The index of the contig which this chunk comes from :param input_variants_local: The input variants for this block - TODO I'm counting on the target and discard regions not being used. They are likely broken with multithreading. - Probably they make more sense after the fact now :param target_regions: Target regions for the run :param discard_regions: discard regions for the run :param mutation_regions: mutation regions (unchecked) for the run - :param errors_for_contig: How many errors this contig will receive + :param errors_per_read: How many errors this contig will receive Ideally this should work for either a file chunk or contig. We'll assume here that we're getting either an entire contig or a file chunk, and that no new subdivisions are needed. @@ -112,16 +110,12 @@ def read_simulator_single( options=local_options, ) - # This gives the percentage of the contig this particular block is - contig_percentage = len(local_seq_record) * bam_header[contig_name] - # How many errors to add to this block - errors_to_add = ceil(contig_percentage * errors_for_contig) if local_options.produce_fastq or local_options.produce_bam: reads_to_write = generate_reads( thread_idx, local_seq_record, seq_error_model, - errors_to_add, + errors_per_read, qual_score_model, fraglen_model, local_variants, diff --git a/neat/read_simulator/utils/generate_reads.py b/neat/read_simulator/utils/generate_reads.py index 49c32d47..6b5cbfba 100644 --- a/neat/read_simulator/utils/generate_reads.py +++ b/neat/read_simulator/utils/generate_reads.py @@ -149,7 +149,7 @@ def generate_reads( thread_index: int, reference: SeqRecord, error_model: SequencingErrorModel, - errors_in_contig: int, + errors_per_read: int, qual_model: TraditionalQualityModel, fraglen_model: FragmentLengthModel, contig_variants: ContigVariants, @@ -167,7 +167,7 @@ def generate_reads( :param thread_index: Index of current thread :param reference: The reference segment that reads will be drawn from. :param error_model: The error model for this run, the forward strand - :param errors_in_contig: Total number of errors to add to contig + :param errors_per_read: Total number of errors to add to contig :param qual_model: The quality score model for this run, forward strand :param fraglen_model: The fragment length model for this run :param contig_variants: An object containing all input and randomly generated variants to be included. @@ -260,8 +260,6 @@ def generate_reads( padding = options.read_len//5 segment = reference[read1[0]: read1[1] + padding].seq - errors_per_read = max((options.read_len // len(reference)) * errors_in_contig, 1) - # if we're at the end of the contig, this may not pick up the full padding actual_padding = len(segment) - options.read_len @@ -292,7 +290,6 @@ def generate_reads( errors_per_read, options.rng ) - errors_in_contig -= num_errors # skip over read 2 for single ended reads. if options.paired_ended: # Padding, as above @@ -331,7 +328,6 @@ def generate_reads( errors_per_read, options.rng ) - errors_in_contig -= num_errors reads_to_write.append((read_1, read_2)) else: reads_to_write.append((read_1, None)) diff --git a/neat/read_simulator/utils/options.py b/neat/read_simulator/utils/options.py index 17abf39e..6aa2ee2b 100644 --- a/neat/read_simulator/utils/options.py +++ b/neat/read_simulator/utils/options.py @@ -286,7 +286,11 @@ def check_and_log_error(keyname, value_to_check, crit1, crit2): if value_to_check not in crit2: _LOG.error(f"Must choose one of {crit2}") sys.exit(1) - elif isinstance(crit1, (int, float)) and isinstance(crit2, (int, float)): + elif not crit1 and not crit2: + # Nothing to check + pass + else: + # Must be a range if not (crit1 <= value_to_check <= crit2): _LOG.error(f'`{keyname}` must be between {crit1} and {crit2} (input: {value_to_check}).') sys.exit(1) From a78be157dd3073235adbdc1f02f59d334125c458 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Thu, 13 Nov 2025 08:12:45 -0600 Subject: [PATCH 3/7] Updating some parameters in error models --- neat/models/error_models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/neat/models/error_models.py b/neat/models/error_models.py index 457c8b60..60da5c28 100644 --- a/neat/models/error_models.py +++ b/neat/models/error_models.py @@ -223,11 +223,11 @@ def get_sequencing_errors( # Deletion error if error_type == Deletion: - deletion_length = self.get_deletion_length() + deletion_length = self.get_deletion_length(rng) if padding - deletion_length < 0: # No space in this read to add this deletion continue - deletion_reference = reference_segment.seq[index: index + deletion_length + 1] + deletion_reference = reference_segment[index: index + deletion_length + 1] deletion_alternate = deletion_reference[0] introduced_errors.append( ErrorContainer(Deletion, index, deletion_length, deletion_reference, deletion_alternate) From 699815261d7895b8d0a61f742636f60507b6b3ae Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Sun, 3 May 2026 08:41:40 -0500 Subject: [PATCH 4/7] Fix off-by-one and index bug in get_sequencing_errors The main loop used <= instead of < causing num_errors+1 errors to be collected. The fallback loop called rng.choice(quality_scores) which returns a score value, not an index, making quality_scores[index] wrong; replaced with rng.integers(len(quality_scores)). Co-Authored-By: Claude Sonnet 4.6 --- neat/models/error_models.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/neat/models/error_models.py b/neat/models/error_models.py index 60da5c28..2eb23082 100644 --- a/neat/models/error_models.py +++ b/neat/models/error_models.py @@ -193,17 +193,17 @@ def get_sequencing_errors( return introduced_errors else: i = len(quality_scores) - while len(error_indexes) <= num_errors and i > 0: + while len(error_indexes) < num_errors and i > 0: index = rng.choice(list(range(len(quality_scores)))) if rng.random() < quality_score_error_rate[quality_scores[index]]: error_indexes.append(index) i -= 1 - # This should fill in any errors to make sure we aren't coming up short + # Fallback: if quality scores are too high to naturally reach num_errors, + # force errors at positions with below-median quality scores median_score = median(quality_scores) while len(error_indexes) < num_errors: - index = rng.choice(quality_scores) - score = quality_scores[index] - if score < median_score: + index = rng.integers(len(quality_scores)) + if quality_scores[index] < median_score: error_indexes.append(index) total_indel_length = 0 From 973714f65401d836c8eae5b8ad55505f1e5fffee Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Sun, 3 May 2026 09:40:09 -0500 Subject: [PATCH 5/7] Update tests for new get_sequencing_errors and errors_per_read signatures - Add num_errors argument to all get_sequencing_errors() calls - Change reference_segment from SeqRecord to Seq in error model tests - Add errors_per_read argument to generate_reads() calls - Add num_errors argument to finalize_read_and_write() calls - Add errors_per_read argument to read_simulator_single() calls - Remove now-unused SeqRecord imports from test_error_and_mut_models and test_seq_error Co-Authored-By: Claude Sonnet 4.6 --- .../test_models/test_error_and_mut_models.py | 15 ++++++------- tests/test_models/test_error_models.py | 22 +++++++++---------- tests/test_models/test_seq_error.py | 5 ++--- .../test_generate_reads.py | 18 +++++++-------- tests/test_read_simulator/test_read.py | 10 ++++----- .../test_read_simulator/test_single_runner.py | 3 +++ 6 files changed, 37 insertions(+), 36 deletions(-) diff --git a/tests/test_models/test_error_and_mut_models.py b/tests/test_models/test_error_and_mut_models.py index 6fd2f83f..3c888620 100644 --- a/tests/test_models/test_error_and_mut_models.py +++ b/tests/test_models/test_error_and_mut_models.py @@ -5,7 +5,6 @@ import numpy as np from numpy.random import default_rng from Bio.Seq import Seq -from Bio.SeqRecord import SeqRecord from neat.models.mutation_model import MutationModel from neat.models.error_models import SequencingErrorModel, TraditionalQualityModel @@ -57,10 +56,10 @@ def test_traditional_quality_model_shapes_and_range(): def test_sequencing_error_model_basic_snvs_only(): rng = default_rng(6) sem = SequencingErrorModel(avg_seq_error=0.05) - ref = SeqRecord(Seq("ACGT" * 20), id="chr1") + ref = Seq("ACGT" * 20) quals = np.array([35] * 80, dtype=int) introduced, pad = sem.get_sequencing_errors( - padding=40, reference_segment=ref, quality_scores=quals, rng=rng + padding=40, reference_segment=ref, quality_scores=quals, num_errors=3, rng=rng ) assert isinstance(introduced, list) assert pad >= 0 @@ -124,17 +123,17 @@ def test_mutation_model_snv_does_not_keep_reference_base(): def test_sequencing_error_model_reproducible_with_seed(): """Error placement should be deterministic given the same RNG state.""" sem = SequencingErrorModel(avg_seq_error=0.05) - ref = SeqRecord(Seq("ACGT" * 30), id="chr1") + ref = Seq("ACGT" * 30) quals = np.array([30] * len(ref), dtype=int) rng1 = default_rng(9) rng2 = default_rng(9) introduced1, pad1 = sem.get_sequencing_errors( - padding=20, reference_segment=ref, quality_scores=quals, rng=rng1 + padding=20, reference_segment=ref, quality_scores=quals, num_errors=3, rng=rng1 ) introduced2, pad2 = sem.get_sequencing_errors( - padding=20, reference_segment=ref, quality_scores=quals, rng=rng2 + padding=20, reference_segment=ref, quality_scores=quals, num_errors=3, rng=rng2 ) assert pad1 == pad2 @@ -150,11 +149,11 @@ def test_sequencing_error_model_nonzero_error_introduces_in_bounds_errors(): """ rng = default_rng(10) sem = SequencingErrorModel(avg_seq_error=0.2) - ref = SeqRecord(Seq("ACGT" * 50), id="chr1") + ref = Seq("ACGT" * 50) quals = np.array([10] * len(ref), dtype=int) introduced, pad = sem.get_sequencing_errors( - padding=20, reference_segment=ref, quality_scores=quals, rng=rng + padding=20, reference_segment=ref, quality_scores=quals, num_errors=3, rng=rng ) # At least one error is expected for these settings. diff --git a/tests/test_models/test_error_models.py b/tests/test_models/test_error_models.py index 5b31708b..3ba927f6 100644 --- a/tests/test_models/test_error_models.py +++ b/tests/test_models/test_error_models.py @@ -149,7 +149,7 @@ def test_sem_zero_error_rate_returns_empty(): m = SequencingErrorModel(avg_seq_error=0.0) rng = np.random.default_rng(0) quality_scores = np.array([40] * 100) - result = m.get_sequencing_errors(20, _SEQ_RECORD, quality_scores, rng) + result = m.get_sequencing_errors(20, _SEQ, quality_scores, 0, rng) assert result == [] @@ -158,7 +158,7 @@ def test_sem_high_error_rate_returns_errors(): m = SequencingErrorModel(avg_seq_error=0.5) rng = np.random.default_rng(0) quality_scores = np.array([1] * 100) # quality 1 → ~79% error rate - result, _ = m.get_sequencing_errors(20, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(20, _SEQ, quality_scores, 3, rng) assert len(result) > 0 @@ -166,7 +166,7 @@ def test_sem_returns_error_container_objects(): m = SequencingErrorModel(avg_seq_error=0.5) rng = np.random.default_rng(0) quality_scores = np.array([1] * 100) - result, _ = m.get_sequencing_errors(20, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(20, _SEQ, quality_scores, 3, rng) for err in result: assert isinstance(err, ErrorContainer) @@ -175,7 +175,7 @@ def test_sem_errors_have_valid_locations(): m = SequencingErrorModel(avg_seq_error=0.5) rng = np.random.default_rng(0) quality_scores = np.array([1] * 100) - result, _ = m.get_sequencing_errors(20, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(20, _SEQ, quality_scores, 3, rng) for err in result: assert 0 <= err.location < len(quality_scores) @@ -184,7 +184,7 @@ def test_sem_snv_errors_have_valid_alt(): m = SequencingErrorModel(avg_seq_error=0.5) rng = np.random.default_rng(0) quality_scores = np.array([1] * 100) - result, _ = m.get_sequencing_errors(20, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(20, _SEQ, quality_scores, 3, rng) snv_errors = [e for e in result if e.error_type == SingleNucleotideVariant] for err in snv_errors: assert err.alt in ("A", "C", "G", "T") @@ -194,7 +194,7 @@ def test_sem_returns_updated_padding(): m = SequencingErrorModel(avg_seq_error=0.5) rng = np.random.default_rng(0) quality_scores = np.array([1] * 100) - _, padding = m.get_sequencing_errors(20, _SEQ_RECORD, quality_scores, rng) + _, padding = m.get_sequencing_errors(20, _SEQ, quality_scores, 3, rng) assert padding >= 0 @@ -202,7 +202,7 @@ def test_sem_high_quality_scores_produce_few_errors(): m = SequencingErrorModel(avg_seq_error=0.009) rng = np.random.default_rng(0) quality_scores = np.array([40] * 151) # q40 → 0.01% error rate - result, _ = m.get_sequencing_errors(20, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(20, _SEQ, quality_scores, 3, rng) # With q40 and length 151, very few errors expected assert len(result) < 10 @@ -250,7 +250,7 @@ def test_sem_deletion_variant_prob_has_no_effect(): ) rng = np.random.default_rng(0) quality_scores = np.array([1] * 151) - result, _ = m.get_sequencing_errors(50, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 3, rng) del_errors = [e for e in result if e.error_type == Del] # No deletions produced due to the unreachable gate (see dead-code comment below) assert len(del_errors) == 0 @@ -265,7 +265,7 @@ def test_sem_insertion_variant_prob_has_no_effect(): ) rng = np.random.default_rng(7) quality_scores = np.array([1] * 151) - result, _ = m.get_sequencing_errors(50, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 3, rng) ins_errors = [e for e in result if e.error_type == Ins] # No insertions produced due to the unreachable gate assert len(ins_errors) == 0 @@ -280,7 +280,7 @@ def test_sem_blacklist_prevents_duplicate_deletion_sites(): ) rng = np.random.default_rng(13) quality_scores = np.array([1] * 151) - result, _ = m.get_sequencing_errors(50, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 3, rng) # All returned errors should be at unique locations (blacklist applied) locations = [e.location for e in result] # Deletions span multiple bases; no two errors at same index @@ -330,7 +330,7 @@ def test_sem_only_snv_errors_produced_regardless_of_variant_probs(): ) rng = np.random.default_rng(0) quality_scores = np.array([1] * 151) - result, _ = m.get_sequencing_errors(50, _SEQ_RECORD, quality_scores, rng) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 3, rng) # All errors are SNVs because the indel branches are unreachable error_types = {e.error_type for e in result} assert error_types == {SingleNucleotideVariant} diff --git a/tests/test_models/test_seq_error.py b/tests/test_models/test_seq_error.py index 041b1aee..6ace4201 100644 --- a/tests/test_models/test_seq_error.py +++ b/tests/test_models/test_seq_error.py @@ -4,7 +4,6 @@ import numpy as np -from Bio.SeqRecord import SeqRecord from Bio.Seq import Seq from neat.models import SequencingErrorModel @@ -17,8 +16,8 @@ def test_get_seq_error_snv_only(): quality_scores = np.full_like(np.arange(10), 36) # Make at least one base very error-prone quality_scores[0] = 0 - reference = SeqRecord(Seq('AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'), id="fake_1", name="fake", description="fake") - errors, _padding = model.get_sequencing_errors(padding=5, reference_segment=reference, quality_scores=quality_scores, rng=rng) + reference = Seq('AAAAAAAAAAAAAAAAAAAAAAAAAAAAAA') + errors, _padding = model.get_sequencing_errors(padding=5, reference_segment=reference, quality_scores=quality_scores, num_errors=3, rng=rng) assert all(e.error_type == SingleNucleotideVariant for e in errors) diff --git a/tests/test_read_simulator/test_generate_reads.py b/tests/test_read_simulator/test_generate_reads.py index 709a3328..184a1b6e 100644 --- a/tests/test_read_simulator/test_generate_reads.py +++ b/tests/test_read_simulator/test_generate_reads.py @@ -401,7 +401,7 @@ def test_generate_reads_single_ended_returns_read_none_pairs(): opts = _make_options(paired=False) cv = ContigVariants() - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, _all_span_targeted(), _nothing_discarded(), opts, None, "chr1", 0, 0) @@ -418,7 +418,7 @@ def test_generate_reads_paired_ended_returns_read_read_pairs(): opts = _make_options(paired=True) cv = ContigVariants() - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, _all_span_targeted(), _nothing_discarded(), opts, None, "chr1", 0, 0) @@ -434,7 +434,7 @@ def test_generate_reads_read_length_matches_options(): opts = _make_options(paired=False) cv = ContigVariants() - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, _all_span_targeted(), _nothing_discarded(), opts, None, "chr1", 0, 0) @@ -454,7 +454,7 @@ def test_generate_reads_targeted_region_flag_false_filters_all(): cv = ContigVariants() no_target = [(0, _SPAN, False)] - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, no_target, _nothing_discarded(), opts, None, "chr1", 0, 0) @@ -469,7 +469,7 @@ def test_generate_reads_discard_region_removes_all(): cv = ContigVariants() discard_all = [(0, _SPAN, True)] - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, _all_span_targeted(), discard_all, opts, None, "chr1", 0, 0) @@ -483,7 +483,7 @@ def test_generate_reads_discard_flag_false_keeps_reads(): opts = _make_options(paired=False) cv = ContigVariants() - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, _all_span_targeted(), _nothing_discarded(), opts, None, "chr1", 0, 0) @@ -509,7 +509,7 @@ def test_generate_reads_variants_populated_on_reads(): ) cv.add_variant(snv) - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, _all_span_targeted(), _nothing_discarded(), opts, None, "chr1", 0, 0) @@ -535,7 +535,7 @@ def test_generate_reads_paired_discard_region_removes_all(): cv = ContigVariants() discard_all = [(0, _SPAN, True)] - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, _all_span_targeted(), discard_all, opts, None, "chr1", 0, 0) @@ -549,7 +549,7 @@ def test_generate_reads_paired_no_discard_produces_read_pairs(): opts = _make_options(paired=True) cv = ContigVariants() - results = generate_reads(0, ref, err, qual, frag, cv, + results = generate_reads(0, ref, err, 3, qual, frag, cv, _all_span_targeted(), _nothing_discarded(), opts, None, "chr1", 0, 0) diff --git a/tests/test_read_simulator/test_read.py b/tests/test_read_simulator/test_read.py index 71d92210..b9d47071 100644 --- a/tests/test_read_simulator/test_read.py +++ b/tests/test_read_simulator/test_read.py @@ -407,7 +407,7 @@ def test_finalize_read_and_write_writes_fastq(): rng = _make_rng() handle = io.StringIO() - r.finalize_read_and_write(err_model, qual_model, handle, 33, True, rng) + r.finalize_read_and_write(err_model, qual_model, handle, 33, True, 3, rng) output = handle.getvalue() assert output.startswith("@test_read") @@ -424,7 +424,7 @@ def test_finalize_read_and_write_reverse_complement(): qual_model = TraditionalQualityModel() rng = _make_rng() - r.finalize_read_and_write(err_model, qual_model, None, 33, False, rng) + r.finalize_read_and_write(err_model, qual_model, None, 33, False, 3, rng) assert len(r.read_sequence) == _READ_LEN @@ -435,7 +435,7 @@ def test_finalize_sets_mapping_quality(): qual_model = TraditionalQualityModel() rng = _make_rng() - r.finalize_read_and_write(err_model, qual_model, None, 33, False, rng) + r.finalize_read_and_write(err_model, qual_model, None, 33, False, 3, rng) assert r.mapping_quality == 70 @@ -450,7 +450,7 @@ def test_make_cigar_all_match(): err_model = SequencingErrorModel(read_length=_READ_LEN) qual_model = TraditionalQualityModel() rng = _make_rng(seed=0) - r.finalize_read_and_write(err_model, qual_model, None, 33, False, rng) + r.finalize_read_and_write(err_model, qual_model, None, 33, False, 3, rng) cigar = r.make_cigar() assert cigar.endswith("M") assert "I" not in cigar or "D" not in cigar # no complex indels for a clean read @@ -462,7 +462,7 @@ def test_make_cigar_reverse_strand(): err_model = SequencingErrorModel(read_length=_READ_LEN) qual_model = TraditionalQualityModel() rng = _make_rng(seed=0) - r.finalize_read_and_write(err_model, qual_model, None, 33, False, rng) + r.finalize_read_and_write(err_model, qual_model, None, 33, False, 3, rng) cigar = r.make_cigar() assert isinstance(cigar, str) assert len(cigar) > 0 \ No newline at end of file diff --git a/tests/test_read_simulator/test_single_runner.py b/tests/test_read_simulator/test_single_runner.py index bc70d2a6..9036700b 100644 --- a/tests/test_read_simulator/test_single_runner.py +++ b/tests/test_read_simulator/test_single_runner.py @@ -304,6 +304,7 @@ def _run(self, tmp_path: Path, *, coverage: int = 2, read_len: int = 50, target_regions, discard_regions, mutation_regions, + 3, ) def test_returns_four_element_tuple(self, tmp_path): @@ -381,6 +382,7 @@ def _run_in(p): return read_simulator_single( 1, 0, opts, None, "chr1", 0, ContigVariants(), [(0, 400, True)], [(0, 400, False)], [(0, 400, 0.01)], + 3, ) _run_in(tmp_a) @@ -465,6 +467,7 @@ def test_record_too_small_logs_and_continues(self, tmp_path, caplog, monkeypatch result = read_simulator_single( 1, 0, opts, None, "chr1", 0, ContigVariants(), [(0, 20, True)], [(0, 20, False)], [(0, 20, 0.01)], + 3, ) assert "Record too small" in caplog.text From c742df45848fd1792a10cb02d2760b977a73ac66 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Sun, 3 May 2026 09:47:26 -0500 Subject: [PATCH 6/7] Add tests for num_errors behaviour and fix fallback infinite-loop Fix fallback loop in get_sequencing_errors: use <= median so uniform quality arrays (all scores equal) always make progress. New tests: - test_sem_num_errors_zero_produces_no_errors - test_sem_num_errors_caps_output - test_sem_fallback_loop_fills_errors_with_uniform_high_quality - test_finalize_read_and_write_returns_error_count - test_errors_per_contig_proportional_to_contig_length - test_errors_per_contig_zero_for_zero_coverage Co-Authored-By: Claude Sonnet 4.6 --- neat/models/error_models.py | 5 ++- tests/test_models/test_error_models.py | 33 ++++++++++++++++ tests/test_read_simulator/test_read.py | 13 +++++++ tests/test_read_simulator/test_runner.py | 39 +++++++++++++++++++ .../test_read_simulator/test_single_runner.py | 2 + 5 files changed, 90 insertions(+), 2 deletions(-) diff --git a/neat/models/error_models.py b/neat/models/error_models.py index 2eb23082..cff8a5d8 100644 --- a/neat/models/error_models.py +++ b/neat/models/error_models.py @@ -199,11 +199,12 @@ def get_sequencing_errors( error_indexes.append(index) i -= 1 # Fallback: if quality scores are too high to naturally reach num_errors, - # force errors at positions with below-median quality scores + # force errors at positions with at-or-below-median quality scores. + # Using <= so that uniform quality arrays (all scores equal) always make progress. median_score = median(quality_scores) while len(error_indexes) < num_errors: index = rng.integers(len(quality_scores)) - if quality_scores[index] < median_score: + if quality_scores[index] <= median_score: error_indexes.append(index) total_indel_length = 0 diff --git a/tests/test_models/test_error_models.py b/tests/test_models/test_error_models.py index 3ba927f6..830af427 100644 --- a/tests/test_models/test_error_models.py +++ b/tests/test_models/test_error_models.py @@ -145,6 +145,39 @@ def test_sem_custom_error_rate(): # SequencingErrorModel — get_sequencing_errors # =========================================================================== +def test_sem_num_errors_zero_produces_no_errors(): + """num_errors=0 short-circuits the error loop even with a high error rate.""" + m = SequencingErrorModel(avg_seq_error=0.9) + rng = np.random.default_rng(0) + quality_scores = np.array([1] * 151) # score 1 → ~79% error rate per base + result, _ = m.get_sequencing_errors(20, _SEQ, quality_scores, 0, rng) + assert result == [] + + +def test_sem_num_errors_caps_output(): + """Error list length never exceeds num_errors.""" + m = SequencingErrorModel(avg_seq_error=0.9) + rng = np.random.default_rng(0) + quality_scores = np.array([1] * 151) + cap = 3 + result, _ = m.get_sequencing_errors(20, _SEQ, quality_scores, cap, rng) + assert len(result) <= cap + + +def test_sem_fallback_loop_fills_errors_with_uniform_high_quality(): + """Fallback loop reaches num_errors when quality is too high for the main loop. + + With quality score 40 (~0.01% error rate) and only 10 iterations, the main + loop almost certainly collects 0 errors. The fallback must make up the deficit. + Using uniform scores verifies the <= median guard prevents an infinite loop. + """ + m = SequencingErrorModel(avg_seq_error=0.5) + rng = np.random.default_rng(0) + quality_scores = np.array([40] * 10) # uniform high quality, 10 iterations max + result, _ = m.get_sequencing_errors(5, _SEQ[:10], quality_scores, 5, rng) + assert len(result) == 5 + + def test_sem_zero_error_rate_returns_empty(): m = SequencingErrorModel(avg_seq_error=0.0) rng = np.random.default_rng(0) diff --git a/tests/test_read_simulator/test_read.py b/tests/test_read_simulator/test_read.py index b9d47071..0b25a884 100644 --- a/tests/test_read_simulator/test_read.py +++ b/tests/test_read_simulator/test_read.py @@ -400,6 +400,19 @@ def test_convert_masking_replaces_ns(): # finalize_read_and_write — produce_fastq=True # --------------------------------------------------------------------------- +def test_finalize_read_and_write_returns_error_count(): + """Return value equals the number of errors actually applied to the read.""" + r = _make_read(reference=_PADDED_REF, padding=20) + err_model = SequencingErrorModel(read_length=_READ_LEN) + qual_model = TraditionalQualityModel() + rng = _make_rng() + + error_count = r.finalize_read_and_write(err_model, qual_model, None, 33, False, 3, rng) + + assert isinstance(error_count, int) + assert error_count == len(r.errors) + + def test_finalize_read_and_write_writes_fastq(): r = _make_read(reference=_PADDED_REF, padding=20) err_model = SequencingErrorModel(read_length=_READ_LEN) diff --git a/tests/test_read_simulator/test_runner.py b/tests/test_read_simulator/test_runner.py index 68a8e8dd..5a00b792 100644 --- a/tests/test_read_simulator/test_runner.py +++ b/tests/test_read_simulator/test_runner.py @@ -141,6 +141,45 @@ def test_filter_bed_regions_returns_list(): assert result == [(0, 200, 0.01)] # only the overlapping region +# =========================================================================== +# errors_per_contig distribution +# =========================================================================== + +def test_errors_per_contig_proportional_to_contig_length(): + """Longer contigs receive proportionally more errors and values sum to >= total_errors.""" + from math import ceil + reference_keys_with_lens = {"chr1": 1000, "chr2": 3000} + average_error = 0.01 + coverage = 10 + + total_reference_length = sum(reference_keys_with_lens.values()) + total_errors = ceil(average_error * total_reference_length * coverage) + normalized_counts = {k: v / total_reference_length for k, v in reference_keys_with_lens.items()} + errors_per_contig = {k: ceil(v * total_errors) for k, v in normalized_counts.items()} + + # chr2 is 3× longer so should receive more errors + assert errors_per_contig["chr2"] > errors_per_contig["chr1"] + # ceil can cause slight overshoot, but sum should be within one per contig of total + assert sum(errors_per_contig.values()) >= total_errors + assert sum(errors_per_contig.values()) <= total_errors + len(reference_keys_with_lens) + + +def test_errors_per_contig_zero_for_zero_coverage(): + """Zero coverage produces zero total errors and zero per contig.""" + from math import ceil + reference_keys_with_lens = {"chr1": 1000, "chr2": 500} + average_error = 0.01 + coverage = 0 + + total_reference_length = sum(reference_keys_with_lens.values()) + total_errors = ceil(average_error * total_reference_length * coverage) + normalized_counts = {k: v / total_reference_length for k, v in reference_keys_with_lens.items()} + errors_per_contig = {k: ceil(v * total_errors) for k, v in normalized_counts.items()} + + assert total_errors == 0 + assert all(v == 0 for v in errors_per_contig.values()) + + # =========================================================================== # Integration test — read_simulator_runner (FASTQ output only) # =========================================================================== diff --git a/tests/test_read_simulator/test_single_runner.py b/tests/test_read_simulator/test_single_runner.py index 9036700b..8b220566 100644 --- a/tests/test_read_simulator/test_single_runner.py +++ b/tests/test_read_simulator/test_single_runner.py @@ -417,6 +417,7 @@ def test_bam_output_written(self, tmp_path): _, _, _, file_dict = read_simulator_single( 1, 0, opts, bam_header, "chr1", 0, ContigVariants(), [(0, 400, True)], [(0, 400, False)], [(0, 400, 0.01)], + 3, ) assert file_dict["bam"] == opts.bam @@ -434,6 +435,7 @@ def test_bam_key_is_none_when_not_requested(self, tmp_path): _, _, _, file_dict = read_simulator_single( 1, 0, opts, None, "chr1", 0, ContigVariants(), [(0, 400, True)], [(0, 400, False)], [(0, 400, 0.01)], + 3, ) assert file_dict["bam"] is None From bcbf78ca57e0c14c6cae73f5b71702e4832ef474 Mon Sep 17 00:00:00 2001 From: Joshua Allen Date: Sun, 3 May 2026 10:03:49 -0500 Subject: [PATCH 7/7] Fix deletion blacklist, indel gate, and quality array dtype bugs - Fix blacklist range to exclude the deletion anchor so deletions don't remove themselves during the cleanup pass - Fix quality array concatenation to use dtype=int so deletion errors (which produce an empty quality slice) don't promote the array to float64 and break the chr() call in finalize_read_and_write - Remove dead num_errors capture from finalize_read_and_write return value - Add tests for indel error paths, errors_per_read gate logic, and score clamping in TraditionalQualityModel Co-Authored-By: Claude Sonnet 4.6 --- neat/models/error_models.py | 6 +- neat/read_simulator/utils/generate_reads.py | 4 +- neat/read_simulator/utils/read.py | 4 +- tests/test_models/test_error_models.py | 96 +++++++++------------ tests/test_read_simulator/test_runner.py | 63 ++++++++++++++ 5 files changed, 113 insertions(+), 60 deletions(-) diff --git a/neat/models/error_models.py b/neat/models/error_models.py index cff8a5d8..f7036337 100644 --- a/neat/models/error_models.py +++ b/neat/models/error_models.py @@ -219,7 +219,7 @@ def get_sequencing_errors( # This is to prevent deletion error collisions and to keep there from being too many indel errors. if 0 < index < self.read_length - max( - self.deletion_len_model) and total_indel_length > self.read_length // 4: + self.deletion_len_model) and total_indel_length <= self.read_length // 4: error_type = rng.choice(a=list(self.variant_probs), p=list(self.variant_probs.values())) # Deletion error @@ -235,11 +235,11 @@ def get_sequencing_errors( ) total_indel_length += deletion_length - del_blacklist.extend(list(range(index, index + deletion_length))) + del_blacklist.extend(list(range(index + 1, index + deletion_length + 1))) padding -= deletion_length elif error_type == Insertion: - insertion_length = self.get_insertion_length() + insertion_length = self.get_insertion_length(rng) insertion_reference = reference_segment[index] insert_string = ''.join(rng.choice(ALLOWED_NUCL, size=insertion_length)) insertion_alternate = insertion_reference + insert_string diff --git a/neat/read_simulator/utils/generate_reads.py b/neat/read_simulator/utils/generate_reads.py index 6b5cbfba..ce15b990 100644 --- a/neat/read_simulator/utils/generate_reads.py +++ b/neat/read_simulator/utils/generate_reads.py @@ -281,7 +281,7 @@ def generate_reads( fastq_handle = ofw.files_to_write[ofw.fq1] else: fastq_handle = None - num_errors = read_1.finalize_read_and_write( + read_1.finalize_read_and_write( error_model, qual_model, fastq_handle, @@ -319,7 +319,7 @@ def generate_reads( fastq_handle = ofw.files_to_write[ofw.fq2] else: fastq_handle = None - num_errors = read_2.finalize_read_and_write( + read_2.finalize_read_and_write( error_model, qual_model, fastq_handle, diff --git a/neat/read_simulator/utils/read.py b/neat/read_simulator/utils/read.py index cd14ea57..94c922a0 100644 --- a/neat/read_simulator/utils/read.py +++ b/neat/read_simulator/utils/read.py @@ -71,7 +71,7 @@ def __init__(self, self.read_sequence: Seq = Seq("") # initialize to empty sequence self.errors: list[ErrorContainer] = [] # initialize self.mutations: dict[int, list] = {} # initialize - self.quality_array: np.ndarray = np.zeros(self.run_read_length) # this will have the correct memory length + self.quality_array: np.ndarray = np.zeros(self.run_read_length, dtype=int) # this will have the correct memory length self.mapping_quality: int = 0 # initialize at 0 self.read_quality_string: str = "" # This will hold the read quality string self.num_ns = 0 @@ -161,7 +161,7 @@ def update_quality_array( # Replace the given quality score with the new one self.quality_array = \ np.concatenate((self.quality_array[:location], - np.array(new_quality_score), + np.array(new_quality_score, dtype=int), self.quality_array[location+ref_length:])) def apply_errors(self, quality_model: TraditionalQualityModel): diff --git a/tests/test_models/test_error_models.py b/tests/test_models/test_error_models.py index 830af427..0e6039ef 100644 --- a/tests/test_models/test_error_models.py +++ b/tests/test_models/test_error_models.py @@ -266,15 +266,14 @@ def test_error_container_insertion_type(): # =========================================================================== -# SequencingErrorModel — indel error paths (lines 209, 213-235) +# SequencingErrorModel — indel error paths # =========================================================================== -def test_sem_deletion_variant_prob_has_no_effect(): - """variant_probs favouring Deletion still produces only SNVs (dead-code bug). +def test_sem_deletion_errors_produced_with_deletion_variant_probs(): + """variant_probs favouring Deletion should produce deletion errors. - The indel gate condition is circular, so deletion errors are never produced - regardless of variant_probs. See test_sem_only_snv_errors_produced_regardless_of_variant_probs - for the full documentation test. + The gate condition `total_indel_length <= read_length // 4` allows indels + until they fill a quarter of the read, then switches back to SNVs. """ from neat.variants import Deletion as Del, Insertion as Ins m = SequencingErrorModel( @@ -283,14 +282,13 @@ def test_sem_deletion_variant_prob_has_no_effect(): ) rng = np.random.default_rng(0) quality_scores = np.array([1] * 151) - result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 3, rng) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 5, rng) del_errors = [e for e in result if e.error_type == Del] - # No deletions produced due to the unreachable gate (see dead-code comment below) - assert len(del_errors) == 0 + assert len(del_errors) > 0, "Expected deletion errors given variant_probs favours them" -def test_sem_insertion_variant_prob_has_no_effect(): - """variant_probs favouring Insertion still produces only SNVs (dead-code bug).""" +def test_sem_insertion_errors_produced_with_insertion_variant_probs(): + """variant_probs favouring Insertion should produce insertion errors.""" from neat.variants import Insertion as Ins, Deletion as Del m = SequencingErrorModel( avg_seq_error=0.9, @@ -298,14 +296,30 @@ def test_sem_insertion_variant_prob_has_no_effect(): ) rng = np.random.default_rng(7) quality_scores = np.array([1] * 151) - result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 3, rng) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 5, rng) ins_errors = [e for e in result if e.error_type == Ins] - # No insertions produced due to the unreachable gate - assert len(ins_errors) == 0 + assert len(ins_errors) > 0, "Expected insertion errors given variant_probs favours them" + + +def test_sem_indel_cap_limits_total_indel_length(): + """Total indel length in errors should not exceed read_length // 4.""" + from neat.variants import Deletion as Del, Insertion as Ins + read_length = 151 + m = SequencingErrorModel( + avg_seq_error=0.9, + read_length=read_length, + variant_probs={Ins: 0.0, Del: 1.0, SingleNucleotideVariant: 0.0}, + ) + rng = np.random.default_rng(0) + quality_scores = np.array([1] * read_length) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 20, rng) + del_errors = [e for e in result if e.error_type == Del] + total_indel_length = sum(e.length for e in del_errors) + assert total_indel_length <= read_length // 4 def test_sem_blacklist_prevents_duplicate_deletion_sites(): - """Errors at blacklisted positions from a deletion are removed.""" + """Errors at positions spanned by a deletion are removed via the blacklist.""" from neat.variants import Deletion as Del, Insertion as Ins m = SequencingErrorModel( avg_seq_error=0.9, @@ -313,13 +327,24 @@ def test_sem_blacklist_prevents_duplicate_deletion_sites(): ) rng = np.random.default_rng(13) quality_scores = np.array([1] * 151) - result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 3, rng) - # All returned errors should be at unique locations (blacklist applied) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 5, rng) locations = [e.location for e in result] - # Deletions span multiple bases; no two errors at same index assert len(locations) == len(set(locations)) +def test_sem_snv_only_probs_produces_no_indels(): + """variant_probs with SNV=1.0 should produce zero indel errors.""" + from neat.variants import Deletion as Del, Insertion as Ins + m = SequencingErrorModel( + avg_seq_error=0.9, + variant_probs={Ins: 0.0, Del: 0.0, SingleNucleotideVariant: 1.0}, + ) + rng = np.random.default_rng(0) + quality_scores = np.array([1] * 151) + result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 5, rng) + indel_errors = [e for e in result if e.error_type in (Del, Ins)] + assert len(indel_errors) == 0 + # =========================================================================== # TraditionalQualityModel — score clamping (line 104) @@ -345,38 +370,3 @@ def test_tqm_score_clamped_to_maximum(): assert all(s == 42 for s in scores) -# =========================================================================== -# SequencingErrorModel — dead-code documentation -# =========================================================================== -# Lines 209-235, 252 (indel error branches) are unreachable because -# `total_indel_length` starts at 0 and is only incremented inside the -# branches that are gated by `total_indel_length > self.read_length // 4`. -# This circular dependency means the variant_probs choice (line 209) is -# never called and deletion/insertion errors are never produced. -# The following test documents this behaviour. - -def test_sem_only_snv_errors_produced_regardless_of_variant_probs(): - """Indel errors are never produced due to the total_indel_length gate.""" - m = SequencingErrorModel( - avg_seq_error=0.9, - variant_probs={Insertion: 0.5, Deletion: 0.5, SingleNucleotideVariant: 0.0}, - ) - rng = np.random.default_rng(0) - quality_scores = np.array([1] * 151) - result, _ = m.get_sequencing_errors(50, _SEQ, quality_scores, 3, rng) - # All errors are SNVs because the indel branches are unreachable - error_types = {e.error_type for e in result} - assert error_types == {SingleNucleotideVariant} - # TODO (post-fix): Once the gate condition is corrected from - # `total_indel_length > self.read_length // 4` - # to - # `total_indel_length <= self.read_length // 4` - # the following assertions should replace the one above: - # - # del_errors = [e for e in result if e.error_type == Deletion] - # ins_errors = [e for e in result if e.error_type == Insertion] - # assert len(del_errors) + len(ins_errors) > 0, \ - # "Expected indel errors given variant_probs favours them" - # # blacklist test: no two errors share a location - # locations = [e.location for e in result] - # assert len(locations) == len(set(locations)) diff --git a/tests/test_read_simulator/test_runner.py b/tests/test_read_simulator/test_runner.py index 5a00b792..d01c71d4 100644 --- a/tests/test_read_simulator/test_runner.py +++ b/tests/test_read_simulator/test_runner.py @@ -180,6 +180,69 @@ def test_errors_per_contig_zero_for_zero_coverage(): assert all(v == 0 for v in errors_per_contig.values()) +def test_errors_per_read_fractional_gate_fires_probabilistically(): + """When block_errors rounds to 0 but is non-zero, the gate may increment errors_per_read. + + Exercises runner.py: + errors_per_read = round(block_errors / estimated_number_of_reads) + if errors_per_read < 1.0 and block_errors > 0: + if rng.random() < average_error: + errors_per_read += 1 + """ + import numpy as np + + block_errors = 0.3 # rounds to 0 + estimated_number_of_reads = 1 + average_error = 0.9 # high rate → gate almost always fires + + errors_per_read = round(block_errors / estimated_number_of_reads) + assert errors_per_read == 0 + + rng = np.random.default_rng(0) + if errors_per_read < 1.0 and block_errors > 0: + if rng.random() < average_error: + errors_per_read += 1 + + assert errors_per_read == 1 # gate fired with high average_error and seed 0 + + +def test_errors_per_read_gate_skipped_when_block_errors_zero(): + """Gate is not entered when block_errors == 0, leaving errors_per_read at 0.""" + import numpy as np + + block_errors = 0.0 + estimated_number_of_reads = 10 + average_error = 0.9 + + errors_per_read = round(block_errors / estimated_number_of_reads) + rng = np.random.default_rng(0) + if errors_per_read < 1.0 and block_errors > 0: + if rng.random() < average_error: + errors_per_read += 1 + + assert errors_per_read == 0 + + +def test_errors_per_read_gate_skipped_when_already_positive(): + """Gate is not entered when errors_per_read rounds to >= 1.""" + import numpy as np + + block_errors = 10.0 + estimated_number_of_reads = 2 # → round(5.0) = 5 + average_error = 0.9 + + errors_per_read = round(block_errors / estimated_number_of_reads) + assert errors_per_read >= 1 + + original = errors_per_read + rng = np.random.default_rng(0) + if errors_per_read < 1.0 and block_errors > 0: + if rng.random() < average_error: + errors_per_read += 1 + + assert errors_per_read == original # unchanged + + # =========================================================================== # Integration test — read_simulator_runner (FASTQ output only) # ===========================================================================