diff --git a/CHANGELOG.md b/CHANGELOG.md index 96b46f2..9778bb5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,17 @@ # Change Log +## [v10.10.0](https://github.com/openvax/varcode/tree/v10.10.0) (2026-09-29) + +- Add `reference_completion_hypotheses` for mapped Exacto and other RNA + observations (#499). Explicitly requested hypotheses prepend/append compatible + reference sequence, retaining the original observation and producer peptide. + Alternative reference isoforms remain separate; evidence distinguishes + observed bases, assumed spans and the original read support. +- Completion respects known RNA ends and predicts from the retained annotated + start using existing genetic-code/selenocysteine rules. Missing starts, + ambiguous coding sequence and missing stops remain explicit. Completion does + not establish RNA coverage, initiation, translation or expression. + ## [v10.9.0](https://github.com/openvax/varcode/tree/v10.9.0) (2026-09-28) - Map Exacto's verified sense exonic runs to reference-transcript segments, diff --git a/docs/api_rna.md b/docs/api_rna.md index 75b8b45..e0a92f1 100644 --- a/docs/api_rna.md +++ b/docs/api_rna.md @@ -19,6 +19,10 @@ See [Transcript models](transcript_models.md) for sequence access and ::: varcode.RNAEvidence +## Reference-completion hypotheses + +::: varcode.reference_completion_hypotheses + ## RNA evidence ::: varcode.RNAEvidenceResolver diff --git a/docs/rna_structures.md b/docs/rna_structures.md index bf3b62d..7f64ee7 100644 --- a/docs/rna_structures.md +++ b/docs/rna_structures.md @@ -169,9 +169,87 @@ fragments and observations without a usable reference frame stay unresolved sequence change; initiation and translation still remain predictions. Mapping does not choose a new ORF or establish initiation at an internal or -upstream ATG. Reference-completed candidates remain future work in -[#499](https://github.com/openvax/varcode/issues/499); initiation-dependent -extensions remain under [#466](https://github.com/openvax/varcode/issues/466). +upstream ATG. Initiation-dependent extensions remain under +[#466](https://github.com/openvax/varcode/issues/466). + +## Make explicit reference-completion hypotheses + +An observed fragment may omit the annotated start or end of a transcript. +Use `reference_completion_hypotheses` to ask what protein a compatible +reference continuation would predict: + +```python +from varcode import reference_completion_hypotheses + +# `rna` is the result of load_exacto_fusions above. Keep these collections +# separate: inferred continuations are not RNA-observed outcomes. +observations = rna.candidates +hypotheses = tuple( + hypothesis + for observation in observations + for hypothesis in reference_completion_hypotheses(observation) +) +for hypothesis in hypotheses: + model = hypothesis.effect.mutant_transcript + spans = hypothesis.evidence["reference_completion"] + start, end = spans["observed_span"] + observed_sequence = model.cdna_sequence[start:end] + print(model.mutant_protein_sequence, spans["assumed_spans"]) +``` + +The helper adds a reference prefix and/or suffix only when that terminal +observed base is mapped to its sense transcript. It preserves every observed +base and internal junction, insertion or mismatch. An unmapped end remains +uncompleted. No candidate is returned when no reference bases can be added. +To request one end, use `ends=("five_prime",)` or `ends=("three_prime",)`. + +By default, each end uses the transcript already named by the observation. +To consider alternative 5′ isoforms, pass +`five_prime_transcripts=observation.effect.five_prime_transcript.gene.transcripts`; +`three_prime_transcripts` works the same way for a sense partner. Only isoforms +that agree with all mapped donor runs in genomic coordinates, sequence and +splice path survive. The reference build, gene, contig and strand must agree. +Different transcript pairs remain distinct even when their proteins match. +An explicitly empty list of choices produces no hypotheses for an eligible end. + +Each result has source `varcode_reference_completion`. Its `reference_completion` +evidence records selected transcript IDs, the half-open `observed_span` in +completed cDNA, and `assumed_spans` with their cDNA and reference coordinates. +Assumed segments are labeled `assumed_reference_five_prime` or +`assumed_reference_three_prime`. Original candidate/model evidence and peptide +are retained in `observed_evidence`, `observed_model_evidence` and +`observed_protein_sequence`; `observed_source` identifies the producer. Read +counts and abundance stay nested with the original observation. They provide +no extra support for a completed isoform, and must not be summed across its +alternative completions. Candidate JSON serialization preserves this provenance. + +RNA producers can set `rna_five_prime_complete=True` and/or +`rna_three_prime_complete=True` in candidate or model evidence when they have +actual transcript-end evidence. Completion never extends those ends. False or +None leaves the end eligible. Exacto's ORF-completeness labels and read endpoints +alone do not establish transcript ends; native rows retain their original scope. + +Translation starts only at the retained annotated 5′ initiator, recorded as +`initiation_status="assumed_reference_start"`. There is no search for a new +internal/upstream ATG. The original Exacto ORF is unchanged, even when its +reading frame differs; `observed_orf_frame` reports its frame relative to the +completed CDS when both starts are known. A mapped start can span observed and +assumed segments. Missing/partial start annotation or ambiguous coding bases +leave the protein None. A missing stop gives `protein_completeness="partial_end"`; +adding a reference 3′ end does not guarantee an in-frame stop after a frameshift. +Existing transcript-specific genetic codes and mapped selenocysteine/3′-UTR +rules apply, including explicit uncertain Sec offsets in the hypothesis evidence. + +For example, the observed fusion fragment `CCC GGG AAA TTT` predicts `PGKF` +in its supplied frame. A compatible reference prefix `ATG AAA` and suffix +`CCC TAA` produce a separate `MKPGKFP` hypothesis. Those added bases remain +assumptions. Its sequence-change flags describe that hypothesis; neither a +start-to-stop prediction nor a changed-protein flag establishes observed +initiation, translation, full-length RNA, tumor specificity or abundance. + +This distinction follows the producer's [observed-sequence ORF selection](https://github.com/pirl-unc/exacto/blob/307c08670d5e706734bddf393bcebc84db497f9f/exacto/exacto-translator/src/algorithms/translation.rs), +the [experimental dependence of initiation on sequence context](https://pubmed.ncbi.nlm.nih.gov/2601709/), +and the [3′-UTR requirement for selenocysteine decoding](https://www.nature.com/articles/353273a0). Not every RNA junction produces a coding fusion. In the osteosarcoma regression fixtures, GABBR1 joins sequence upstream of SLC29A1, OTUD7A joins diff --git a/tests/test_reference_completion.py b/tests/test_reference_completion.py new file mode 100644 index 0000000..af8a265 --- /dev/null +++ b/tests/test_reference_completion.py @@ -0,0 +1,376 @@ +"""Completion assumptions never replace an RNA observation or its ORF.""" + +from copy import deepcopy +from dataclasses import replace +from types import SimpleNamespace + +import pytest + +from varcode import EffectCandidate, ReferenceSegment, reference_completion_hypotheses +from varcode.effects import GeneFusion +from varcode.mutant_transcript import _AssembledAllele + +from . import test_exacto_reference_mapping as mapping +from .test_exacto_reference_mapping import load, rows_for, transcript + +context = mapping.context + + +def fragment(context, start=6, end=18, translate=False): + _, first, second = context + return load(context, rows_for([ + (first, start, 9, "match", None), (second, 9, end, "match", None)]), translate=translate) + + +@pytest.mark.parametrize("ends,expected,span,completeness", [ + (("five_prime", "three_prime"), "ATGAAACCCGGGAAATTTCCCTAA", [6, 18], "start_to_stop"), + (("five_prime",), "ATGAAACCCGGGAAATTT", [6, 18], "partial_end"), + (("three_prime",), "CCCGGGAAATTTCCCTAA", [0, 12], "unknown"), +]) +def test_completion_preserves_observed_bases_and_marks_assumptions(context, ends, expected, span, completeness): + candidate = fragment(context, translate=True) + original = candidate.effect.mutant_transcript + evidence = deepcopy(candidate.evidence) + hypothesis, = reference_completion_hypotheses(candidate, ends=ends) + model = hypothesis.effect.mutant_transcript + assert model.cdna_sequence == expected + assert model.cdna_sequence[slice(*span)] == original.cdna_sequence + assert hypothesis.evidence["reference_completion"]["observed_span"] == span + assert hypothesis.evidence["protein_completeness"] == completeness + assert hypothesis.source == model.annotator_name == "varcode_reference_completion" + assert hypothesis.evidence == model.evidence + assert candidate.effect.mutant_transcript is original + assert candidate.evidence == evidence + assert original.cdna_sequence == "CCCGGGAAATTT" + assert original.mutant_protein_sequence == "PGKF" + assert model.mutant_protein_sequence == { + "start_to_stop": "MKPGKFP", "partial_end": "MKPGKF", "unknown": None}[completeness] + if completeness != "unknown": + assert hypothesis.evidence["observed_orf_frame"] == "in_frame" + assert hypothesis.evidence["initiation_status"] == "assumed_reference_start" + for assumed in hypothesis.evidence["reference_completion"]["assumed_spans"]: + assert assumed["end"] <= span[0] or assumed["start"] >= span[1] + assert "".join(s.source.sequence[s.start:s.end] for s in model.reference_segments) == expected + + +def test_no_added_bases_or_requested_ends_means_no_duplicate(context): + _, first, second = context + candidate = load(context, rows_for([(first, 0, 9, "match", None), + (second, 9, 24, "match", None)])) + assert reference_completion_hypotheses(candidate) == () + assert reference_completion_hypotheses(fragment(context), ends=()) == () + + +@pytest.mark.parametrize("end", ["five_prime", "three_prime"]) +def test_real_rna_end_evidence_prevents_extension_but_read_ends_do_not(context, end): + candidate = fragment(context) + candidate = replace(candidate, evidence=dict(candidate.evidence, **{ + "rna_" + end + "_complete": True, "read_count": 12, "read_end": 200})) + result, = reference_completion_hypotheses(candidate) + assumed = result.evidence["reference_completion"]["assumed_spans"] + assert len(assumed) == 1 + assert assumed[0]["label"] != "assumed_reference_" + end + assert "read_count" not in result.evidence + assert result.evidence["observed_evidence"]["read_count"] == 12 + result.evidence["observed_evidence"]["exacto_structure"][0]["sequence"] = "changed" + assert candidate.evidence["exacto_structure"][0]["sequence"] != "changed" + + +def test_model_end_evidence_is_respected(context): + candidate = fragment(context) + model = candidate.effect.mutant_transcript + candidate.effect.mutant_transcript = replace(model, evidence=dict( + model.evidence, rna_five_prime_complete=True, rna_three_prime_complete=True)) + assert reference_completion_hypotheses(candidate) == () + + +def alternate(source, identifier, prepend="", append=""): + """Explicit extra terminal exons; preserve the observed exon coordinates.""" + positions = list(source.positions) + if prepend: + positions = (list(range(20, 20 + len(prepend))) if source.strand == "+" else + list(range(400 + len(prepend) - 1, 399, -1))) + positions + if append: + positions += (list(range(400, 400 + len(append))) if source.strand == "+" else + list(range(20 + len(append) - 1, 19, -1))) + return SimpleNamespace(**dict(vars(source), id=identifier, + sequence=prepend + source.sequence + append, + positions=positions, spliced_offset=positions.index, + start_codon_spliced_offsets=[p + len(prepend) for p in source.start_codon_spliced_offsets], + stop_codon_spliced_offsets=[p + len(prepend) for p in source.stop_codon_spliced_offsets])) + + +def test_isoform_pairs_remain_distinct_even_with_identical_proteins(context): + _, first, second = context + upstream = alternate(first, "T1_alt", prepend="CCA") + downstream = alternate(second, "T2_alt", append="GGG") + results = reference_completion_hypotheses( + fragment(context), five_prime_transcripts=[first, upstream, upstream], + three_prime_transcripts=[second, downstream]) + assert len(results) == 4 + assert {r.effect.mutant_protein_sequence for r in results} == {"MKPGKFP"} + assert {(r.effect.transcript.id, r.effect.partner_transcript.id) for r in results} == { + ("T1", "T2"), ("T1_alt", "T2"), ("T1", "T2_alt"), ("T1_alt", "T2_alt")} + for result in results: + model = result.effect.mutant_transcript + assert "".join(s.source.sequence[s.start:s.end] for s in model.reference_segments) == model.cdna_sequence + assert result.evidence["reference_completion"]["observed_span"] == ( + [9, 21] if result.effect.transcript.id == "T1_alt" else [6, 18]) + + +@pytest.mark.parametrize("bad", ["gene", "strand", "contig", "genome", "sequence", "splice", "missing"]) +def test_incompatible_alternatives_are_not_completed(context, bad): + _, first, _ = context + other = alternate(first, "bad") + if bad in ("gene", "strand", "contig", "genome"): + setattr(other, "gene_id" if bad == "gene" else bad, "wrong") + elif bad == "sequence": + other.sequence = first.sequence[:7] + "T" + first.sequence[8:] + elif bad == "splice": + # Add an exon between two observed runs; terminal matching alone + # would accept it, but the observed splice path rules it out. + other.spliced_offset = lambda p: first.spliced_offset(p) + (3 if first.spliced_offset(p) >= 7 else 0) + else: + other.sequence = None + assert reference_completion_hypotheses(fragment(context), five_prime_transcripts=[other]) == () + + +@pytest.mark.parametrize("side", ["five_prime", "three_prime"]) +def test_unmapped_terminal_bases_prevent_only_that_end_completion(context, side): + candidate = fragment(context) + model = candidate.effect.mutant_transcript + segments = list(model.reference_segments) + index = 0 if side == "five_prime" else -1 + old = segments[index] + sequence = old.source.sequence[old.start:old.end] + segments[index] = ReferenceSegment(_AssembledAllele(sequence), 0, len(sequence)) + candidate.effect.mutant_transcript = replace(model, reference_segments=tuple(segments)) + result, = reference_completion_hypotheses(candidate) + spans = result.evidence["reference_completion"]["assumed_spans"] + assert len(spans) == 1 and spans[0]["label"] != "assumed_reference_" + side + + +def test_split_annotated_start_codon_is_retained(context): + candidate = fragment(context, start=1) + result, = reference_completion_hypotheses(candidate) + assert result.effect.mutant_protein_sequence == "MKPGKFP" + assert result.evidence["cds_start"] == 0 + + +def test_out_of_frame_internal_atg_is_a_separate_producer_hypothesis(context): + from varcode import load_exacto_fusions + from .test_exacto_fusions import link, table + _, first, second = context + first.sequence = first.coding_sequence = "ATGAAACATGAATTTCCCAAATAA" + first.protein_sequence = "MKHEFPK" + rows = rows_for([(first, 6, 12, "match", None), (second, 12, 18, "match", None)]) + sequence = "CATGAAAAATTT" + primary = [] + origins = [r["index"] for r in rows for _ in r["sequence"]] + for position in range(1, 10): + offset = position - 1 + primary.append(dict(peptide_id="internal", primary_structure_index=str(offset), + type="base", amino_acid="MKN"[offset // 3], amino_acid_index=str(offset // 3), + codon_index=str(offset % 3), nucleotide=sequence[position], transcript_model_id="1", + reference_transcript_ids="T1,T2", transcript_structure_index=origins[position], + read_start=str(position), read_end=str(position))) + candidate, = load_exacto_fusions(table(rows), table([link()]), + variants_by_id={"D1": context[0]}, primary_structures_path=table(primary)).candidates + result, = reference_completion_hypotheses(candidate) + assert candidate.effect.mutant_protein_sequence == "MKN" + assert result.evidence["observed_protein_sequence"] == "MKN" + assert result.effect.mutant_protein_sequence == "MKHEKFP" + assert result.evidence["observed_orf_frame"] == "out_of_frame" + + +def test_unchanged_completed_sequence_does_not_reclassify_observed_fragment(context): + _, first, second = context + first.sequence = first.coding_sequence = second.sequence + first.protein_sequence = second.protein_sequence + candidate = fragment(context, translate=True) + assert candidate.effect.modifies_protein_sequence is None + result, = reference_completion_hypotheses(candidate) + assert result.effect.modifies_protein_sequence is False + assert result.effect.modifies_coding_sequence is False + assert candidate.effect.modifies_protein_sequence is None + + +def test_inframe_internal_atg_suffix_remains_an_observation_alongside_completion(context): + _, first, second = context + for t in (first, second): + t.sequence = t.coding_sequence = "ATGAAAATGGGGTTTCCCAAATAA" + t.protein_sequence = "MKMGFPK" + candidate = fragment(context, end=24, translate=True) + assert candidate.effect.mutant_protein_sequence == "MGFPK" + assert candidate.evidence["protein_completeness"] == "start_to_stop" + assert candidate.effect.modifies_protein_sequence is None + result, = reference_completion_hypotheses(candidate) + assert result.effect.mutant_protein_sequence == "MKMGFPK" + assert result.effect.modifies_protein_sequence is False + assert result.evidence["observed_orf_frame"] == "in_frame" + assert candidate.effect.mutant_protein_sequence == "MGFPK" + assert candidate.effect.modifies_protein_sequence is None + + +def test_observed_insertion_is_never_filled_or_removed(context): + _, first, second = context + candidate = load(context, rows_for([(first, 6, 9, "match", None), + (first, 9, 9, "insertion", "A"), (first, 9, 15, "match", None), + (second, 15, 18, "match", None)])) + result, = reference_completion_hypotheses(candidate) + assert result.effect.mutant_transcript.cdna_sequence == "ATGAAACCCAGGGTTTTTTCCCTAA" + # ATG AAA CCC AGG GTT TTT TCC CTA, followed by one incomplete base. + assert result.effect.mutant_protein_sequence == "MKPRVFSL" + assert result.evidence["protein_completeness"] == "partial_end" + assert result.effect.modifies_protein_sequence is True + + +@pytest.mark.parametrize("mutation", ["layout", "end_flag", "recursive", "ends"]) +def test_invalid_or_recursive_completion_is_rejected(context, mutation): + candidate = fragment(context) + options = {} + if mutation == "layout": + candidate.effect.mutant_transcript = replace(candidate.effect.mutant_transcript, cdna_sequence="ACGT") + elif mutation == "end_flag": + candidate = replace(candidate, evidence=dict(candidate.evidence, rna_five_prime_complete="yes")) + elif mutation == "recursive": + candidate, = reference_completion_hypotheses(candidate) + else: + options["ends"] = "five_prime" + with pytest.raises(ValueError): + reference_completion_hypotheses(candidate, **options) + + +def test_selenocysteine_uses_retained_reference_utr_in_the_hypothesis(): + # UGA is Sec only under the existing mapped-UTR policy. The assumed UTR + # is visible as reference completion, never claimed as observed coverage. + first = transcript("T1", "+", "1", "ATGAAACCCGGGTTTCCCAAATAA") + second = transcript("T2", "+", "2", "ATGTGAAAATAACCC") + second.protein_sequence = "MUK" + second.coding_sequence = second.sequence[:12] + second.stop_codon_spliced_offsets = [9, 10, 11] + from varcode import MutantTranscript + segments = (ReferenceSegment(first, 6, 9), ReferenceSegment(second, 3, 9)) + model = MutantTranscript(reference_transcript=first, reference_segments=segments, + cdna_sequence="CCCTGAAAA") + candidate = EffectCandidate(GeneFusion(None, first, second, mutant_transcript=model), source="rna") + result, = reference_completion_hypotheses(candidate) + assert result.effect.mutant_protein_sequence == "MKPUK" + assert result.evidence["selenocysteine_decoded_offsets"] == [9] + result, = reference_completion_hypotheses(candidate, ends=("five_prime",)) + assert result.effect.mutant_protein_sequence == "MKP" + partial_utr = replace(model, + reference_segments=(segments[0], ReferenceSegment(second, 3, 13)), + cdna_sequence="CCCTGAAAATAAC") + partial = replace(candidate, effect=GeneFusion(None, first, second, mutant_transcript=partial_utr)) + result, = reference_completion_hypotheses(partial, ends=("five_prime",)) + assert result.effect.mutant_protein_sequence == "MKPUK" + assert result.evidence["selenocysteine_uncertain_offsets"] == [9] + + +def test_unannotated_partner_still_allows_five_prime_completion(context): + from varcode.effects import TranslocationToIntergenic + candidate = fragment(context) + model = candidate.effect.mutant_transcript + candidate = replace(candidate, effect=TranslocationToIntergenic( + candidate.effect.variant, context[1], mutant_transcript=model)) + result, = reference_completion_hypotheses(candidate) + assert isinstance(result.effect, TranslocationToIntergenic) + assert result.effect.mutant_transcript.cdna_sequence == "ATGAAACCCGGGAAATTT" + assert result.evidence["protein_completeness"] == "partial_end" + assert result.evidence["reference_completion"]["three_prime_transcript_id"] is None + + +def test_ambiguous_observed_bases_remain_ambiguous(context): + _, first, second = context + candidate = load(context, rows_for([(first, 6, 9, "match", None), + (first, 9, 12, "mismatch", "NNN"), (second, 12, 18, "match", None)])) + result, = reference_completion_hypotheses(candidate) + assert result.effect.mutant_transcript.cdna_sequence == "ATGAAACCCNNNAAATTTCCCTAA" + assert result.effect.mutant_protein_sequence is None + assert result.evidence["protein_status"] == "unresolved_ambiguous_coding_sequence" + assert result.effect.modifies_protein_sequence is None + + +def test_partial_end_request_accepts_an_iterator(context): + result, = reference_completion_hypotheses(fragment(context), ends=iter(["five_prime"])) + assert result.effect.mutant_protein_sequence == "MKPGKF" + + +def test_original_three_prime_partner_perspective_uses_the_five_prime_start(context): + candidate = fragment(context) + _, first, second = context + candidate = replace(candidate, effect=GeneFusion(context[0], second, first, + mutant_transcript=candidate.effect.mutant_transcript, + five_prime_transcript=first, three_prime_transcript=second)) + result, = reference_completion_hypotheses(candidate) + assert result.effect.five_prime_transcript is first + assert result.effect.three_prime_transcript is second + assert result.effect.mutant_protein_sequence == "MKPGKFP" + + +def test_isoform_with_a_different_upstream_cds_generates_its_own_protein(context): + _, first, _ = context + upstream = alternate(first, "T1_long", prepend="ATG") + upstream.start_codon_spliced_offsets = [0, 1, 2] + upstream.coding_sequence = upstream.sequence + upstream.protein_sequence = "M" + first.protein_sequence + results = reference_completion_hypotheses(fragment(context), five_prime_transcripts=[first, upstream]) + assert [r.effect.mutant_protein_sequence for r in results] == ["MKPGKFP", "MMKPGKFP"] + + +@pytest.mark.parametrize("starts", [[], [0, 1], [0, 2, 3]]) +def test_missing_or_partial_start_annotation_does_not_guess_an_orf(context, starts): + candidate = fragment(context) + context[1].start_codon_spliced_offsets = starts + result, = reference_completion_hypotheses(candidate) + assert result.effect.mutant_protein_sequence is None + assert result.evidence["protein_completeness"] == "unknown" + + +def test_pyensembl_missing_start_annotation_is_an_unresolved_protein(context): + class Noncoding(SimpleNamespace): + @property + def start_codon_spliced_offsets(self): + raise ValueError("No start_codon features found") + + _, first, _ = context + candidate = fragment(context) + noncoding = Noncoding(**{k: v for k, v in vars(first).items() if k != "start_codon_spliced_offsets"}) + result, = reference_completion_hypotheses(candidate, five_prime_transcripts=[noncoding]) + assert result.effect.mutant_protein_sequence is None + + +def test_observed_stop_is_not_overridden_by_three_prime_completion(context): + _, first, second = context + candidate = load(context, rows_for([(first, 6, 9, "match", None), + (first, 9, 12, "mismatch", "TAA"), (second, 12, 18, "match", None)])) + result, = reference_completion_hypotheses(candidate) + assert result.effect.mutant_protein_sequence == "MKP" + assert result.evidence["cds_end"] == 12 + assert result.evidence["protein_completeness"] == "start_to_stop" + + +def test_real_transcripts_and_candidate_json_roundtrip(): + from pyensembl import cached_release + from varcode import MutantTranscript, StructuralVariant + genome = cached_release(81) + first = genome.transcript_by_id("ENST00000003084") + second = genome.transcript_by_id("ENST00000357654") + left = min(first.start_codon_spliced_offsets) + 6 + right = min(second.start_codon_spliced_offsets) + 9 + observed = first.sequence[left:left + 9] + second.sequence[right:right + 9] + model = MutantTranscript(reference_transcript=first, + reference_segments=(ReferenceSegment(first, left, left + 9), ReferenceSegment(second, right, right + 9)), + cdna_sequence=observed, evidence={"read_count": 7}) + variant = StructuralVariant(first.contig, first.start, "BND", genome=genome, + mate_contig=second.contig, mate_start=second.start) + candidate = EffectCandidate(GeneFusion(variant, first, second, mutant_transcript=model), source="rna") + result, = reference_completion_hypotheses(candidate) + restored = EffectCandidate.from_json(result.to_json()) + assert restored.source == result.source + assert restored.evidence == result.evidence + assert restored.effect.mutant_protein_sequence == result.effect.mutant_protein_sequence + completed = restored.effect.mutant_transcript + assert completed.cdna_sequence == first.sequence[:left + 9] + second.sequence[right:] + assert completed.cdna_sequence[slice(*restored.evidence["reference_completion"]["observed_span"])] == observed diff --git a/varcode/__init__.py b/varcode/__init__.py index 49fa928..9888fa6 100644 --- a/varcode/__init__.py +++ b/varcode/__init__.py @@ -79,6 +79,7 @@ make_rna_outcome, ) from .exacto import load_exacto_fusions +from .reference_completion import reference_completion_hypotheses from .realized_effects import predict_realized_effect from .transcript_model import predict_transcript_model_effect from .splice_outcomes import ( @@ -197,6 +198,7 @@ "make_rna_outcome", "make_fusion_outcome", "load_exacto_fusions", + "reference_completion_hypotheses", "EffectAnnotator", "FastEffectAnnotator", diff --git a/varcode/reference_completion.py b/varcode/reference_completion.py new file mode 100644 index 0000000..6ec95cc --- /dev/null +++ b/varcode/reference_completion.py @@ -0,0 +1,263 @@ +"""Explicit reference-completion hypotheses for mapped RNA observations.""" + +from copy import deepcopy +from dataclasses import replace +from itertools import product + +from .effect_candidates import EffectCandidate +from .effects.codon_tables import codon_table_for_transcript, translate_sequence +from .effects.effect_classes import GeneFusion, TranslocationToIntergenic +from .effects.selenocysteine import segment_selenocysteine +from .effects.sequence_change import _retained_start +from .mutant_transcript import ReferenceSegment +from .nucleotides import reverse_complement + + +_SOURCE = "varcode_reference_completion" +_ENDS = ("five_prime", "three_prime") + + +def _genomic_pieces(transcript, start, end): + """Inclusive genomic pieces for a cDNA interval, in transcript order.""" + offset = 0 + for lo, hi in sorted(transcript.exon_intervals, reverse=transcript.strand == "-"): + length = hi - lo + 1 + a, b = max(start, offset) - offset, min(end, offset + length) - offset + if a < b: + yield ((lo + a, lo + b - 1) if transcript.strand == "+" + else (hi - a, hi - b + 1)) + offset += length + + +def _remap(segments, source, target): + """A compatible isoform keeps one cDNA shift across all mapped donor runs.""" + if any(getattr(source, key, None) != getattr(target, key, None) + for key in ("gene_id", "contig", "strand", "genome")): + return None + sequence = getattr(target, "sequence", None) + if not sequence: + return None + shift = None + mapped = [] + for segment in segments: + if segment.source != source: + mapped.append(segment) + continue + if segment.strand != "+": + return None + if target == source: + start = segment.start + else: + offsets = [] + try: + for first, last in _genomic_pieces(source, segment.start, segment.end): + a, b = target.spliced_offset(first), target.spliced_offset(last) + if b - a != abs(last - first): + return None + if offsets and a != offsets[-1][1] + 1: + return None + offsets.append((a, b)) + except (AttributeError, ValueError): + return None + if not offsets or offsets[-1][1] - offsets[0][0] + 1 != segment.length: + return None + start = offsets[0][0] + delta = start - segment.start + if shift is not None and shift != delta: + return None + shift = delta + end = start + segment.length + if (start < 0 or end > len(sequence) + or sequence[start:end].upper() != source.sequence[segment.start:segment.end].upper() + or set(sequence[start:end].upper()) - set("ACGT")): + return None + mapped.append(replace(segment, source=target, start=start, end=end)) + return tuple(mapped) + + +def _end_options(model, source, choices, end, enabled): + segments = model.reference_segments + terminal = segments[0 if end == "five_prime" else -1] + if not enabled or source is None or terminal.source != source or terminal.strand != "+": + return [(source, segments, None)] + options, seen = [], set() + for transcript in (source,) if choices is None else choices: + identifier = getattr(transcript, "id", None) + if identifier is None: + raise ValueError("Completion transcripts must have an id") + if identifier in seen: + continue + seen.add(identifier) + mapped = _remap(segments, source, transcript) + if mapped is None: + continue + edge = mapped[0 if end == "five_prime" else -1] + start, stop = ((0, edge.start) if end == "five_prime" + else (edge.end, len(transcript.sequence))) + extension = (ReferenceSegment(transcript, start, stop, + label="assumed_reference_" + end) + if start < stop else None) + options.append((transcript, mapped, extension)) + return options + + +def _translate(model, transcript): + """Use the retained annotated start, never an ORF search over added bases.""" + evidence = dict(model.evidence) + evidence.update(protein_status="not_determined", protein_completeness="unknown") + try: + starts = sorted(getattr(transcript, "start_codon_spliced_offsets", ())) + except ValueError: + # PyEnsembl raises when start-codon annotation is absent or partial. + starts = [] + if len(starts) != 3 or starts != list(range(starts[0], starts[0] + 3)): + return replace(model, evidence=evidence) + start = _retained_start(model, transcript) + sequence = model.cdna_sequence + table = codon_table_for_transcript(transcript) + if start is None or sequence[start:start + 3] not in table.start_codons: + return replace(model, evidence=evidence) + decoded, uncertain, _ = segment_selenocysteine(model) + sec = decoded | uncertain + stop = next((i for i in range(start + 3, len(sequence) - 2, 3) + if sequence[i:i + 3] in table.stop_codons and i not in sec), None) + end = len(sequence) if stop is None else stop + 3 + coding = sequence[start:end - (end - start) % 3] + if set(coding) - set("ACGT"): + evidence["protein_status"] = "unresolved_ambiguous_coding_sequence" + return replace(model, evidence=evidence) + protein = translate_sequence(coding, codon_table=table, + selenocysteine={p - start for p in sec}) + protein = "M" + protein[1:] + observed_start = evidence["observed_model_evidence"].get("cds_start") + relation = "unknown" + if isinstance(observed_start, int) and not isinstance(observed_start, bool): + delta = observed_start + evidence["reference_completion"]["observed_span"][0] - start + relation = "in_frame" if delta % 3 == 0 else "out_of_frame" + evidence.update( + cds_start=start, cds_end=end, initiation_status="assumed_reference_start", + protein_status="predicted_from_reference_completion", + protein_completeness="partial_end" if stop is None else "start_to_stop", + protein_translation_table=table.ncbi_id, + observed_orf_frame=relation, + selenocysteine_decoded_offsets=sorted(decoded), + selenocysteine_uncertain_offsets=sorted(uncertain)) + return replace(model, mutant_protein_sequence=protein, evidence=evidence) + + +def reference_completion_hypotheses( + candidate, *, five_prime_transcripts=None, three_prime_transcripts=None, + ends=("five_prime", "three_prime")): + """Complete mapped RNA ends as separate, explicitly assumed hypotheses. + + Parameters + ---------- + candidate : EffectCandidate + Observed GeneFusion or TranslocationToIntergenic with materialized + cDNA and reference segments. The candidate is never modified. + five_prime_transcripts, three_prime_transcripts : iterable or None + Reference isoforms to consider at each end. None uses that end's + annotated transcript. Alternatives must agree with the observed mapped + donor runs in genomic coordinates, sequence and splice path. An empty + iterable yields no hypotheses when that end can be completed. + ends : iterable of str + Ends to complete: ``five_prime``, ``three_prime``, or both (default). + + Returns + ------- + tuple of EffectCandidate + New candidates with source ``varcode_reference_completion``. At least + one reference base is added per hypothesis. Unmapped ends are left + alone. Explicit ``rna_five_prime_complete`` / ``rna_three_prime_complete`` + evidence set to True prevents extension of that end. ORF completeness + and read ends do not establish RNA ends. + + Evidence records half-open observed/assumed spans in completed cDNA, + selected transcripts, and nested original evidence/protein. Read counts + are not promoted to support for the inferred sequence. Translation uses + the retained annotated 5' start and existing selenocysteine rules; it + does not establish biological initiation, translation or expression. + Return these alongside observations, never as RNA-observed outcomes. + + Raises + ------ + ValueError + For unsupported effects, malformed layouts, invalid end evidence, or + an already reference-completed input. + """ + if isinstance(ends, str): + raise ValueError("ends must contain five_prime and/or three_prime") + ends = set(ends) + if ends - set(_ENDS): + raise ValueError("ends must contain five_prime and/or three_prime") + effect = candidate.effect + if not isinstance(effect, (GeneFusion, TranslocationToIntergenic)): + raise ValueError("Reference completion requires an observed fusion outcome") + model = effect.mutant_transcript + if (candidate.source == _SOURCE or (model and (model.evidence or {}).get( + "sequence_status") == "reference_completion_hypothesis")): + raise ValueError("Cannot complete an already reference-completed hypothesis") + if (model is None or not model.cdna_sequence or not model.reference_segments + or model.edits): + raise ValueError("Reference completion requires a materialized segment layout without edits") + rendered = [] + for segment in model.reference_segments: + sequence = getattr(segment.source, "sequence", None) + if sequence is None or segment.length <= 0 or segment.end > len(sequence): + raise ValueError("Invalid observed reference segment") + piece = sequence[segment.start:segment.end] + rendered.append(piece if segment.strand == "+" else reverse_complement(piece)) + if "".join(rendered).upper() != model.cdna_sequence.upper(): + raise ValueError("Reference segments do not render the observed cDNA") + enabled = {} + for end in _ENDS: + flags = [(evidence or {}).get("rna_" + end + "_complete") + for evidence in (candidate.evidence, model.evidence)] + if any(flag is not None and not isinstance(flag, bool) for flag in flags): + raise ValueError("RNA end completeness must be bool or None") + enabled[end] = end in ends and True not in flags + first = getattr(effect, "five_prime_transcript", effect.transcript) + last = getattr(effect, "three_prime_transcript", None) + left = _end_options(model, first, five_prime_transcripts, "five_prime", enabled["five_prime"]) + right = _end_options(model, last, three_prime_transcripts, "three_prime", enabled["three_prime"]) + hypotheses = [] + for (five, left_segments, prefix), (three, right_segments, suffix) in product(left, right): + if prefix is None and suffix is None: + continue + # Both remappings apply to different genes, preserving all opaque bases. + observed = tuple(r if s.source == last else l for s, l, r in zip( + model.reference_segments, left_segments, right_segments)) + segments = ((prefix,) if prefix else ()) + observed + ((suffix,) if suffix else ()) + offset = prefix.length if prefix else 0 + observed_end = offset + len(model.cdna_sequence) + assumed = [] + for segment, start, end in ((prefix, 0, offset), + (suffix, observed_end, observed_end + (suffix.length if suffix else 0))): + if segment: + assumed.append(dict(start=start, end=end, transcript_id=segment.source.id, + reference_start=segment.start, reference_end=segment.end, + label=segment.label)) + evidence = dict( + sequence_status="reference_completion_hypothesis", + observed_source=candidate.source, + observed_evidence=deepcopy(dict(candidate.evidence)), + observed_model_evidence=deepcopy(dict(model.evidence or {})), + observed_protein_sequence=model.mutant_protein_sequence, + reference_completion=dict( + observed_span=[offset, observed_end], assumed_spans=assumed, + five_prime_transcript_id=five.id, + three_prime_transcript_id=three.id if three else None)) + sequence = ((prefix.source.sequence[prefix.start:prefix.end] if prefix else "") + + model.cdna_sequence + + (suffix.source.sequence[suffix.start:suffix.end] if suffix else "")) + completed = _translate(replace( + model, reference_transcript=five, reference_segments=segments, + cdna_sequence=sequence.upper(), mutant_protein_sequence=None, + annotator_name=_SOURCE, evidence=evidence), five) + if isinstance(effect, GeneFusion): + result = GeneFusion(effect.variant, five, three, mutant_transcript=completed) + else: + result = TranslocationToIntergenic(effect.variant, five, mutant_transcript=completed) + result.candidate_evidence = completed.evidence + hypotheses.append(EffectCandidate(result, source=_SOURCE, evidence=completed.evidence)) + return tuple(hypotheses) diff --git a/varcode/version.py b/varcode/version.py index 3afb665..e6f9023 100644 --- a/varcode/version.py +++ b/varcode/version.py @@ -1 +1 @@ -__version__ = "10.9.0" +__version__ = "10.10.0"