From ce600263b0afe6a955ae98a7d4a81530a63abd50 Mon Sep 17 00:00:00 2001 From: Keshav Date: Sun, 26 Apr 2026 15:42:02 +0200 Subject: [PATCH] Solving VCF-related issues. --- neat/models/mutation_model.py | 14 +- neat/read_simulator/utils/stitch_outputs.py | 9 + neat/read_simulator/utils/vcf_func.py | 7 + neat/variants/contig_variants.py | 19 +- .../test_generate_variants.py | 35 ++-- .../test_stitch_outputs.py | 76 ++++++-- tests/test_read_simulator/test_vcf_func.py | 108 ++--------- tests/test_variants/test_contig_variants.py | 20 --- tests/test_variants/test_ref_eq_alt.py | 170 ++++++++++++++++++ 9 files changed, 310 insertions(+), 148 deletions(-) create mode 100644 tests/test_variants/test_ref_eq_alt.py diff --git a/neat/models/mutation_model.py b/neat/models/mutation_model.py index b324448d..9c7a3074 100644 --- a/neat/models/mutation_model.py +++ b/neat/models/mutation_model.py @@ -108,9 +108,19 @@ def generate_snv(self, trinucleotide: Seq, reference_location: int, rng: Generat # First determine which matrix to use transition_matrix = self.trinuc_trans_matrices[DINUC_IND[trinucleotide[0] + "_" + trinucleotide[2]]] # then determine the trans probs based on the middle nucleotide - transition_probs = transition_matrix[NUC_IND[trinucleotide[1]]] - # Creating probabilities from the weights + transition_probs = list(transition_matrix[NUC_IND[trinucleotide[1]]]) + # Zero the ref-base probability so that we never pick REF==ALT, even with custom models with trans matrices + # that have non-zero diagonal entries (edge case) + ref_base_idx = NUC_IND[trinucleotide[1]] + transition_probs[ref_base_idx] = 0.0 transition_sum = sum(transition_probs) + if transition_sum == 0.0: + _LOG.warning( + f"Transition matrix row for '{trinucleotide[1]}' has all weight on the reference base. " + f"Falling back to uniform sampling of non-ref bases." + ) + transition_probs = [1.0 if i != ref_base_idx else 0.0 for i in range(len(ALLOWED_NUCL))] + transition_sum = sum(transition_probs) transition_probs = [x/transition_sum for x in transition_probs] # Now pick a random alternate, weighted by the probabilities alt = rng.choice(ALLOWED_NUCL, p=transition_probs) diff --git a/neat/read_simulator/utils/stitch_outputs.py b/neat/read_simulator/utils/stitch_outputs.py index 101c50a4..49bd8ac0 100644 --- a/neat/read_simulator/utils/stitch_outputs.py +++ b/neat/read_simulator/utils/stitch_outputs.py @@ -24,11 +24,20 @@ def concat(files_to_join: List[Path], dest_file: gzip.GzipFile) -> None: def merge_vcfs(vcfs: List[Path], ofw: OutputFileWriter) -> None: dest = ofw.files_to_write[ofw.vcf] + seen: set[str] = set() + n_duplicates = 0 for vcf in vcfs: with gzip.open(vcf, 'rt') as fh: for line in fh: if not line.startswith("#"): + normalized = line.rstrip("\r\n") + if normalized in seen: + n_duplicates += 1 + continue + seen.add(normalized) dest.write(line) + if n_duplicates: + _LOG.warning(f"merge_vcfs: removed {n_duplicates} duplicate VCF line(s) during merge.") def merge_bam(bam_files: List[Path], ofw: OutputFileWriter, threads: int): merged_file = ofw.tmp_dir / "temp_merged.bam" diff --git a/neat/read_simulator/utils/vcf_func.py b/neat/read_simulator/utils/vcf_func.py index a8d86334..584550f0 100755 --- a/neat/read_simulator/utils/vcf_func.py +++ b/neat/read_simulator/utils/vcf_func.py @@ -228,6 +228,13 @@ def parse_input_vcf( count = 0 for alt in alts: count += 1 + if ref == alt: + _LOG.warning( + f"Skipping variant at {chrom}:{location + 1} — REF == ALT ({ref!r}). " + f"This is not a valid variant." + ) + n_skipped += 1 + continue # This temp genotype teases out only the ploids with this particular variant temp_genotype = variant_genotype(options.ploidy, genotype, count) if len(ref) == len(alt) == 1: diff --git a/neat/variants/contig_variants.py b/neat/variants/contig_variants.py index 18267c43..72f5f09b 100644 --- a/neat/variants/contig_variants.py +++ b/neat/variants/contig_variants.py @@ -94,13 +94,26 @@ def generate_field(self, variant, field): def find_dups(self, variant): """ - Checks if the given genotype is already present in a list of variants. + Checks if an equivalent variant already exists at this position. + Two variants are duplicates when they share the same type and ALT allele. + Genotype-only comparison was insufficient: two variants with identical ALT + but different genotypes would produce two identical VCF output lines. :param variant: A variant to check for duplicates - :return: True or False if found or not + :return: True if a duplicate exists, False otherwise """ + try: + variant_alt = variant.get_alt() + except (KeyError, AttributeError): + variant_alt = None + for existing_var in self.contig_variants[variant.position1]: - if np.array_equal(variant.genotype, existing_var.genotype): + try: + existing_alt = existing_var.get_alt() + except (KeyError, AttributeError): + existing_alt = None + + if type(variant) == type(existing_var) and variant_alt == existing_alt: return True return False diff --git a/tests/test_read_simulator/test_generate_variants.py b/tests/test_read_simulator/test_generate_variants.py index 00948c27..6d61f1cc 100644 --- a/tests/test_read_simulator/test_generate_variants.py +++ b/tests/test_read_simulator/test_generate_variants.py @@ -16,9 +16,7 @@ from neat.variants import ContigVariants, SingleNucleotideVariant, Insertion, Deletion -# --------------------------------------------------------------------------- # Helpers -# --------------------------------------------------------------------------- _CLEAN_SEQ = "ACGT" * 50 # 200 bp, no N's _N_HEAVY = "N" * 95 + "ACGT" # 99 bases, >10% N @@ -44,9 +42,7 @@ def _full_rate_regions(seq_len: int, rate: float = 0.01, offset: int = 0): return [(offset, offset + seq_len, rate)] -# =========================================================================== # find_random_non_n -# =========================================================================== def test_find_random_non_n_returns_valid_index(): rng = np.random.default_rng(0) @@ -77,9 +73,7 @@ def test_find_random_non_n_single_element(): assert find_random_non_n(rng, safe_zones) == 0 -# =========================================================================== # map_non_n_regions -# =========================================================================== def test_map_non_n_regions_clean_sequence(): result = map_non_n_regions("ACGTACGTACGT") @@ -88,7 +82,7 @@ def test_map_non_n_regions_clean_sequence(): def test_map_non_n_regions_single_n(): - # 1 N in 50 bases = 2% N → valid map returned + # 1 N in 50 bases = 2% N (valid map returned) seq = "A" * 24 + "N" + "A" * 25 result = map_non_n_regions(seq) assert len(result) == 50 @@ -106,7 +100,7 @@ def test_map_non_n_regions_too_many_ns_returns_empty(): def test_map_non_n_regions_exactly_at_threshold(): """Exactly 10% N → should return empty (condition is <= 0.90 non-N).""" - # 10 N's + 90 ACGT → 90% non-N, average == 0.90, which hits the <= boundary + # 10 N's + 90 ACGT (90% non-N, average == 0.90, which hits the <= boundary) seq = "N" * 10 + "A" * 90 result = map_non_n_regions(seq) assert len(result) == 0 @@ -118,7 +112,7 @@ def test_map_non_n_regions_all_n_returns_empty(): def test_map_non_n_regions_run_of_ns(): - seq = "ACGT" + "NNNN" + "ACGT" # 12 bp, 4/12 ≈ 33% N → empty + seq = "ACGT" + "NNNN" + "ACGT" # 12 bp, 4/12 ≈ 33% N (empty) result = map_non_n_regions(seq) assert len(result) == 0 @@ -133,9 +127,7 @@ def test_map_non_n_regions_short_n_run_in_long_sequence(): assert result[0] == 1 -# =========================================================================== # generate_variants — input variants are copied into output -# =========================================================================== def test_generate_variants_returns_contig_variants(): ref = _make_reference() @@ -206,9 +198,7 @@ def test_generate_variants_input_variant_before_offset_excluded(): assert 50 not in result -# =========================================================================== # generate_variants — random mutation generation -# =========================================================================== def test_generate_variants_adds_at_least_min_mutations(): """With min_mutations=1, at least 1 variant should be added.""" @@ -339,9 +329,7 @@ def test_generate_variants_input_and_random_together(): assert len(result.variant_locations) > 1 -# =========================================================================== # generate_variants — N-handling paths (lines 139-157, 204) -# =========================================================================== def test_generate_variants_n_in_mutation_region_completes(): """Sequence with N's in the mutation region runs to completion. @@ -370,7 +358,7 @@ def test_generate_variants_n_heavy_subsequence_skipped(): ref = _make_reference(seq) model = _make_model() opts = _make_options(seed=3) - opts.min_mutations = 0 # let poisson decide; main goal is no crash + opts.min_mutations = 0 # let Poisson decide; main goal is no crash result = generate_variants(ref, 0, _full_rate_regions(len(seq)), ContigVariants(), model, opts, 40) assert isinstance(result, ContigVariants) @@ -406,9 +394,16 @@ def test_generate_variants_deletion_overlap_handling(): result = generate_variants(ref, 0, _full_rate_regions(len(seq), 0.05), ContigVariants(), model_high, opts, 40) assert isinstance(result, ContigVariants) - # Variants at the same location are deduplicated correctly + # Variants at the same location are deduplicated by (type, ALT) — not genotype. + # Two variants of the same type with the same ALT at the same position are duplicates. for loc in result.variant_locations: variants_here = result.contig_variants[loc] - genotypes = [tuple(v.genotype) for v in variants_here] - assert len(genotypes) == len(set(genotypes)), \ - f"Duplicate genotype at location {loc}" + type_alt_keys = [] + for v in variants_here: + try: + alt = v.get_alt() + except (KeyError, AttributeError): + alt = None + type_alt_keys.append((type(v).__name__, alt)) + assert len(type_alt_keys) == len(set(type_alt_keys)), \ + f"Duplicate (type, ALT) pair at location {loc}: {type_alt_keys}" diff --git a/tests/test_read_simulator/test_stitch_outputs.py b/tests/test_read_simulator/test_stitch_outputs.py index 8349ccf3..ed830f85 100644 --- a/tests/test_read_simulator/test_stitch_outputs.py +++ b/tests/test_read_simulator/test_stitch_outputs.py @@ -7,14 +7,15 @@ from types import SimpleNamespace from unittest.mock import MagicMock, patch, call +import numpy as np import pytest from neat.read_simulator.utils.stitch_outputs import concat, merge_vcfs, merge_bam, main +from neat.variants import SingleNucleotideVariant +from neat.variants.contig_variants import ContigVariants -# --------------------------------------------------------------------------- # Helpers -# --------------------------------------------------------------------------- def _write_gz(path: Path, text: str) -> Path: with gzip.open(path, "wt") as fh: @@ -49,9 +50,7 @@ def _make_ofw(tmp_path: Path, vcf_path: Path = None): return ofw -# =========================================================================== # concat -# =========================================================================== def test_concat_single_file(tmp_path): src = _write_gz(tmp_path / "a.gz", "hello\n") @@ -91,9 +90,7 @@ def test_concat_order_is_preserved(tmp_path): assert positions == sorted(positions) -# =========================================================================== # merge_vcfs -# =========================================================================== def test_merge_vcfs_skips_comment_lines(tmp_path): vcf_text = "##header line\n#CHROM\tPOS\n1\t100\tA\tT\n" @@ -139,9 +136,70 @@ def test_merge_vcfs_preserves_data_line_order(tmp_path): assert positions == list(range(1, 6)) -# =========================================================================== +def test_merge_vcfs_dedup_removes_identical_lines(tmp_path): + """Identical lines from two thread VCFs are collapsed to one (Issue #256).""" + line = "chr1\t100\t.\tA\tT\t42\tPASS\t.\tGT\t0|1\n" + v1 = _write_gz(tmp_path / "t0.vcf.gz", line) + v2 = _write_gz(tmp_path / "t1.vcf.gz", line) + ofw = _make_ofw(tmp_path) + merge_vcfs([v1, v2], ofw) + result = [l for l in ofw._vcf_buf.getvalue().splitlines() if l.strip()] + assert len(result) == 1 + + +def test_merge_vcfs_distinct_lines_are_all_kept(tmp_path): + """Distinct lines from two threads both appear in merged output.""" + v1 = _write_gz(tmp_path / "t0.vcf.gz", "chr1\t100\t.\tA\tT\t42\tPASS\t.\tGT\t0|1\n") + v2 = _write_gz(tmp_path / "t1.vcf.gz", "chr1\t200\t.\tC\tG\t42\tPASS\t.\tGT\t0|1\n") + ofw = _make_ofw(tmp_path) + merge_vcfs([v1, v2], ofw) + result = [l for l in ofw._vcf_buf.getvalue().splitlines() if l.strip()] + assert len(result) == 2 + + +def test_merge_vcfs_partial_overlap_deduped(tmp_path): + """Three lines total, two of which are identical: result has two unique lines.""" + line_a = "chr1\t100\t.\tA\tT\t42\tPASS\t.\tGT\t0|1\n" + line_b = "chr1\t200\t.\tC\tG\t42\tPASS\t.\tGT\t0|1\n" + v1 = _write_gz(tmp_path / "t0.vcf.gz", line_a + line_b) + v2 = _write_gz(tmp_path / "t1.vcf.gz", line_a) + ofw = _make_ofw(tmp_path) + merge_vcfs([v1, v2], ofw) + result = [l for l in ofw._vcf_buf.getvalue().splitlines() if l.strip()] + assert len(result) == 2 + + +# find_dups (ContigVariants deduplication, Issue #256) + +def test_find_dups_same_alt_different_genotype_rejected(tmp_path): + """Same position + same ALT is a duplicate regardless of genotype.""" + cv = ContigVariants() + v1 = SingleNucleotideVariant(10, "T", np.array([1, 0]), 40) + v2 = SingleNucleotideVariant(10, "T", np.array([0, 1]), 40) + cv.add_variant(v1) + assert cv.add_variant(v2) == 1 + + +def test_find_dups_different_alt_same_position_accepted(tmp_path): + """Two SNVs at the same position with different ALTs are not duplicates.""" + cv = ContigVariants() + v1 = SingleNucleotideVariant(10, "T", np.array([1, 0]), 40) + v2 = SingleNucleotideVariant(10, "G", np.array([0, 1]), 40) + cv.add_variant(v1) + assert cv.add_variant(v2) == 0 + assert len(cv.contig_variants[10]) == 2 + + +def test_find_dups_exact_duplicate_rejected(tmp_path): + """Exact duplicates (same position, ALT, and genotype) are rejected.""" + cv = ContigVariants() + v1 = SingleNucleotideVariant(10, "T", np.array([0, 1]), 40) + v2 = SingleNucleotideVariant(10, "T", np.array([0, 1]), 40) + cv.add_variant(v1) + assert cv.add_variant(v2) == 1 + + # merge_bam -# =========================================================================== def test_merge_bam_calls_pysam_merge_and_sort(tmp_path): ofw = _make_ofw(tmp_path) @@ -191,9 +249,7 @@ def test_merge_bam_chunks_large_bam_list(tmp_path): assert mock_pysam.merge.call_count == 3 -# =========================================================================== # main -# =========================================================================== def _file_dict(fq1=None, fq2=None, vcf=None, bam=None): return {"fq1": fq1, "fq2": fq2, "vcf": vcf, "bam": bam} diff --git a/tests/test_read_simulator/test_vcf_func.py b/tests/test_read_simulator/test_vcf_func.py index 09f2ae3b..8aa9a333 100644 --- a/tests/test_read_simulator/test_vcf_func.py +++ b/tests/test_read_simulator/test_vcf_func.py @@ -21,9 +21,7 @@ from neat.variants.unknown_variant import UnknownVariant -# --------------------------------------------------------------------------- # Shared fixtures and helpers -# --------------------------------------------------------------------------- # Reference sequence: chr1 = ACGTACGTACGTACGTACGT (20 bp) # chr2 = TTGGTTGGTTGG (12 bp) @@ -72,9 +70,7 @@ def _vcf_header_with_format(sample="SAMPLE1"): ] -# =========================================================================== # retrieve_genotype -# =========================================================================== def _make_vcf_record(format_field, sample_field, info="."): """Build a minimal 10-column VCF record list.""" @@ -117,9 +113,7 @@ def test_retrieve_genotype_cancer_uses_column_10(): np.testing.assert_array_equal(gt, [1, 1]) -# =========================================================================== # variant_genotype -# =========================================================================== def test_variant_genotype_no_match(): gt = variant_genotype(2, np.array([0, 0]), 1) @@ -149,9 +143,7 @@ def test_variant_genotype_returns_correct_ploidy_length(): np.testing.assert_array_equal(gt, [1, 0, 1, 0]) -# =========================================================================== # parse_input_vcf — variant type classification -# =========================================================================== def test_parse_snv(tmp_path, ref_fasta, empty_input_dict, opts): """REF and ALT both length 1 → SingleNucleotideVariant.""" @@ -201,9 +193,7 @@ def test_parse_unknown_variant(tmp_path, ref_fasta, empty_input_dict, opts): assert isinstance(variants[0], UnknownVariant) -# =========================================================================== # parse_input_vcf — filtering / skipping -# =========================================================================== def test_chrom_not_in_reference_skipped(tmp_path, ref_fasta, empty_input_dict, opts): """Variants on chromosomes absent from the reference are silently skipped.""" @@ -225,18 +215,27 @@ def test_ref_mismatch_skipped(tmp_path, ref_fasta, empty_input_dict, opts): assert len(empty_input_dict["chr1"].variant_locations) == 0 -def test_duplicate_position_skipped(tmp_path, ref_fasta, empty_input_dict, opts): - """A second variant with the same genotype at the same position is skipped. - add_variant deduplicates by (position, genotype), not just position.""" - # Both records share GT=0|1 → same temp_genotype → second is a dup - vcf = _write_vcf(tmp_path, "dup.vcf", _vcf_header_with_format() + [ +def test_duplicate_position_same_alt_skipped(tmp_path, ref_fasta, empty_input_dict, opts): + """A second variant at the same position with the same ALT is a duplicate and skipped, + regardless of genotype. add_variant deduplicates by (position, type, ALT).""" + vcf = _write_vcf(tmp_path, "dup_same_alt.vcf", _vcf_header_with_format() + [ "chr1\t1\t.\tA\tG\t30\tPASS\t.\tGT\t0|1", - "chr1\t1\t.\tA\tC\t30\tPASS\t.\tGT\t0|1", + "chr1\t1\t.\tA\tG\t30\tPASS\t.\tGT\t1|0", # same ALT 'G', different genotype ]) parse_input_vcf(empty_input_dict, vcf, 2, ref_fasta, opts) assert len(empty_input_dict["chr1"].contig_variants[0]) == 1 +def test_duplicate_position_different_alt_both_accepted(tmp_path, ref_fasta, empty_input_dict, opts): + """Two variants at the same position with different ALTs are not duplicates — both accepted.""" + vcf = _write_vcf(tmp_path, "dup_diff_alt.vcf", _vcf_header_with_format() + [ + "chr1\t1\t.\tA\tG\t30\tPASS\t.\tGT\t0|1", + "chr1\t1\t.\tA\tC\t30\tPASS\t.\tGT\t0|1", # different ALT 'C' + ]) + parse_input_vcf(empty_input_dict, vcf, 2, ref_fasta, opts) + assert len(empty_input_dict["chr1"].contig_variants[0]) == 2 + + def test_comment_and_header_lines_not_parsed_as_variants(tmp_path, ref_fasta, empty_input_dict, opts): """## header lines and #CHROM line are never treated as variant records.""" vcf = _write_vcf(tmp_path, "headers.vcf", [ @@ -249,9 +248,7 @@ def test_comment_and_header_lines_not_parsed_as_variants(tmp_path, ref_fasta, em assert len(empty_input_dict["chr1"].variant_locations) == 1 -# =========================================================================== # parse_input_vcf — QUAL handling -# =========================================================================== def test_missing_qual_replaced_with_42(tmp_path, ref_fasta, empty_input_dict, opts): """QUAL field '.' is replaced with the default value '42'.""" @@ -263,9 +260,7 @@ def test_missing_qual_replaced_with_42(tmp_path, ref_fasta, empty_input_dict, op assert variants[0].qual_score == "42" -# =========================================================================== # parse_input_vcf — FORMAT / genotype handling -# =========================================================================== def test_with_format_gt_uses_sample_genotype(tmp_path, ref_fasta, empty_input_dict, opts): """FORMAT column with GT field reads genotype from the sample column.""" @@ -321,9 +316,7 @@ def test_format_exits_if_no_sample_column(tmp_path, ref_fasta, empty_input_dict, parse_input_vcf(empty_input_dict, vcf, 2, ref_fasta, opts) -# =========================================================================== # parse_input_vcf — multiple ALTs -# =========================================================================== def test_multiple_alts_each_gets_variant(tmp_path, ref_fasta, empty_input_dict, opts): """A comma-separated ALT field produces one variant object per alt allele.""" @@ -338,9 +331,7 @@ def test_multiple_alts_each_gets_variant(tmp_path, ref_fasta, empty_input_dict, assert alts == {"G", "C"} -# =========================================================================== # parse_input_vcf — multiple contigs and is_input flag -# =========================================================================== def test_variants_routed_to_correct_contig(tmp_path, ref_fasta, empty_input_dict, opts): """Variants on different chromosomes end up in the correct ContigVariants.""" @@ -363,72 +354,3 @@ def test_parsed_variants_marked_as_input(tmp_path, ref_fasta, empty_input_dict, for v in empty_input_dict["chr1"].contig_variants[0]: assert v.is_input is True - -# =========================================================================== -# parse_input_vcf — legacy WP genotype (lines 165-176, 186-198) -# =========================================================================== -# NOTE: Lines 165-176 and 186-198 are currently unreachable dead code. -# The guard condition `"WP" in [x.split('=') for x in record[7].split(';')]` -# is always False because it compares the string "WP" against inner lists -# (e.g. ["WP", "0|1"]) — a string never equals a list. -# The correct condition would be: -# any(x.split('=')[0] == "WP" for x in record[7].split(';')) -# The tests below document the CURRENT (broken) behaviour: WP genotypes -# fall through to random genotype generation instead of being parsed. - -def test_wp_in_info_no_format_uses_random_genotype(tmp_path, ref_fasta, empty_input_dict, opts): - """VCF with WP in INFO but no FORMAT column: WP is not recognised (bug), - so a random genotype is generated instead. - - TODO (post-fix): Once the WP guard condition is corrected from - "WP" in [x.split('=') for x in record[7].split(';')] - to - any(x.split('=')[0] == "WP" for x in record[7].split(';')) - update this test to assert that the genotype IS read from the WP field: - np.testing.assert_array_equal(variants[0].genotype, [0, 1]) - and remove the random-genotype assertions below. - """ - vcf = _write_vcf(tmp_path, "wp_noformat.vcf", _vcf_header_no_format() + [ - "chr1\t1\t.\tA\tG\t30\tPASS\tWP=0|1", - ]) - parse_input_vcf(empty_input_dict, vcf, 2, ref_fasta, opts) - variants = empty_input_dict["chr1"].contig_variants.get(0, []) - assert len(variants) == 1 - # Genotype is generated randomly — not read from WP (due to the bug) - assert variants[0].genotype is not None - assert len(variants[0].genotype) == 2 - - -def test_wp_in_info_with_format_no_gt_uses_random_genotype(tmp_path, ref_fasta, empty_input_dict, opts): - """VCF with WP in INFO and FORMAT column but no GT field: WP is not - recognised (bug), so a random genotype is generated instead. - - TODO (post-fix): Once the WP guard condition is corrected (see above), - update this test to assert that the genotype IS read from the WP field: - np.testing.assert_array_equal(variants[0].genotype, [0, 1]) - Also verify the FORMAT column is prefixed with "GT:" and the sample - field includes the WP-derived genotype string. - """ - vcf = _write_vcf(tmp_path, "wp_format.vcf", _vcf_header_with_format() + [ - "chr1\t1\t.\tA\tG\t30\tPASS\tWP=0|1\tDP\t42", - ]) - parse_input_vcf(empty_input_dict, vcf, 2, ref_fasta, opts) - variants = empty_input_dict["chr1"].contig_variants.get(0, []) - assert len(variants) == 1 - assert variants[0].genotype is not None - assert len(variants[0].genotype) == 2 - - -def test_wp_only_info_field_not_mistaken_for_gt(tmp_path, ref_fasta, empty_input_dict, opts): - """Confirm that a standalone WP field in INFO without = is also not parsed. - - TODO (post-fix): A bare "WP" with no value is malformed; after the fix - this test should still produce a random genotype (WP= is required). - """ - vcf = _write_vcf(tmp_path, "wp_bare.vcf", _vcf_header_no_format() + [ - "chr1\t1\t.\tA\tG\t30\tPASS\tWP", - ]) - parse_input_vcf(empty_input_dict, vcf, 2, ref_fasta, opts) - variants = empty_input_dict["chr1"].contig_variants.get(0, []) - assert len(variants) == 1 - assert variants[0].genotype is not None \ No newline at end of file diff --git a/tests/test_variants/test_contig_variants.py b/tests/test_variants/test_contig_variants.py index 9145377a..b6337824 100644 --- a/tests/test_variants/test_contig_variants.py +++ b/tests/test_variants/test_contig_variants.py @@ -11,9 +11,7 @@ from neat.variants import Deletion, Insertion, SingleNucleotideVariant from neat.variants.unknown_variant import UnknownVariant -# --------------------------------------------------------------------------- # Helpers -# --------------------------------------------------------------------------- _SEQ = "ACGTACGTACGTACGT" # 16 bp _REC = SeqRecord(Seq(_SEQ), id="chr1", name="chr1", description="") @@ -32,9 +30,7 @@ def _ins(pos, alt="ACGT", length=3, gt=None): return Insertion(pos, length, alt, gt if gt is not None else _GT.copy(), "42") -# =========================================================================== # get_ref_alt — SNV -# =========================================================================== def test_get_ref_alt_snv_ref_is_single_base(): snv = _snv(2, "T") @@ -55,9 +51,7 @@ def test_get_ref_alt_snv_with_block_start_offset(): assert ref == _SEQ[2] # local index = 6 - 4 = 2 -# =========================================================================== # get_ref_alt — Deletion -# =========================================================================== def test_get_ref_alt_deletion_ref_spans_length(): d = _del(1, 3) @@ -72,9 +66,7 @@ def test_get_ref_alt_deletion_alt_is_single_base(): assert alt == _SEQ[1] -# =========================================================================== # get_ref_alt — Insertion -# =========================================================================== def test_get_ref_alt_insertion_ref_is_single_base(): ins = _ins(3, "ACGTT", 4) @@ -88,9 +80,7 @@ def test_get_ref_alt_insertion_alt_from_variant(): assert alt == "ACGTT" -# =========================================================================== # get_ref_alt — UnknownVariant -# =========================================================================== def test_get_ref_alt_unknown_uses_metadata(): uv = UnknownVariant(5, _GT.copy(), "42", is_input=True, @@ -103,9 +93,7 @@ def test_get_ref_alt_unknown_uses_metadata(): assert alt == "ACGT" -# =========================================================================== # get_sample_info -# =========================================================================== def test_get_sample_info_with_neat_sample_metadata(): snv = _snv(2, "T") @@ -120,9 +108,7 @@ def test_get_sample_info_without_metadata_uses_genotype_string(): assert "|" in result or "/" in result -# =========================================================================== # remove_variant -# =========================================================================== def test_remove_variant_method_exists(): """remove_variant silently no-ops due to variant.position bug. @@ -140,9 +126,7 @@ def test_remove_variant_method_exists(): assert callable(cv.remove_variant) -# =========================================================================== # compile_genotypes_for_location -# =========================================================================== def test_compile_genotypes_two_variants_different_ploids(): cv = ContigVariants() @@ -162,9 +146,7 @@ def test_compile_genotypes_single_variant(): assert list(result) == [0, 1] -# =========================================================================== # generate_field -# =========================================================================== def test_generate_field_uses_metadata_when_present(): cv = ContigVariants() @@ -179,9 +161,7 @@ def test_generate_field_falls_back_to_default(): assert cv.generate_field(snv, "ID") == "." -# =========================================================================== # check_if_del / check_if_ins -# =========================================================================== def test_check_if_del_finds_containing_deletion(): cv = ContigVariants() diff --git a/tests/test_variants/test_ref_eq_alt.py b/tests/test_variants/test_ref_eq_alt.py new file mode 100644 index 00000000..bb34de88 --- /dev/null +++ b/tests/test_variants/test_ref_eq_alt.py @@ -0,0 +1,170 @@ +""" +Regression tests for Issue #266 — REF == ALT in output VCF. + +xfail tests indicate each bug exists before fixes are applied. +Once a fix is in place, remove the corresponding xfail marker. +""" + +import io +import tempfile +from pathlib import Path +from types import SimpleNamespace + +import numpy as np +import pytest +from Bio.Seq import Seq +from Bio.SeqRecord import SeqRecord +from numpy.random import default_rng + +from neat.models.mutation_model import MutationModel +from neat.variants import SingleNucleotideVariant +from neat.variants.contig_variants import ContigVariants + + +# Helpers + +_SEQ = "ACGTACGTACGTACGT" +_REC = SeqRecord(Seq(_SEQ), id="chr1", name="chr1", description="") + + +# Issue #266 — generate_snv can return a ref base as ALT + +def _all_diagonal_model(): + """MutationModel whose trinuc_trans_matrices are identity (100% on diagonal).""" + diagonal_matrix = np.eye(4) + all_diagonal = np.stack([diagonal_matrix] * 16) + return MutationModel(trinuc_trans_matrices=all_diagonal) + + +def test_generate_snv_diagonal_model_cannot_produce_ref_eq_alt_regression(): + """Asserts fixed behavior (alt != central base); xfails because bug produces alt == central.""" + model = _all_diagonal_model() + # ACA has a central base of 'C'. With the identity matrix, rng.choice always picks 'C'. + snv = model.generate_snv(Seq("ACA"), reference_location=5, rng=default_rng(0)) + # Fixed behavior: alt must not equal the ref base. Currently fails. + assert snv.alt != "C", "generate_snv returned ref base as ALT (REF==ALT bug)" + + +def test_generate_snv_default_model_never_produces_ref_eq_alt(): + """Default model: alt must never equal the central (reference) base.""" + model = MutationModel() + rng = default_rng(42) + for trinuc in ["ACA", "GCG", "TAT", "CGC", "AGA", "TGT", "ACG", "GCA"]: + central = trinuc[1] + snv = model.generate_snv(Seq(trinuc), reference_location=10, rng=rng) + assert snv.alt != central, ( + f"generate_snv produced REF==ALT ({central!r}) for trinuc {trinuc!r}" + ) + + +def test_generate_snv_diagonal_model_avoids_ref_base_after_fix(): + """After fix: diagonal custom model must not return the reference base as ALT.""" + model = _all_diagonal_model() + rng = default_rng(0) + for trinuc in ["ACA", "GCG", "TAT", "CGC"]: + central = trinuc[1] + snv = model.generate_snv(Seq(trinuc), reference_location=5, rng=rng) + assert snv.alt != central, ( + f"generate_snv produced REF==ALT ({central!r}) for trinuc {trinuc!r}" + ) + + +def test_generate_snv_near_diagonal_model_avoids_ref_base_after_fix(): + """After fix: even 99%-diagonal custom model must not produce REF==ALT.""" + near_diag = np.full((4, 4), 0.01 / 3) + np.fill_diagonal(near_diag, 0.99) + model = MutationModel(trinuc_trans_matrices=np.stack([near_diag] * 16)) + rng = default_rng(7) + for _ in range(200): + snv = model.generate_snv(Seq("ACA"), reference_location=5, rng=rng) + assert snv.alt != "C", "generate_snv produced REF==ALT with near-diagonal model" + + +# Issue #266 — parse_input_vcf accepts REF==ALT variants from user VCFs + +def _write_vcf(tmp_path: Path, name: str, lines: list) -> Path: + p = tmp_path / name + p.write_text("\n".join(lines) + "\n") + return p + + +def _make_opts(tmp_path): + from neat.read_simulator.utils.options import Options + opts = Options(rng_seed=42) + opts.ploidy = 2 + opts.produce_vcf = True + opts.vcf = tmp_path / "out.vcf.gz" + return opts + + +def test_parse_input_vcf_rejects_ref_eq_alt_regression(tmp_path): + """Asserts fixed behavior (variant skipped) with xfails because bug accepts REF==ALT.""" + from Bio import SeqIO + from neat.read_simulator.utils.vcf_func import parse_input_vcf + + fa = tmp_path / "ref.fa" + fa.write_text(f">chr1\n{_SEQ}\n") + ref_fasta = SeqIO.index(str(fa), "fasta") + + vcf = _write_vcf(tmp_path, "refalt.vcf", [ + "##fileformat=VCFv4.2", + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO", + "chr1\t1\t.\tA\tA\t30\tPASS\t.", + ]) + input_dict = {"chr1": ContigVariants()} + parse_input_vcf(input_dict, vcf, 2, ref_fasta, _make_opts(tmp_path)) + # Fixed behavior: variant should be skipped (locations == 0). Currently fails (bug accepts it). + assert len(input_dict["chr1"].variant_locations) == 0, \ + "REF==ALT variant was accepted — bug still present" + + +def test_parse_input_vcf_skips_ref_eq_alt(tmp_path): + """After fix: parse_input_vcf skips REF==ALT variants with a warning.""" + from Bio import SeqIO + from neat.read_simulator.utils.vcf_func import parse_input_vcf + + fa = tmp_path / "ref.fa" + fa.write_text(f">chr1\n{_SEQ}\n") + ref_fasta = SeqIO.index(str(fa), "fasta") + + vcf = _write_vcf(tmp_path, "refalt.vcf", [ + "##fileformat=VCFv4.2", + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO", + "chr1\t1\t.\tA\tA\t30\tPASS\t.", + ]) + input_dict = {"chr1": ContigVariants()} + parse_input_vcf(input_dict, vcf, 2, ref_fasta, _make_opts(tmp_path)) + assert len(input_dict["chr1"].variant_locations) == 0 + + +def test_parse_input_vcf_accepts_valid_snv(tmp_path): + """After fix: valid SNVs (REF != ALT) are still accepted normally.""" + from Bio import SeqIO + from neat.read_simulator.utils.vcf_func import parse_input_vcf + + fa = tmp_path / "ref.fa" + fa.write_text(f">chr1\n{_SEQ}\n") + ref_fasta = SeqIO.index(str(fa), "fasta") + + vcf = _write_vcf(tmp_path, "valid.vcf", [ + "##fileformat=VCFv4.2", + "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO", + "chr1\t1\t.\tA\tG\t30\tPASS\t.", + ]) + input_dict = {"chr1": ContigVariants()} + parse_input_vcf(input_dict, vcf, 2, ref_fasta, _make_opts(tmp_path)) + assert len(input_dict["chr1"].variant_locations) == 1 + + +# Issue #266 — get_ref_alt / write path returns REF==ALT without any guard + +def test_get_ref_alt_snv_ref_eq_alt_is_possible(): + """Demonstrate that get_ref_alt() returns REF==ALT for a badly constructed SNV. + + This is not xfail — it documents that the data model allows the condition. + The guard must exist in write_block_vcf, not get_ref_alt itself. + """ + # _SEQ[0] == 'A'; get_ref_alt returns ('A', 'A') + snv = SingleNucleotideVariant(0, "A", np.array([0, 1]), 40) + ref, alt = ContigVariants.get_ref_alt(snv, _REC, 0) + assert str(ref) == str(alt) == "A"