diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index ec6d3c6..f9d497b 100644 --- a/.github/workflows/tests.yml +++ b/.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: diff --git a/requirements.txt b/requirements.txt index c0eef08..09ec848 100644 --- a/requirements.txt +++ b/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 diff --git a/tests/test_mouse.py b/tests/test_mouse.py index 319e916..97e585a 100644 --- a/tests/test_mouse.py +++ b/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) diff --git a/tests/test_rna_haplotype_sequences.py b/tests/test_rna_haplotype_sequences.py index 9f916a6..be3c665 100644 --- a/tests/test_rna_haplotype_sequences.py +++ b/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() diff --git a/tests/test_rna_read_phasing_source.py b/tests/test_rna_read_phasing_source.py index ac8d5a5..ea50b80 100644 --- a/tests/test_rna_read_phasing_source.py +++ b/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() diff --git a/varcode/rna_read_phasing.py b/varcode/rna_read_phasing.py index ef1cf4b..ba112b1 100644 --- a/varcode/rna_read_phasing.py +++ b/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,7 +270,7 @@ 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): @@ -276,12 +278,13 @@ def _matches_haplotype(self, read, left0, right0, sequence): 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): diff --git a/varcode/version.py b/varcode/version.py index e6f9023..1a94747 100644 --- a/varcode/version.py +++ b/varcode/version.py @@ -1 +1 @@ -__version__ = "10.10.0" +__version__ = "10.10.2"