diff --git a/CHANGELOG.md b/CHANGELOG.md index 070bb52..68cfe1d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # Change Log +## [v10.11.1](https://github.com/openvax/varcode/tree/v10.11.1) (2026-09-30) + +- Compose mixed inherited/somatic cis groups in the experimental transcript + model (#500). Inherited alleles anchor the patient baseline and are edited + once; only novel alleles contribute new somatic edits. Repeated normalized + alleles are deduplicated without losing the supplied group membership. +- Resolve germline phase through any member of a known-cis group, preserving + uncertainty for unlinked alleles and phase provenance on each hypothesis. + ## [v10.11.0](https://github.com/openvax/varcode/tree/v10.11.0) (2026-09-29) - `Genome(native_genome)` inherits optional native PyEnsembl reference DNA diff --git a/docs/experimental_annotators.md b/docs/experimental_annotators.md index 7738e3b..d4c4f8d 100644 --- a/docs/experimental_annotators.md +++ b/docs/experimental_annotators.md @@ -48,6 +48,27 @@ the phase cap too when `predict_transcript_model_effect` is called directly. Exceeding either returns a [`HypothesisLimit`](germline.md#phase-enumeration-limit) rather than a partial set of candidates. +### Mixed inherited and somatic haplotypes + +`predict_transcript_model_effect(variants, transcript, germline_variants=...)` +and `annotate_haplotype` take a known-cis group. Members matching the supplied +germline are included once in the patient baseline. Only novel alleles are +applied to that baseline to obtain the mutant; `effect.variant` is the first +novel member and `effect.variants` retains the supplied group. Matching uses +reference assembly, contig, normalized position, REF and ALT. A group with no +novel allele returns `GermlineAlleleOverlap`, without inferring LOH. + +The group itself supplies cis constraints. Phase can propagate from any member +to other germline alleles; unlinked alleles still produce separate hypotheses. +Homozygous alleles are on both haplotypes and do not bridge their phase. +As with other phase constraints, conflicting resolver answers are logged and +ignored. Each hypothesis retains its germline cis/trans assignments, phase +state, resolver source and novel `somatic_variants` in its evidence. + +`VariantCollection.effects` constructs these groups only when the resolver +supports cis. Unknown or trans relationships retain individual predictions; +listing alleles in a collection alone does not assert that they are in cis. + ## Transcript-model results Use the [sequence accessors](transcript_models.md#sequence-access) for a single diff --git a/tests/test_germline_effects.py b/tests/test_germline_effects.py index e00cfd0..e5b4a04 100644 --- a/tests/test_germline_effects.py +++ b/tests/test_germline_effects.py @@ -456,6 +456,33 @@ def _germline_vcf(genotype): class TestPartitionGermlineByPhase: + def test_cis_group_anchors_germline_through_any_member(self): + somatic = _somatic_at() + partner = _germline_at(102) + germ = (_germline_at(101), _germline_at(103)) + phase = partition_germline_by_phase( + somatic, germ, _PartialPhaseResolver({}, {(101, 102): False}), + cis_variants=(somatic, partner)) + assert phase.trans == germ[:1] + assert phase.unphased == germ[1:] + + def test_homozygous_member_does_not_anchor_other_germline(self): + somatic = _somatic_at() + germ = (_germline_at(101), _germline_at(102)) + phase = partition_germline_by_phase( + somatic, germ, _PartialPhaseResolver({}, {(101, 102): True}), + homozygous=germ[:1], cis_variants=(somatic, germ[0])) + assert phase.cis == germ[:1] + assert phase.unphased == germ[1:] + + def test_cis_group_constraints_precede_conflicting_resolver_answers(self, caplog): + germ = (_germline_at(101),) + phase = partition_germline_by_phase( + _somatic_at(), germ, _PartialPhaseResolver({101: False}), + cis_variants=germ) + assert phase.cis == germ + assert "contradicts" in caplog.text + def test_without_resolver_everything_is_unphased(self): germ = [_germline_at(101), _germline_at(102)] phase = partition_germline_by_phase(_somatic_at(), germ) diff --git a/tests/test_germline_overlap.py b/tests/test_germline_overlap.py index 7ba7bcc..530679f 100644 --- a/tests/test_germline_overlap.py +++ b/tests/test_germline_overlap.py @@ -7,7 +7,7 @@ Completeness, EffectCollection, GermlineAlleleOverlap, GermlineContext, Variant, detect_germline_overlap, ) -from varcode.effects import Unresolved +from varcode.effects import Substitution from varcode.transcript_model import predict_transcript_model_effect @@ -62,6 +62,9 @@ def test_mixed_group_does_not_apply_inherited_allele_twice(allele): novel = Variant("7", 117531100, "T", "A", genome=81) effect = predict_transcript_model_effect( (novel, allele), transcript, germline_variants=(allele,)) - assert isinstance(effect, Unresolved) - assert effect.mechanism == "germline_overlap_haplotype" - assert effect.is_germline_overlap and effect.is_loh is None + assert isinstance(effect, Substitution) + assert effect.short_description == "p.S159T" + outcome, = effect.candidates[0].outcomes + assert outcome.hypothesis.phase == ((allele, "cis"),) + assert outcome.baseline.protein_sequence[158] == "S" + assert outcome.mutant.protein_sequence[158] == "T" diff --git a/tests/test_mixed_haplotypes.py b/tests/test_mixed_haplotypes.py new file mode 100644 index 0000000..6b846ee --- /dev/null +++ b/tests/test_mixed_haplotypes.py @@ -0,0 +1,208 @@ +"""Mixed inherited/somatic groups use one patient baseline (#500).""" + +import pytest +from pyensembl import cached_release + +from varcode import ( + EffectCollection, GermlineAlleleOverlap, GermlineContext, HypothesisLimit, + TranscriptModelEffectAnnotator, Variant, VariantCollection, + predict_transcript_model_effect, +) + + +class Phase: + phase_source = "mixed_haplotype_test" + + def __init__(self, *links): + self.links = {frozenset((left, right)): answer + for left, right, answer in links} + + def in_cis(self, left, right, transcript=None): + return self.links.get(frozenset((left, right))) + + +@pytest.fixture +def transcript(): + return cached_release(81).transcript_by_id("ENST00000003084") + + +@pytest.fixture(params=["C", "TGCC", ""], ids=["snv", "insertion", "deletion"]) +def alleles(request, transcript): + novel = Variant("7", 117531100, "T", "A", genome=81) + inherited = Variant("7", 117531101, "T", request.param, genome=81) + offset = transcript.spliced_offset(117531101) + baseline = (transcript.sequence[:offset] + request.param + + transcript.sequence[offset + 1:]) + # The novel edit precedes the inherited edit, so this offset is unchanged. + offset = transcript.spliced_offset(novel.start) + mutant = baseline[:offset] + "A" + baseline[offset + 1:] + return novel, inherited, baseline, mutant + + +def outcomes(effect): + return [outcome for candidate in effect.candidates + for outcome in candidate.outcomes] + + +@pytest.mark.parametrize("inherited_first", [False, True]) +@pytest.mark.parametrize("resolved", [False, True]) +def test_mixed_edits_are_applied_once(alleles, transcript, inherited_first, resolved): + novel, inherited, baseline, mutant = alleles + members = ((inherited, novel, inherited) if inherited_first + else (novel, inherited, inherited)) + resolver = Phase((novel, inherited, True)) if resolved else None + effect = predict_transcript_model_effect( + members, transcript, germline_variants=(inherited, inherited), + phase_resolver=resolver, max_phase_hypotheses=1) + + assert effect.variant == novel + assert effect.variants == members + outcome, = outcomes(effect) + assert outcome.baseline.cdna_sequence == baseline + assert outcome.mutant.cdna_sequence == mutant + assert outcome.hypothesis.phase == ((inherited, "cis"),) + assert outcome.hypothesis.evidence["somatic_variants"] == (novel,) + assert outcome.hypothesis.evidence["phase_state"] == "phased" + assert outcome.hypothesis.evidence["phase_source"] == ( + resolver.phase_source if resolved else None) + + +@pytest.mark.parametrize("relation", [True, False, None]) +def test_other_germline_phase_is_preserved(alleles, transcript, relation): + novel, inherited, baseline, mutant = alleles + position = 117531107 + offset = transcript.spliced_offset(position) + ref = transcript.sequence[offset] + other = Variant("7", position, ref, "A" if ref != "A" else "C", genome=81) + # Evidence through the inherited member must propagate to the novel edit. + resolver = Phase((inherited, other, relation)) + effect = predict_transcript_model_effect( + (novel, inherited), transcript, germline_variants=(inherited, other), + phase_resolver=resolver) + + realized = outcomes(effect) + assert len(realized) == (2 if relation is None else 1) + expected_sides = {"cis", "trans"} if relation is None else { + "cis" if relation else "trans"} + assert {dict(o.hypothesis.phase)[other] for o in realized} == expected_sides + # Locate the other site after the inherited indel has shifted it. + shifted = offset + len(baseline) - len(transcript.sequence) + for outcome in realized: + phase = dict(outcome.hypothesis.phase) + assert phase[inherited] == "cis" + base, mut = baseline, mutant + if phase[other] == "cis": + base = base[:shifted] + other.alt + base[shifted + 1:] + mut = mut[:shifted] + other.alt + mut[shifted + 1:] + assert outcome.baseline.cdna_sequence == base + assert outcome.mutant.cdna_sequence == mut + assert outcome.hypothesis.phase_probability == (0.5 if relation is None else 1) + assert outcome.hypothesis.evidence["phase_state"] == ( + "unknown" if relation is None else "phased") + + +def test_haplotype_hook_keeps_repeated_inherited_alleles(alleles, transcript): + novel, inherited, baseline, mutant = alleles + members = (inherited, inherited, novel) + effect = TranscriptModelEffectAnnotator().annotate_haplotype( + members, transcript, germline_ctx=GermlineContext.from_variants([inherited])) + assert effect.variant == novel + assert effect.variants == members + outcome, = outcomes(effect) + assert outcome.baseline.cdna_sequence == baseline + assert outcome.mutant.cdna_sequence == mutant + + +@pytest.mark.parametrize("alt", ["TGAA", ""], ids=["insertion", "deletion"]) +def test_novel_indels_use_the_same_patient_baseline(alleles, transcript, alt): + _, inherited, baseline, _ = alleles + novel = Variant("7", 117531100, "T", alt, genome=81) + effect = predict_transcript_model_effect( + (inherited, novel, inherited), transcript, germline_variants=(inherited,)) + outcome, = outcomes(effect) + offset = transcript.spliced_offset(117531100) + assert outcome.baseline.cdna_sequence == baseline + assert outcome.mutant.cdna_sequence == baseline[:offset] + alt + baseline[offset + 1:] + + +@pytest.mark.parametrize("relation", [True, False, None]) +def test_collection_only_composes_supported_cis_group(alleles, transcript, relation): + novel, inherited, baseline, mutant = alleles + resolver = Phase((novel, inherited, relation)) + effects = VariantCollection([inherited, novel, inherited]).effects( + annotator="transcript_model", phase_resolver=resolver, + germline=GermlineContext.from_variants([inherited, inherited])) + local = [e for e in effects if e.transcript == transcript] + joints = [e for e in local if len(getattr(e, "variants", ())) > 1] + assert any(isinstance(e, GermlineAlleleOverlap) for e in local) + if relation is not True: + assert not joints + individual, = [e for e in local if e.variant == novel] + assert {dict(o.hypothesis.phase)[inherited] for o in outcomes(individual)} == ( + {"trans"} if relation is False else {"cis", "trans"}) + return + joint, = joints + assert joint.variant == novel + assert joint.phase_source == resolver.phase_source + outcome, = outcomes(joint) + assert outcome.baseline.cdna_sequence == baseline + assert outcome.mutant.cdna_sequence == mutant + restored = EffectCollection.from_json(EffectCollection([joint]).to_json())[0] + assert restored.variants == joint.variants + assert restored.variant == novel + assert restored.phase_source == resolver.phase_source + + +def test_all_inherited_members_establish_no_new_allele(alleles, transcript): + _, inherited, _, _ = alleles + other = Variant("7", 117531100, "T", "A", genome=81) + members = (inherited, other, inherited) + effect = predict_transcript_model_effect( + members, transcript, germline_variants=(inherited, other, inherited)) + assert isinstance(effect, GermlineAlleleOverlap) + assert effect.variants == members + assert effect.is_loh is None + assert effect.modifies_coding_sequence is False + assert effect.modifies_protein_sequence is False + + +def test_normalized_equivalent_alleles_share_baseline(transcript): + novel = Variant("7", 117531100, "T", "A", genome=81) + inherited = Variant("7", 117531101, "T", "TGCC", genome=81) + unpadded = Variant("7", 117531101, "", "GCC", genome=81) + effect = predict_transcript_model_effect( + (novel, unpadded, inherited), transcript, + germline_variants=(inherited, unpadded)) + outcome, = outcomes(effect) + assert outcome.hypothesis.phase == ((inherited, "cis"),) + assert outcome.hypothesis.evidence["somatic_variants"] == (novel,) + assert len(outcome.baseline.cdna_sequence) == len(transcript.sequence) + 3 + assert len(outcome.mutant.cdna_sequence) == len(transcript.sequence) + 3 + + +def test_mixed_haplotype_retains_phase_limit(transcript): + novel = Variant("7", 117531100, "T", "A", genome=81) + inherited = Variant("7", 117531101, "T", "C", genome=81) + other = Variant("7", 117531102, "G", "A", genome=81) + effect = predict_transcript_model_effect( + (inherited, novel, inherited), transcript, + germline_variants=(inherited, other), max_phase_hypotheses=1) + assert isinstance(effect, HypothesisLimit) + assert effect.variant == novel + assert effect.phase.cis == (inherited,) + assert effect.phase.unphased == (other,) + assert effect.reference_effect is not None + + +def test_reverse_strand_mixed_haplotype(): + transcript = cached_release(81).transcript_by_id("ENST00000357654") + inherited = Variant("17", 43082570, "C", "A", genome=81) + novel = Variant("17", 43082563, "T", "A", genome=81) + effect = predict_transcript_model_effect( + (inherited, novel, inherited), transcript, germline_variants=(inherited,)) + outcome, = outcomes(effect) + baseline = list(transcript.sequence) + baseline[transcript.spliced_offset(inherited.start)] = "T" + assert outcome.baseline.cdna_sequence == "".join(baseline) + baseline[transcript.spliced_offset(novel.start)] = "T" + assert outcome.mutant.cdna_sequence == "".join(baseline) diff --git a/varcode/germline.py b/varcode/germline.py index 10fda66..cf15e2c 100644 --- a/varcode/germline.py +++ b/varcode/germline.py @@ -888,7 +888,7 @@ def _hypothesis(self, cis, phase_state): def partition_germline_by_phase( somatic_variant, germline_variants, phase_resolver=None, - homozygous=()) -> PhasePartition: + homozygous=(), cis_variants=()) -> PhasePartition: """Split ``germline_variants`` by their phase relative to ``somatic_variant``. Variants in ``homozygous`` sit on both haplotypes and are always cis, @@ -899,16 +899,30 @@ def partition_germline_by_phase( germline variant. Variants whose phase is known only relative to each other become ``linked`` groups. An answer that contradicts earlier ones is ignored with a warning. + + ``cis_variants`` supplies additional variants known to be in cis with + the somatic variant. They anchor the phase graph before resolver answers, + so germline phase can follow through any member of that group. + Homozygous alleles do not link the two haplotypes through this graph. """ germline = tuple(germline_variants) + homozygous = tuple(homozygous) if str(somatic_variant.contig).lstrip("chr").upper() in ("M", "MT", "Y"): return PhasePartition(germline=germline, cis=germline, implicit=True) rest = [g for g in germline if g not in homozygous] - # Union-find over the somatic variant (node 0) and ``rest``; parity is - # 0 for the same haplotype as the parent node and 1 for the other one. + # Union-find over the somatic variant (node 0), ``rest``, and any + # additional cis-group members. Parity is 0 for the same haplotype + # as the parent node and 1 for the other one. nodes = [somatic_variant] + rest + cis_variants = tuple(cis_variants) + for variant in cis_variants: + if not detect_germline_overlap(variant, nodes + list(homozygous)): + nodes.append(variant) parent = list(range(len(nodes))) parity = [0] * len(nodes) + for i, variant in enumerate(nodes): + if detect_germline_overlap(variant, cis_variants): + parent[i] = 0 def find(i): if parent[i] != i: diff --git a/varcode/transcript_model.py b/varcode/transcript_model.py index 8ca2f79..9ec623f 100644 --- a/varcode/transcript_model.py +++ b/varcode/transcript_model.py @@ -172,6 +172,15 @@ def _attach_candidate_set(ordered): return top_effect +def _unique_alleles(variants): + """Keep the first occurrence of each normalized allele, in input order.""" + unique = [] + for variant in variants: + if not detect_germline_overlap(variant, unique): + unique.append(variant) + return tuple(unique) + + def predict_transcript_model_effect( variants, transcript, germline_variants=(), phase_resolver=None, sequence_provider=None, max_hypotheses=64, max_phase_hypotheses=None, @@ -189,6 +198,11 @@ def predict_transcript_model_effect( phase/splice outcomes, gives a :class:`~varcode.HypothesisLimit` instead of a partial candidate set. ``homozygous_germline`` lists germline variants on both haplotypes, which are always cis. + + ``variants`` is a known-cis group. Members also present in the germline + anchor the patient baseline; only novel alleles are applied to that + baseline to obtain the mutant. Repeated normalized alleles are edited + once. Other germline alleles retain their phase uncertainty. """ max_hypotheses = _check_max_hypotheses(max_hypotheses) max_phase_hypotheses = max_hypotheses if max_phase_hypotheses is None else ( @@ -196,18 +210,15 @@ def predict_transcript_model_effect( variants = tuple(variants) if not variants: raise ValueError("predict_transcript_model_effect requires a somatic variant") - germline_variants = tuple(germline_variants) - primary = variants[0] - if any(detect_germline_overlap(v, germline_variants) for v in variants): - if len(variants) == 1: - return GermlineAlleleOverlap(primary, transcript) - result = Unresolved( - primary, transcript, mechanism="germline_overlap_haplotype", - reason="A mixed inherited/somatic group requires an allele-aware baseline") - result.is_germline_overlap = True - result.is_loh = None - result.loh_status = "not_assessed" + germline_variants = _unique_alleles(germline_variants) + somatic_variants = tuple( + v for v in _unique_alleles(variants) + if not detect_germline_overlap(v, germline_variants)) + if not somatic_variants: + result = GermlineAlleleOverlap(variants[0], transcript) + result.variants = variants return result + primary = somatic_variants[0] if (getattr(primary, "sv_type", None) in ("DUP", "INV") and primary.alt_assembly): if len(variants) != 1 or germline_variants: @@ -240,11 +251,13 @@ def predict_transcript_model_effect( phase = partition_germline_by_phase( primary, germline_variants, phase_resolver, - homozygous=homozygous_germline) + homozygous=homozygous_germline, cis_variants=variants) try: outcomes = _realize_outcomes( - variants, transcript, phase.hypotheses(max_phase_hypotheses), - sequence_provider, max_hypotheses) + variants, somatic_variants, transcript, + phase.hypotheses(max_phase_hypotheses), + sequence_provider, max_hypotheses, + phase_source=getattr(phase_resolver, "phase_source", None)) except NonlocalStructuralEdit as error: return Unresolved( primary, transcript, mechanism="nonlocal_structural_variant", @@ -262,21 +275,21 @@ def predict_transcript_model_effect( def _realize_outcomes( - variants, transcript, phase_hypotheses, sequence_provider, - max_hypotheses): + variants, somatic_variants, transcript, phase_hypotheses, + sequence_provider, max_hypotheses, phase_source=None): """Classify every phase hypothesis × splice plan. Raises :class:`~varcode.HypothesisLimitError` past ``max_hypotheses`` outcomes and ``NonlocalStructuralEdit`` for edits the layout cannot hold. """ - primary = variants[0] + primary = somatic_variants[0] outcomes = [] enumeration_index = 0 for phase_hypothesis in phase_hypotheses: reference = GenomicLayout.from_transcript( transcript, flank=50, sequence_provider=sequence_provider) baseline_layout = reference.apply_variants(phase_hypothesis.cis) - mutant_layout = baseline_layout.apply_variants(variants) + mutant_layout = baseline_layout.apply_variants(somatic_variants) baseline_runs, baseline_statuses = _status_map( transcript, baseline_layout) mutant_runs, mutant_statuses = _status_map(transcript, mutant_layout) @@ -313,6 +326,8 @@ def _realize_outcomes( phase_hypotheses, phase_hypothesis), evidence={ "phase_state": phase_hypothesis.phase_state, + "phase_source": phase_source, + "somatic_variants": somatic_variants, "shared_splice_sites": tuple(sorted(shared_keys)), }, enumeration_index=enumeration_index) diff --git a/varcode/version.py b/varcode/version.py index 120f813..217d4a5 100644 --- a/varcode/version.py +++ b/varcode/version.py @@ -1 +1 @@ -__version__ = "10.11.0" +__version__ = "10.11.1"