Merge pull request #267 from ncsa/vcf-fixes

Solving VCF-related issues.

Author: joshfactorial

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)

neat/read_simulator/utils/stitch_outputs.py

@@ -29,11 +29,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"

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:

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

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}"

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}