Reject contradictory RNA fragment alleles in phasing and support (10.10.2)

.github/workflows/tests.yml

@@ -6,7 +6,10 @@
 # - get a badge for tests passing
 # - download binary dependencies from conda
 name: Tests
-on: [push, pull_request]
+on:
+  push:
+    branches: [main]
+  pull_request:
 
 jobs:
   build:

requirements.txt

@@ -1,6 +1,6 @@
 numpy>=2.0.0,<3.0.0
-pandas>=2.0.0,<3.0.0
+pandas>=2.0.0,<4.0.0
 pyensembl>=2.6.4
 memoized_property>=1.0.2
 serializable>=1.1.0
-sercol>=0.1.4
+sercol>=1.0.3

tests/test_mouse.py

@@ -18,12 +18,12 @@
 from .data import data_path
 
 MOUSE_ENSEMBL_RELEASE = 95
-SERVER = "ftp://ftp.ensembl.org"
+SERVER = "https://ftp.ensembl.org"
 MOUSE_GTF_PATH = \
-    SERVER + "/pub/release-%d/gtf/mus_musculus/Mus_musculus.GRCm38.%d.gtf.gz" % (
+    SERVER + "/pub/release-%d/gtf/mus_musculus/Mus_musculus.GRCm38.%d.chr_patch_hapl_scaff.gtf.gz" % (
         MOUSE_ENSEMBL_RELEASE, MOUSE_ENSEMBL_RELEASE)
 MOUSE_TRANSCRIPT_FASTA_PATH = \
-    SERVER + "/pub/release-%d/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz"
+    SERVER + "/pub/release-%d/fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz" % MOUSE_ENSEMBL_RELEASE
 MOUSE_PROTEIN_FASTA_PATH = \
     SERVER + "/pub/release-%d/fasta/mus_musculus/pep/Mus_musculus.GRCm38.pep.all.fa.gz" % (
         MOUSE_ENSEMBL_RELEASE)

tests/test_rna_haplotype_sequences.py

@@ -143,3 +143,30 @@ def test_no_transitive_join_of_registered_haplotypes(tmp_path):
         assert s.in_cis(a, b) is None
     finally:
         s.close()
+
+
+@pytest.mark.parametrize("quality,expected", [(40, None), (5, True)])
+def test_anchored_match_cannot_override_an_informative_contradictory_mate(tmp_path, variants, quality, expected):
+    matching = read("9M28N6M")
+    contradictory = read("%dM" % len(REFERENCE), REFERENCE, quality=quality)
+    matching.flag, contradictory.flag = 99, 147
+    s = source(tmp_path, [matching, contradictory])
+    try:
+        s.register_haplotype(variants)
+        assert s.in_cis(variants[0], variants[-1]) is expected
+        assert s.supports_variant(variants[-1]) == (1 if expected else 0)
+    finally:
+        s.close()
+
+
+def test_independent_anchored_fragment_rescues_a_conflicted_pair(tmp_path, variants):
+    matching = read("9M28N6M")
+    contradictory = read("9M28N6M", "AAGACCTGGCCCATT")
+    matching.flag, contradictory.flag = 99, 147
+    s = source(tmp_path, [matching, contradictory, read("9M28N6M", name="clean")])
+    try:
+        s.register_haplotype(variants)
+        assert s.in_cis(variants[0], variants[-1]) is True
+        assert s.supports_variant(variants[-1]) == 1
+    finally:
+        s.close()

tests/test_rna_read_phasing_source.py

@@ -541,3 +541,37 @@ def test_splice_junction_skip_does_not_support_intronic_variant(tmp_path):
         assert source.has_evidence(variant) is False
     finally:
         source.close()
+
+
+@pytest.mark.parametrize("locus", [120, 140])
+@pytest.mark.parametrize("conflict", ["A", "C"])
+@pytest.mark.parametrize("clean_count", [0, 1, 2])
+def test_conflicting_mates_are_not_definitive_fragments(tmp_path, locus, conflict, clean_count):
+    variants = [Variant("1", 120, "A", "T"), Variant("1", 140, "A", "G")]
+    reads = []
+    for i in range(2):
+        name = "conflict%d" % i
+        reads.append(_snv_read(name, 100, {120: "T", 140: "G"}, flag=99))
+        calls = {120: "T", 140: "G", locus: conflict}
+        reads.append(_snv_read(name, 100, calls, flag=147))
+    reads.extend(_snv_read("clean%d" % i, 100, {120: "T", 140: "G"}) for i in range(clean_count))
+    source = RNAReadPhasingSource(_write_bam(tmp_path, reads), max_distance_from_read_edge=0)
+    try:
+        assert source.supports_variant(variants[locus == 140]) == clean_count
+        assert source.in_cis(*variants) is (True if clean_count >= 2 else None)
+        assert source.in_cis(*reversed(variants)) is (True if clean_count >= 2 else None)
+    finally:
+        source.close()
+
+
+def test_low_quality_disagreement_does_not_override_clean_mate(tmp_path):
+    v1, v2 = Variant("1", 120, "A", "T"), Variant("1", 140, "A", "G")
+    reads = [r for i in range(2) for r in (
+        _snv_read("q%d" % i, 100, {120: "T", 140: "G"}, flag=99),
+        _snv_read("q%d" % i, 100, {120: "T", 140: "A"}, low_quality_pos=140, flag=147))]
+    source = RNAReadPhasingSource(_write_bam(tmp_path, reads), max_distance_from_read_edge=0)
+    try:
+        assert source.in_cis(v1, v2) is True
+        assert source.supports_variant(v2) == 2
+    finally:
+        source.close()

varcode/rna_read_phasing.py

@@ -228,6 +228,8 @@ def register_haplotype(self, variants, *, flanking_bases=5):
         alignment, with the configured base-quality and read-edge filters.
         Equivalent D/N/split-gap encodings can then support the same known
         sequence. This does not establish a genomic deletion from an RNA skip.
+        A fully covered high-quality mate with a different sequence vetoes that
+        fragment's support. Missing coverage and low-quality context do not.
         Reference is fetched from the variants' genome; unavailable sequence,
         overlapping edits and reference mismatches raise ValueError. No BAM
         bases are corrected, no missing bases filled, and no mates assembled.
@@ -268,20 +270,21 @@ def register_haplotype(self, variants, *, flanking_bases=5):
             self._support_cache.clear()
             self._phase_cache.clear()
 
-    def _matches_haplotype(self, read, left0, right0, sequence):
+    def _haplotype_allele(self, read, left0, right0, sequence):
         anchors = {}
         for operation, start, end, query, _ in self._cigar_events(read):
             if operation in (_CIGAR_MATCH, _CIGAR_EQUAL, _CIGAR_DIFF):
                 for position in (left0, right0):
                     if start <= position < end:
                         anchors[position] = query + position - start
         if left0 not in anchors or right0 not in anchors:
-            return False
+            return None
         left, right = anchors[left0], anchors[right0]
-        if right < left or right - left + 1 != len(sequence):
-            return False
-        return (self._base_calls_ok(read, range(left, right + 1))
-                and read.query_sequence[left:right + 1].upper() == sequence)
+        if right < left or not self._base_calls_ok(read, range(left, right + 1)):
+            return None
+        # A fully observed, high-quality nonmatch can contradict another mate,
+        # but must never be interpreted as reference/trans evidence.
+        return "alt" if read.query_sequence[left:right + 1].upper() == sequence else "other"
 
     def _register_variant(self, variant):
         key = self._variant_key(variant)
@@ -389,7 +392,7 @@ def _substitution_allele(self, read, variant):
             return "alt"
         if observed.upper() == variant.ref.upper():
             return "ref"
-        return None
+        return "other"
 
     def _has_exact_deletion(self, read, start0, end0):
         for operation, ref_start, ref_end, _, _ in self._cigar_events(read):
@@ -421,7 +424,7 @@ def _deletion_allele(self, read, variant):
         observed = "".join(read.query_sequence[pos] for pos in query_positions)
         if observed.upper() == variant.ref.upper():
             return "ref"
-        return None
+        return "other"
 
     def _inserted_positions_after_anchor(self, read, anchor0):
         inserted = []
@@ -449,7 +452,7 @@ def _insertion_allele(self, read, variant):
             inserted = "".join(read.query_sequence[pos] for pos in inserted_positions)
             if inserted.upper() == variant.alt.upper():
                 return "alt"
-            return None
+            return "other"
         if self._base_calls_ok(read, [anchor_query]):
             return "ref"
         return None
@@ -458,14 +461,12 @@ def _read_allele(self, read, variant):
         if not self._passes_read_filters(read) or read.query_sequence is None:
             return None
         key = self._variant_key(variant)
-        registered = False
-        for keys, left, right, sequence in self._haplotypes:
-            if key in keys:
-                registered = True
-                if self._matches_haplotype(read, left, right, sequence):
-                    return "alt"
-        if registered:
-            return None
+        haplotype_calls = [self._haplotype_allele(read, left, right, sequence)
+                           for keys, left, right, sequence in self._haplotypes if key in keys]
+        if haplotype_calls:
+            if "alt" in haplotype_calls:
+                return "alt"
+            return "other" if "other" in haplotype_calls else None
         if variant.is_insertion:
             return self._insertion_allele(read, variant)
         if variant.is_deletion:
@@ -495,12 +496,13 @@ def supports_variant(self, variant) -> Optional[int]:
         if reads is None:
             self._support_cache[key] = None
             return None
-        supporting_fragments = set()
+        calls_by_fragment = defaultdict(set)
         for read in reads:
-            if self._read_allele(read, variant) == "alt":
-                supporting_fragments.add((
-                    read.get_tag("RG") if read.has_tag("RG") else "", read.query_name))
-        count = len(supporting_fragments)
+            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)
+        count = sum(calls == {"alt"} for calls in calls_by_fragment.values())
         self._support_cache[key] = count
         return count
 
@@ -532,27 +534,19 @@ def _fragment_alleles_for_pair(self, v1, v2):
         anchored = any(self._variant_key(v) in keys
                        for v in (v1, v2) for keys, _, _, _ in self._haplotypes)
         for fragment_reads in grouped.values():
-            if anchored:
-                # Never assemble a new anchored haplotype from mates or from
-                # different registered combinations through a shared allele.
-                pairs = {tuple(self._read_allele(read, v) for v in (v1, v2))
-                         for read in fragment_reads}
-                fragments.append(("alt", "alt") if ("alt", "alt") in pairs
-                                 else (None, None))
-                continue
+            pairs = [tuple(self._read_allele(read, v) for v in (v1, v2))
+                     for read in fragment_reads]
             alleles = []
-            for variant in (v1, v2):
-                calls = [
-                    self._read_allele(read, variant)
-                    for read in fragment_reads
-                ]
-                if "alt" in calls:
-                    alleles.append("alt")
-                elif "ref" in calls:
-                    alleles.append("ref")
-                else:
-                    alleles.append(None)
-            fragments.append(tuple(alleles))
+            for calls in zip(*pairs):
+                informative = set(calls) - {None}
+                alleles.append(next(iter(informative)) if len(informative) == 1 else None)
+            if anchored:
+                # Require one full-context alignment, and no informative mate
+                # contradicting either allele. Never assemble partial contexts.
+                unambiguous_match = tuple(alleles) == ("alt", "alt") and ("alt", "alt") in pairs
+                fragments.append(("alt", "alt") if unambiguous_match else (None, None))
+            else:
+                fragments.append(tuple(alleles))
         return fragments
 
     def _phase_counts(self, v1, v2):

varcode/version.py

@@ -1 +1 @@
-__version__ = "10.10.0"
+__version__ = "10.10.2"