From a8b632ded60f8f2bfdceae4a359466388d51e1f6 Mon Sep 17 00:00:00 2001 From: Alex Rubinsteyn Date: Tue, 6 Oct 2026 12:51:53 -0400 Subject: [PATCH] Fix RNA variant support cache keys (#575) --- CHANGELOG.md | 7 +++ tests/test_osteosarc_rna_phasing.py | 5 +++ tests/test_rna_haplotype_sequences.py | 34 +++++++++++++++ tests/test_rna_read_phasing_source.py | 62 +++++++++++++++++++++++++++ varcode/rna_read_phasing.py | 15 ++++--- varcode/version.py | 2 +- 6 files changed, 117 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 575431b..2a309db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # Change Log +## [v10.11.7](https://github.com/openvax/varcode/tree/v10.11.7) (2026-10-06) + +- Cache RNA variant support counts under the variant key, keeping it separate + from sequencing-template keys (#575). Verify count accuracy, reuse of zero + and missing-contig results, haplotype cache invalidation, and caching on real + Sid RNA templates without changing the allele or fragment filtering rules. + ## [v10.11.6](https://github.com/openvax/varcode/tree/v10.11.6) (2026-10-05) - Exercise real Sid RNA evidence through `VariantCollection.effects()` for diff --git a/tests/test_osteosarc_rna_phasing.py b/tests/test_osteosarc_rna_phasing.py index adfcbfd..d8260a9 100644 --- a/tests/test_osteosarc_rna_phasing.py +++ b/tests/test_osteosarc_rna_phasing.py @@ -9,6 +9,7 @@ import hashlib import json from pathlib import Path +from unittest.mock import Mock import pytest from pyensembl import cached_release @@ -136,6 +137,8 @@ def test_real_sid_rna_fragment_evidence( # Q20/MAPQ20, five-base edge exclusion and SAM flag filters remain # at their defaults. Allow a 5% allele error rate for the ONT data. phasing_error_rate=0.01 if member == SHORT_NTF3 else 0.05) + fetch = Mock(wraps=source._fetch_reads_for_variant) + source._fetch_reads_for_variant = fetch try: assert source._phase_counts(*variants) == counts assert tuple(source.supports_variant(v) for v in variants) == support @@ -146,6 +149,8 @@ def test_real_sid_rna_fragment_evidence( assert resolver.in_cis(*variants) is phase for variant, partner in (variants, tuple(reversed(variants))): assert tuple(resolver.phased_partners(variant)) == ((partner,) if phase is True else ()) + assert tuple(source.supports_variant(v) for v in reversed(variants)) == support[::-1] + assert fetch.call_count == len(variants) finally: source.close() diff --git a/tests/test_rna_haplotype_sequences.py b/tests/test_rna_haplotype_sequences.py index be3c665..04f451d 100644 --- a/tests/test_rna_haplotype_sequences.py +++ b/tests/test_rna_haplotype_sequences.py @@ -1,6 +1,7 @@ """Known MAP2 RNA haplotypes are sequence observations, not gap mechanisms.""" from types import SimpleNamespace +from unittest.mock import Mock import pytest @@ -62,6 +63,39 @@ def test_map2_identical_sequence_across_gap_decompositions(tmp_path, variants, c s.close() +@pytest.mark.parametrize("cigar,sequence,before,after,phase", [ + ("9M28N6M", COMBINED, (1, 1, 0), (1, 1, 1), True), + ("9M28D6M", "TAGACAGGGCCCATT", (1, 1, 1), (0, 0, 0), None), +]) +def test_haplotype_registration_invalidates_support_and_phase_caches( + tmp_path, variants, cigar, sequence, before, after, phase): + s = source(tmp_path, [read(cigar, sequence)]) + fetch = Mock(wraps=s._fetch_reads_for_variant) + s._fetch_reads_for_variant = fetch + fetch_pair = Mock(wraps=s._fetch_reads_for_pair) + s._fetch_reads_for_pair = fetch_pair + try: + for _ in range(2): + assert tuple(s.supports_variant(v) for v in variants) == before + assert s.in_cis(variants[0], variants[-1]) is (True if before[-1] else None) + assert fetch.call_count == 3 + assert fetch_pair.call_count == 1 + s.register_haplotype(variants) + for _ in range(2): + assert tuple(s.supports_variant(v) for v in variants) == after + assert s.in_cis(variants[0], variants[-1]) is phase + assert fetch.call_count == 6 + assert fetch_pair.call_count == 2 + # Re-registering the same hypothesis does not change allele calls. + s.register_haplotype(variants) + assert tuple(s.supports_variant(v) for v in variants) == after + assert s.in_cis(variants[0], variants[-1]) is phase + assert fetch.call_count == 6 + assert fetch_pair.call_count == 2 + finally: + s.close() + + @pytest.mark.parametrize("quality", [None, 0, 19]) @pytest.mark.parametrize("cigar", ["9M28D6M", "9M28N6M", "5M22D2M6D8M"]) def test_registered_haplotype_does_not_bypass_quality_filters(tmp_path, variants, quality, cigar): diff --git a/tests/test_rna_read_phasing_source.py b/tests/test_rna_read_phasing_source.py index ea50b80..5f4880c 100644 --- a/tests/test_rna_read_phasing_source.py +++ b/tests/test_rna_read_phasing_source.py @@ -10,6 +10,8 @@ # See the License for the specific language governing permissions and # limitations under the License. +from unittest.mock import Mock + import pytest from varcode import ( @@ -114,6 +116,66 @@ def test_supports_variant_counts_alt_fragments_with_contig_normalization(tmp_pat source.close() +@pytest.mark.parametrize("reverse", [False, True]) +@pytest.mark.parametrize("conflict", ["A", "C"]) +def test_support_cache_keeps_variant_counts_separate(tmp_path, reverse, conflict): + reads = [] + # Two templates share a QNAME across libraries. Each has two overlapping + # supporting mates, so they contribute two fragments, not four reads. + for group in ("library1", "library2"): + for flag in (99, 147): + read = _snv_read("shared", 100, {120: "T"}, flag=flag) + read.set_tag("RG", group) + reads.append(read) + reads.extend([ + _snv_read("other_alt", 100, {120: "G"}), + _snv_read("conflict", 100, {120: "T"}, flag=99), + _snv_read("conflict", 100, {120: conflict}, flag=147), + _snv_read("low_quality_mate", 100, {120: "T"}, flag=99), + _snv_read("low_quality_mate", 100, {120: "A"}, + flag=147, low_quality_pos=120), + _snv_read("low_quality_only", 180, {200: "T"}, low_quality_pos=200), + ]) + for name, flag in (("duplicate", 1024), ("secondary", 256), + ("supplementary", 2048), ("improper_pair", 65)): + reads.append(_snv_read(name, 100, {120: "T"}, flag=flag)) + low_mapping_quality = _snv_read("low_mapping_quality", 100, {120: "T"}) + low_mapping_quality.mapping_quality = 19 + reads.append(low_mapping_quality) + variants = [ + (Variant("1", 120, "A", "T"), 3), + (Variant("1", 120, "A", "G"), 1), + (Variant("1", 160, "A", "G"), 0), # Only reference calls. + (Variant("1", 200, "A", "T"), 0), # Only low-quality calls. + (Variant("1", 500, "A", "G"), 0), # No overlapping records. + (Variant("2", 120, "A", "T"), None), # Absent contig. + ] + if reverse: + variants.reverse() + source = RNAReadPhasingSource( + _write_bam(tmp_path, reads), max_distance_from_read_edge=0) + fetch = Mock(wraps=source._fetch_reads_for_variant) + source._fetch_reads_for_variant = fetch + classify = Mock(wraps=source._read_allele) + source._read_allele = classify + try: + for variant, expected in variants: + assert source.supports_variant(variant) == expected + assert fetch.call_count == len(variants) + classifications = classify.call_count + # Interleave queries in the opposite order, using equivalent new + # Variant instances to ensure caching is by allele, not object identity. + for variant, expected in reversed(variants): + equivalent = Variant(variant.contig, variant.start, variant.ref, + variant.alt, genome=variant.genome) + assert source.supports_variant(equivalent) == expected + assert source.has_evidence(equivalent) is (expected is not None and expected >= 2) + assert fetch.call_count == len(variants) + assert classify.call_count == classifications + finally: + source.close() + + def test_read_phase_resolver_reports_cis_from_same_read_alt_cooccurrence(tmp_path): v1 = Variant("1", 120, "A", "T") v2 = Variant("1", 140, "A", "G") diff --git a/varcode/rna_read_phasing.py b/varcode/rna_read_phasing.py index ba112b1..42b5421 100644 --- a/varcode/rna_read_phasing.py +++ b/varcode/rna_read_phasing.py @@ -487,23 +487,24 @@ def supports_variant(self, variant) -> Optional[int]: """Count quality-filtered reads/fragments supporting ``variant.alt``. Returns ``None`` when the variant's contig is absent from the BAM. + Results are cached until a new haplotype is registered. """ self._register_variant(variant) - key = self._variant_key(variant) - if key in self._support_cache: - return self._support_cache[key] + variant_key = self._variant_key(variant) + if variant_key in self._support_cache: + return self._support_cache[variant_key] reads = self._fetch_reads_for_variant(variant) if reads is None: - self._support_cache[key] = None + self._support_cache[variant_key] = None return None calls_by_fragment = defaultdict(set) for read in reads: call = self._read_allele(read, variant) if call is not None: - key = (read.get_tag("RG") if read.has_tag("RG") else "", read.query_name) - calls_by_fragment[key].add(call) + fragment_key = (read.get_tag("RG") if read.has_tag("RG") else "", read.query_name) + calls_by_fragment[fragment_key].add(call) count = sum(calls == {"alt"} for calls in calls_by_fragment.values()) - self._support_cache[key] = count + self._support_cache[variant_key] = count return count def has_evidence(self, variant) -> bool: diff --git a/varcode/version.py b/varcode/version.py index 900edd1..12c43bc 100644 --- a/varcode/version.py +++ b/varcode/version.py @@ -1 +1 @@ -__version__ = "10.11.6" +__version__ = "10.11.7"