RNA and transcripts API

See Transcript models for sequence access and RNA imports for observed structures and Exacto predictions.

Mutant transcripts

varcode.MutantTranscript(reference_transcript: Optional[object] = None, edits: Tuple[TranscriptEdit, ...] = tuple(), reference_segments: Optional[Tuple[ReferenceSegment, ...]] = None, cdna_sequence: Optional[str] = None, mutant_protein_sequence: Optional[str] = None, annotator_name: str = 'unknown', evidence: Optional[dict] = None) dataclass

Bases: DataclassSerializable

A reference transcript (or assembled set of reference segments) with zero or more variant-derived edits applied, optionally carrying the mutated cDNA and protein sequences.

Producers (the protein-diff annotator, RNA-evidence importers, the splice-outcomes rewrite, germline-aware annotation, structural-variant annotators) construct this once per (transcript-or-segments, variant-set, context) and hand it to downstream consumers. Each consumer reads the fields it cares about — edits for provenance, cdna_sequence / mutant_protein_sequence for protein-level analysis.

Two shapes:

• Point-variant shape (reference_transcript is set, reference_segments is None) — the mutant is derived from a single reference transcript by applying zero or more :class:TranscriptEdit objects. This is the shape used by the protein-diff annotator for SNVs, MNVs, and simple indels.

• Structural-variant shape (reference_segments is set, reference_transcript may be None or the primary / 5'-partner transcript) — the mutant is assembled by concatenating :class:ReferenceSegment slices in order. A gene fusion is two segments from two transcripts; a translocation to intergenic is one transcript segment plus a genomic-interval segment; an inversion is three forward/ reverse/forward segments. edits may still be populated for point-variant edits layered on top of the assembled segments, but typically an SV carries no extra edits.

Sequence fields are Optional[str] because not every producer computes them eagerly. Callers that require the protein check or compute it themselves; the protein-diff annotator guarantees it for point variants.

Forward-looking hooks (not implemented here; documented so new integrations know where to plug in):

• Personalized / full-genome reference — pass a :class:ReferenceSegment whose source is a patient-specific contig object. varcode's translation logic reads source.sequence; it doesn't care whether that's GRCh38 or a custom assembly.
• Long-read resolution — when an SV has an :attr:alt_assembly on the :class:StructuralVariant, the SV annotator can wrap that sequence as a single synthetic segment.
• SV outcomes ambiguity — a translocation producing many candidate ORFs that only RNA can resolve should return List[MutantTranscript] or wrap it in a :class:MultiOutcomeEffect per #299, each with its own evidence dict capturing the disambiguator.

reference_transcript: Optional[object] = None class-attribute instance-attribute

The :class:pyensembl.Transcript this mutant is derived from (point-variant shape), or the primary / 5'-partner transcript (SV shape). None when the SV has no canonical primary transcript (e.g. intergenic-to-intergenic BNDs). Not typed tightly here so :mod:pyensembl isn't a hard import dependency for anyone who just wants the dataclass.

edits: Tuple[TranscriptEdit, ...] = field(default_factory=tuple) class-attribute instance-attribute

Edits applied to produce this mutant, sorted by :attr:TranscriptEdit.cdna_start. Empty tuple means the mutant is identical to the reference (or, for SV shape, the assembled segments carry the rearrangement directly without further point-level edits).

reference_segments: Optional[Tuple[ReferenceSegment, ...]] = None class-attribute instance-attribute

Ordered tuple of :class:ReferenceSegment objects that, when concatenated in order (applying reverse-complement to - strand segments), produce the mutant cDNA. None for the point-variant shape; a fusion's segments would be (5p_partner_segment, 3p_partner_segment). Coordinates are in each segment's own reference system. A partial structural model may describe only a retained reference fragment: its producer identifies this in evidence and leaves the full cdna_sequence unknown. Concatenating those partial segments does not establish a full allele.

cdna_sequence: Optional[str] = None class-attribute instance-attribute

The mutated spliced mRNA, when computed. None if the producer hasn't materialized it yet.

mutant_protein_sequence: Optional[str] = None class-attribute instance-attribute

The translated mutant protein, stopping at the first stop codon. None if not yet translated, or if the edit set doesn't produce a coherent ORF (e.g. start-codon loss). Callers that need a guaranteed-present protein should use the protein-diff annotator once it lands.

annotator_name: str = 'unknown' class-attribute instance-attribute

Name of the :class:EffectAnnotator (or other producer) that created this MutantTranscript. Used as provenance in serialization and for A/B comparisons.

evidence: Optional[dict] = None class-attribute instance-attribute

Optional producer-specific evidence (RNA read counts, Isovar fragment ids, SpliceAI scores, long-read assembly metadata). Shape is annotator-specific and not part of the stable contract; consumers that care about a particular evidence shape should type-check it at the call site.

is_identical_to_reference: bool property

True if no edits were applied AND there are no reference-rearranging segments. Does NOT check cdna_sequence / mutant_protein_sequence — a producer can legitimately carry an identical sequence with zero edits and a single identity segment.

is_structural: bool property

True when this mutant was assembled from :attr:reference_segments (SV shape) rather than applying :attr:edits to a single reference transcript.

total_length_delta: int property

Sum of :attr:TranscriptEdit.length_delta across all edits — how much longer or shorter the mutant cDNA is than the reference (point-variant shape). For SV shape, returns 0; the length of an assembled cDNA is the sum of segment lengths, not a delta against a single reference.

from_sequence(sequence, *, reference_transcript=None, mutant_protein_sequence=None, annotator_name='unknown', evidence=None, label='observed_sequence') classmethod

Wrap an oriented external sequence without reconstructing it.

Sequence is already in transcript 5'-to-3' order. This constructor does not infer an ORF, splice structure, or sequence completeness. Producers must describe partial observations in evidence.

Source code in varcode/mutant_transcript.py

@classmethod
def from_sequence(cls, sequence, *, reference_transcript=None,
                  mutant_protein_sequence=None, annotator_name="unknown",
                  evidence=None, label="observed_sequence"):
    """Wrap an oriented external sequence without reconstructing it.

    Sequence is already in transcript 5'-to-3' order. This constructor
    does not infer an ORF, splice structure, or sequence completeness.
    Producers must describe partial observations in ``evidence``.
    """
    return cls(
        reference_transcript=reference_transcript,
        reference_segments=(ReferenceSegment(
            source=_AssembledAllele(sequence), start=0, end=len(sequence),
            label=label),),
        cdna_sequence=sequence,
        mutant_protein_sequence=mutant_protein_sequence,
        annotator_name=annotator_name,
        evidence=evidence)

varcode.apply_variant_to_transcript(variant, transcript)

Construct a :class:MutantTranscript by applying variant to transcript's spliced cDNA.

Returns a :class:MutantTranscript whose cdna_sequence is populated, plus mutant_protein_sequence when the variant lies after the start codon (so translation from the canonical start is well-defined). The codon table is selected from the transcript's contig — mitochondrial transcripts use NCBI table 2 automatically (see :func:varcode.effects.codon_tables.codon_table_for_transcript).

Returns None when the variant can't be cleanly applied:

• Transcript is not protein-coding or is incomplete.
• Variant doesn't overlap the transcript at all.
• Variant spans more than one exon (splice-junction-crossing variants need the splice-aware path; not handled here).
• Reference allele doesn't match the transcript's cDNA at the computed offset.

Callers that get None should fall back to the fast :class:EffectAnnotator. The forthcoming protein-diff annotator layers effect classification on top of this builder.

Source code in varcode/mutant_transcript.py

def apply_variant_to_transcript(variant, transcript):
    """Construct a :class:`MutantTranscript` by applying ``variant``
    to ``transcript``'s spliced cDNA.

    Returns a :class:`MutantTranscript` whose ``cdna_sequence`` is
    populated, plus ``mutant_protein_sequence`` when the variant
    lies after the start codon (so translation from the canonical
    start is well-defined). The codon table is selected from the
    transcript's contig — mitochondrial transcripts use NCBI table
    2 automatically (see :func:`varcode.effects.codon_tables.codon_table_for_transcript`).

    Returns ``None`` when the variant can't be cleanly applied:

    * Transcript is not protein-coding or is incomplete.
    * Variant doesn't overlap the transcript at all.
    * Variant spans more than one exon (splice-junction-crossing
      variants need the splice-aware path; not handled here).
    * Reference allele doesn't match the transcript's cDNA at the
      computed offset.

    Callers that get ``None`` should fall back to the fast
    :class:`EffectAnnotator`. The forthcoming protein-diff annotator
    layers effect classification on top of this builder.
    """
    from pyensembl import Transcript

    if not isinstance(transcript, Transcript):
        return None
    if not transcript.is_protein_coding or not transcript.complete:
        return None

    full_sequence = str(transcript.sequence)
    resolved = _resolve_variant_edit(variant, transcript, full_sequence)
    if resolved is None:
        return None
    edit, cdna_offset = resolved
    mutant_cdna = (
        full_sequence[:edit.cdna_start]
        + edit.alt_bases
        + full_sequence[edit.cdna_end:]
    )
    mutant_protein = None
    cds_start = min(transcript.start_codon_spliced_offsets)
    if cdna_offset >= cds_start:
        mutant_protein = _translate_from_cds(mutant_cdna, transcript, (edit,))
    return MutantTranscript(
        reference_transcript=transcript,
        edits=(edit,),
        cdna_sequence=mutant_cdna,
        mutant_protein_sequence=mutant_protein,
        annotator_name="protein_diff",
    )

varcode.apply_variants_to_transcript(variants, transcript)

Apply a list of variants to a single transcript, yielding one :class:MutantTranscript that carries all the resulting edits (#269). Used for haplotype-aware joint effect prediction — cis variants on the same transcript become one combined mutant rather than N independent per-variant mutants.

Edits are applied in cDNA-coordinate order (highest offset first, so earlier offsets aren't shifted) to transcript's spliced cDNA. Returns None when any of the usual single-variant preconditions fail (non-coding, incomplete, etc.), or when the provided variants conflict — i.e. claim to edit overlapping cDNA ranges. The caller is responsible for falling back to per-variant effects in that case.

mutant_protein_sequence is populated when at least one edit lands after the canonical CDS start; the joint cDNA is translated from there to the first stop.

Order of variants doesn't matter — edits are sorted by cDNA offset internally.

Source code in varcode/mutant_transcript.py

def apply_variants_to_transcript(variants, transcript):
    """Apply a list of variants to a single transcript, yielding one
    :class:`MutantTranscript` that carries all the resulting edits
    (#269). Used for haplotype-aware joint effect prediction — cis
    variants on the same transcript become one combined mutant
    rather than N independent per-variant mutants.

    Edits are applied in cDNA-coordinate order (highest offset
    first, so earlier offsets aren't shifted) to ``transcript``'s
    spliced cDNA. Returns ``None`` when any of the usual
    single-variant preconditions fail (non-coding, incomplete, etc.),
    or when the provided variants conflict — i.e. claim to edit
    overlapping cDNA ranges. The caller is responsible for falling
    back to per-variant effects in that case.

    ``mutant_protein_sequence`` is populated when at least one edit
    lands after the canonical CDS start; the joint cDNA is translated
    from there to the first stop.

    Order of ``variants`` doesn't matter — edits are sorted by cDNA
    offset internally.
    """
    from pyensembl import Transcript

    if not isinstance(transcript, Transcript):
        return None
    if not transcript.is_protein_coding or not transcript.complete:
        return None

    full_sequence = str(transcript.sequence)
    resolved = []
    for variant in variants:
        r = _resolve_variant_edit(variant, transcript, full_sequence)
        if r is None:
            return None
        resolved.append(r)
    # Sort by cDNA start so overlap detection is a simple linear scan
    # and we can apply high-to-low without shifting earlier edits.
    resolved.sort(key=lambda pair: pair[0].cdna_start)
    for i in range(len(resolved) - 1):
        e1 = resolved[i][0]
        e2 = resolved[i + 1][0]
        # Range overlap.
        if e1.cdna_end > e2.cdna_start:
            return None
        # An insertion at the same cdna_start as the next edit is
        # order-dependent (whether the inserted bases survive a
        # following deletion, or where they sit relative to a
        # substitution at the same offset, depends on apply order).
        # In contrast, adjacent edits at e1.cdna_end == e2.cdna_start
        # where e1 is NOT an insertion are unambiguous and allowed:
        # high-to-low application applies e2 first, leaving the slice
        # below e2.cdna_start untouched, so e1 then operates on its
        # original range and the two edits compose to a single
        # deterministic result regardless of input order.
        if e1.cdna_end == e2.cdna_start and e1.cdna_start == e1.cdna_end:
            return None
    # Apply edits from the highest cDNA offset down so earlier offsets
    # aren't affected by upstream edits' length changes.
    mutant_cdna = full_sequence
    for edit, _ in sorted(
            resolved, key=lambda pair: pair[0].cdna_start, reverse=True):
        mutant_cdna = (
            mutant_cdna[:edit.cdna_start]
            + edit.alt_bases
            + mutant_cdna[edit.cdna_end:]
        )
    mutant_protein = None
    cds_start = min(transcript.start_codon_spliced_offsets)
    if any(anchor >= cds_start for _, anchor in resolved):
        mutant_protein = _translate_from_cds(
            mutant_cdna, transcript, [edit for edit, _ in resolved])
    return MutantTranscript(
        reference_transcript=transcript,
        edits=tuple(edit for edit, _ in resolved),
        cdna_sequence=mutant_cdna,
        mutant_protein_sequence=mutant_protein,
        annotator_name="protein_diff",
    )

Observed RNA import

varcode.load_exacto_fusions(structures_path, integrated_path, *, variants_by_id, cds_starts=None, primary_structures_path=None)

Import selected SV-linked RNA models as an RNAEvidence resolver.

PARAMETER	DESCRIPTION
structures_path	Exacto transcript-structures and integrated-variants TSVs (optionally gzip-compressed). Structure index order is transcript order; sequences are already oriented and are NOT reverse-complemented again. TYPE: path or text stream
integrated_path	Exacto transcript-structures and integrated-variants TSVs (optionally gzip-compressed). Structure index order is transcript order; sequences are already oriented and are NOT reverse-complemented again. TYPE: path or text stream
variants_by_id	Exacto DNA call IDs to existing StructuralVariant objects, normally loaded from VCF. Only these IDs are selected from the integration table. The explicit join avoids inventing a DNA breakpoint from an RNA splice. TYPE: mapping
cds_starts	Optional, explicitly chosen zero-based ORF starts in assembled sequence, keyed by (transcript_model_id, tuple(sorted(reference_transcript_ids))). No ORF is chosen automatically. See make_fusion_outcome. TYPE: mapping or None DEFAULT: None
primary_structures_path	Optional native Exacto primary-structures TSV (including peptide_id). Import each peptide separately, validate its codons against the observed RNA, and retain partial-protein status and per-base provenance. Cannot be combined with cds_starts. These are sequence predictions, not protein expression evidence. Models without peptide rows remain untranslated. TYPE: path or text stream or None DEFAULT: None

RETURNS	DESCRIPTION
RNAEvidence	.candidates retains every selected variant/model/reference-ID group. Pass the result to effects(rna_resolver=...). Original structure and integration rows, including RNA and DNA call IDs, remain in evidence.

Notes

Supports linear two-locus models with an annotated sense 5' anchor. A missing or antisense 3' partner remains TranslocationToIntergenic, never a guessed coding fusion. Incomplete sequence, unknown transcript IDs, circular paths and multi-gene (>2) models raise rather than silently losing structure. This does not import all Exacto variant types. Exacto's native primary structures use the standard genetic code; that producer choice is recorded. Read/model completeness and read support are not inferred from row counts.

Source code in varcode/exacto.py

def load_exacto_fusions(structures_path, integrated_path, *, variants_by_id,
                       cds_starts=None, primary_structures_path=None):
    """Import selected SV-linked RNA models as an RNAEvidence resolver.

    Parameters
    ----------
    structures_path, integrated_path : path or text stream
        Exacto transcript-structures and integrated-variants TSVs (optionally
        gzip-compressed). Structure ``index`` order is transcript order;
        sequences are already oriented and are NOT reverse-complemented again.
    variants_by_id : mapping
        Exacto DNA call IDs to existing StructuralVariant objects, normally
        loaded from VCF. Only these IDs are selected from the integration table.
        The explicit join avoids inventing a DNA breakpoint from an RNA splice.
    cds_starts : mapping or None
        Optional, explicitly chosen zero-based ORF starts in assembled sequence,
        keyed by ``(transcript_model_id, tuple(sorted(reference_transcript_ids)))``.
        No ORF is chosen automatically. See ``make_fusion_outcome``.
    primary_structures_path : path or text stream or None
        Optional native Exacto primary-structures TSV (including peptide_id).
        Import each peptide separately, validate its codons against the observed
        RNA, and retain partial-protein status and per-base provenance. Cannot
        be combined with cds_starts. These are sequence predictions, not protein
        expression evidence. Models without peptide rows remain untranslated.

    Returns
    -------
    RNAEvidence
        ``.candidates`` retains every selected variant/model/reference-ID group.
        Pass the result to ``effects(rna_resolver=...)``. Original structure and
        integration rows, including RNA and DNA call IDs, remain in evidence.

    Notes
    -----
    Supports linear two-locus models with an annotated sense 5' anchor. A
    missing or antisense 3' partner remains TranslocationToIntergenic, never a
    guessed coding fusion. Incomplete sequence, unknown transcript IDs, circular
    paths and multi-gene (>2) models raise rather than silently losing structure.
    This does not import all Exacto variant types. Exacto's native primary
    structures use the standard genetic code; that producer choice is recorded.
    Read/model completeness and read support are not inferred from row counts.
    """
    if primary_structures_path is not None and cds_starts is not None:
        raise ValueError("Choose primary structures or explicit cds_starts, not both")
    variants = {str(key): value for key, value in variants_by_id.items()}
    if len(variants) != len(variants_by_id):
        raise ValueError("DNA call IDs collide after conversion to strings")
    if any(not getattr(v, "is_structural", False) for v in variants.values()):
        raise ValueError("variants_by_id must contain only StructuralVariant objects")
    links = defaultdict(list)
    for row in _rows(integrated_path, ["transcript_model_id", "reference_transcript_ids",
                                       "rna_variant_call_id", "dna_variant_call_id"]):
        if row["dna_variant_call_id"] in variants:
            key = (*_model_key(row), row["dna_variant_call_id"])
            if row not in links[key]:
                links[key].append(row)
    selected = {key[:2] for key in links}
    models = defaultdict(list)
    fields = ["transcript_model_id", "reference_transcript_ids", "index",
              "read_start", "read_end", "sequence", "type", "kind", "context",
              "chromosome_1", "position_1", "strand_1", "chromosome_2",
              "position_2", "strand_2", "transcript_id_1", "transcript_id_2"]
    for row in _rows(structures_path, fields):
        key = _model_key(row)
        if key in selected:
            models[key].append(row)
    if selected - set(models):
        raise ValueError("Missing structures for linked transcript models: %s" %
                         sorted(selected - set(models)))
    peptides = defaultdict(lambda: defaultdict(list))
    peptide_models = {}
    if primary_structures_path is not None:
        fields = ["peptide_id", "primary_structure_index", "type", "amino_acid",
                  "amino_acid_index", "codon_index", "nucleotide", "transcript_model_id",
                  "reference_transcript_ids", "transcript_structure_index", "read_start", "read_end"]
        for row in _rows(primary_structures_path, fields):
            key = _model_key(row)
            if key not in selected:
                continue
            peptide = row["peptide_id"]
            if not peptide or peptide_models.setdefault(peptide, key) != key:
                raise ValueError("Missing peptide_id or peptide assigned to multiple RNA models")
            peptides[key][peptide].append(row)
    candidates = []
    for (model_id, refs, dna_id), integration_rows in sorted(links.items()):
        variant = variants[dna_id]
        rows, bases, sequence = _sequence(models[(model_id, refs)])
        _validate_path(bases, variant.genome)
        if any("circular" in row["kind"].lower() or "circular" in row["context"].lower()
               for row in rows):
            raise ValueError("Circular RNA requires a separate model, not a linear fusion")
        first, last = bases[0], bases[-1]
        transcript = _transcript(variant.genome, first["transcript_id_1"])
        if transcript is None or first["strand_1"] != transcript.strand:
            raise ValueError("An annotated sense 5-prime transcript anchor is required")
        partner = _transcript(variant.genome, last["transcript_id_2"])
        transcript_ids = {row[k] for row in bases
                          for k in ["transcript_id_1", "transcript_id_2"] if row[k]}
        genes = {_transcript(variant.genome, tid).gene_id for tid in transcript_ids}
        if len(genes) > 2:
            raise ValueError("More than two genes in an observed fusion path")
        partner_status = "sense" if partner is not None else "unannotated"
        if partner is not None and last["strand_2"] != partner.strand:
            partner_status, partner = "antisense", None
        if partner is not None and partner.gene_id == transcript.gene_id:
            raise ValueError("Selected model does not identify a two-gene fusion")
        evidence = dict(
            dna_variant_call_id=dna_id,
            rna_variant_call_ids=sorted({r["rna_variant_call_id"] for r in integration_rows}),
            reference_transcript_ids=list(refs), partner_status=partner_status,
            exacto_structure=rows, exacto_integration=integration_rows,
            sequence_status="observed_model_completeness_unknown")
        for primary_rows in peptides[(model_id, refs)].values() or [None]:
            candidate = make_fusion_outcome(
                variant, transcript, sequence=sequence, transcript_model_id=model_id,
                partner_transcript=partner, source="exacto",
                cds_start=(cds_starts or {}).get((model_id, refs)),
                extra_evidence=evidence)
            if primary_rows is not None:
                protein, primary_evidence = _primary_protein(primary_rows, rows)
                combined = dict(candidate.evidence, **primary_evidence)
                candidate.effect.mutant_transcript = replace(
                    candidate.effect.mutant_transcript,
                    mutant_protein_sequence=protein, evidence=combined)
                candidate = replace(candidate, evidence=combined)
            candidates.append(candidate)
    return RNAEvidence(candidates)

varcode.make_fusion_outcome(variant, transcript, *, sequence, transcript_model_id, partner_transcript=None, cds_start=None, source='rna', read_count=None, extra_evidence=None)

Import one RNA junction/model using existing structural effect classes.

sequence is the observed 5'-to-3' sequence, never genomic-forward sequence. transcript is its annotated 5' anchor. Supply a 3' partner only when the observed path joins that transcript in sense orientation; otherwise leave it None and retain the locus/orientation in evidence. An absent partner gives TranslocationToIntergenic, not a coding GeneFusion.

No reference exons are appended and no ORF is guessed. An explicit zero-based cds_start requests translation of a complete start-to-stop ORF in the supplied sequence; invalid or incomplete ORFs raise ValueError. Such a protein is sequence-predicted, not evidence of translation.

Source code in varcode/rna_evidence.py

def make_fusion_outcome(
        variant, transcript, *, sequence, transcript_model_id,
        partner_transcript=None, cds_start=None, source="rna", read_count=None,
        extra_evidence=None):
    """Import one RNA junction/model using existing structural effect classes.

    ``sequence`` is the observed 5'-to-3' sequence, never genomic-forward
    sequence. ``transcript`` is its annotated 5' anchor. Supply a 3' partner
    only when the observed path joins that transcript in sense orientation;
    otherwise leave it None and retain the locus/orientation in evidence.
    An absent partner gives TranslocationToIntergenic, not a coding GeneFusion.

    No reference exons are appended and no ORF is guessed. An explicit
    zero-based ``cds_start`` requests translation of a complete start-to-stop
    ORF in the supplied sequence; invalid or incomplete ORFs raise ValueError.
    Such a protein is sequence-predicted, not evidence of translation.
    """
    from .effects.effect_classes import GeneFusion, TranslocationToIntergenic
    from .effects.codon_tables import codon_table_for_transcript, translate_sequence
    from .mutant_transcript import MutantTranscript

    if not getattr(variant, "is_structural", False):
        raise ValueError("RNA fusion import requires a StructuralVariant")
    if not isinstance(sequence, str) or not sequence or set(sequence.upper()) - set("ACGTN"):
        raise ValueError("Observed sequence must be nonempty DNA (ACGTN)")
    if not transcript_model_id:
        raise ValueError("transcript_model_id is required")
    if read_count is not None and (isinstance(read_count, bool)
                                   or not isinstance(read_count, int) or read_count < 0):
        raise ValueError("read_count must be a non-negative integer or None")
    sequence = sequence.upper()
    protein = None
    evidence = dict(extra_evidence or {})
    evidence.update(sequence_status=evidence.get("sequence_status", "unknown"),
                    protein_status="not_determined")
    if cds_start is not None:
        if isinstance(cds_start, bool) or not isinstance(cds_start, int) or cds_start < 0:
            raise ValueError("cds_start must be a non-negative integer")
        coding = sequence[cds_start:]
        table = codon_table_for_transcript(transcript)
        if coding[:3] not in table.start_codons:
            raise ValueError("cds_start does not identify a start codon")
        stop = next((i for i in range(3, len(coding) - 2, 3)
                     if coding[i:i + 3] in table.stop_codons), None)
        if stop is None or "N" in coding[:stop]:
            raise ValueError("A complete, unambiguous start-to-stop ORF is required")
        protein = "M" + translate_sequence(coding[3:stop], codon_table=table)
        evidence.update(cds_start=cds_start, cds_end=cds_start + stop + 3,
                        protein_status="predicted_from_observed_rna")
    evidence.update(transcript_model_id=str(transcript_model_id))
    if read_count is not None:
        evidence["read_count"] = read_count
    model = MutantTranscript.from_sequence(
        sequence, reference_transcript=transcript,
        mutant_protein_sequence=protein, annotator_name=source, evidence=evidence)
    if partner_transcript is None:
        effect = TranslocationToIntergenic(variant, transcript, mutant_transcript=model)
    else:
        if transcript.gene_id == partner_transcript.gene_id:
            raise ValueError("A gene fusion requires two different genes")
        effect = GeneFusion(variant, transcript, partner_transcript, mutant_transcript=model)
    return make_rna_outcome(effect, source=source, extra_evidence=evidence)

varcode.RNAEvidence(candidates=())

Imported candidates, also usable as an RNAEvidenceResolver.

Keeps each observed model separately; read counts do not become likelihoods. A fusion is available on both explicitly identified partner transcripts.

Source code in varcode/rna_evidence.py

def __init__(self, candidates=()):
    self.candidates = tuple(candidates)

RNA evidence

varcode.RNAEvidenceResolver

Bases: Protocol

Source of RNA-observed outcomes for a (variant, transcript) pair.

Implementers return zero or more :class:~varcode.effect_candidates.EffectCandidate objects describing isoforms, fusions, or RNA-level events that were actually observed in reads. An empty sequence means "no evidence for this pair" — the existing DNA-predicted outcomes are left alone.

Returned outcomes should set source to a producer-specific string (the name of the RNA assembler, long-read caller, fusion detector, etc.) and populate evidence with whatever shape that producer natively emits (transcript model IDs, junction read counts, etc.). See :func:make_rna_outcome for a convenience factory that fills the common fields.

observed_outcomes(variant, transcript) -> Sequence[EffectCandidate]

Return RNA-observed outcomes for variant on transcript, or an empty sequence when no evidence is available. Must not raise on unknown (variant, transcript) pairs — return an empty sequence instead.

Source code in varcode/rna_evidence.py

def observed_outcomes(self, variant, transcript) -> Sequence[EffectCandidate]:
    """Return RNA-observed outcomes for ``variant`` on
    ``transcript``, or an empty sequence when no evidence is
    available. Must not raise on unknown ``(variant, transcript)``
    pairs — return an empty sequence instead."""
    ...

varcode.NullRNAEvidenceResolver

No-op resolver that always reports "no evidence".

Useful as a default in pipelines where an RNA resolver is optional and as a baseline in tests. apply_rna_evidence_to_effects is safe to call with this resolver — it's a no-op walk.

varcode.apply_rna_evidence_to_effects(effects: Iterable, resolver) -> Iterable

Attach RNA-observed candidates from resolver to each effect.

Walks effects and, for any effect with a resolvable (variant, transcript), asks resolver.observed_outcomes for RNA-observed candidates.

Splice mechanism sets use RNA evidence as a reconciliation signal: a new set replaces the old one, retaining an audit trail of raw RNA evidence, added candidates, excluded DNA-predicted candidates, and per-current-candidate RNA support. Other multi-outcome effects keep the additive side-channel behavior: observed candidates are stashed on _extra_candidates and exposed through .candidates.

Single-outcome point-variant effects (Missense, FrameShift, etc.) are left untouched even when the resolver has evidence — those classes don't expose a multi-candidate view, and replacing them with a multi-outcome wrapper would break downstream isinstance checks. Producers that need to surface RNA observations on point variants should report them as a separate :class:MultiOutcomeEffect rather than mutating an existing single-outcome one. (The point-variant diff is generally already correct from DNA, so this is rarely an issue in practice.) Single-outcome structural calls (e.g. Intronic) are wrapped in StructuralVariantEffect when observations exist, retaining the DNA classification alongside the imported models.

Safe to call on a mixed collection where only some variants have RNA evidence; no-op when resolver is None or doesn't implement the protocol.

Returns effects for chaining convenience.

Source code in varcode/rna_evidence.py

def apply_rna_evidence_to_effects(effects: Iterable, resolver) -> Iterable:
    """Attach RNA-observed candidates from ``resolver`` to each effect.

    Walks ``effects`` and, for any effect with a resolvable
    ``(variant, transcript)``, asks ``resolver.observed_outcomes``
    for RNA-observed candidates.

    Splice mechanism sets use RNA evidence as a reconciliation signal:
    a new set replaces the old one, retaining an audit trail of raw RNA
    evidence, added candidates, excluded DNA-predicted candidates, and
    per-current-candidate RNA support. Other multi-outcome effects keep
    the additive side-channel behavior: observed candidates are stashed
    on ``_extra_candidates`` and exposed through ``.candidates``.

    Single-outcome point-variant effects (Missense, FrameShift, etc.) are left
    untouched even when the resolver has evidence — those classes
    don't expose a multi-candidate view, and replacing them with a
    multi-outcome wrapper would break downstream ``isinstance`` checks.
    Producers that need to surface RNA observations on point variants
    should report them as a separate :class:`MultiOutcomeEffect` rather
    than mutating an existing single-outcome one. (The point-variant
    diff is generally already correct from DNA, so this is rarely an
    issue in practice.) Single-outcome structural calls (e.g. Intronic) are
    wrapped in StructuralVariantEffect when observations exist, retaining the
    DNA classification alongside the imported models.

    Safe to call on a mixed collection where only some variants have
    RNA evidence; no-op when ``resolver`` is None or doesn't implement
    the protocol.

    Returns ``effects`` for chaining convenience.
    """
    if resolver is None:
        return effects
    if not hasattr(resolver, "observed_outcomes"):
        return effects

    # Lazy import to avoid an import cycle (effect_classes imports
    # from varcode.effect_candidates, which sits below us).
    from .effects.effect_classes import MultiOutcomeEffect, StructuralVariantEffect

    replacements = []
    changed = False
    for effect in effects:
        is_structural = getattr(getattr(effect, "variant", None), "is_structural", False)
        if not isinstance(effect, MultiOutcomeEffect) and not is_structural:
            replacements.append(effect)
            continue
        variant = getattr(effect, "variant", None)
        transcript = getattr(effect, "transcript", None)
        if variant is None or transcript is None:
            replacements.append(effect)
            continue
        observed = resolver.observed_outcomes(variant, transcript)
        if not observed:
            replacements.append(effect)
            continue
        # Coerce to tuple eagerly so we don't keep a generator that
        # would silently exhaust on a second candidates() read.
        observed_tuple = tuple(observed)
        if not observed_tuple:
            replacements.append(effect)
            continue
        if not isinstance(effect, MultiOutcomeEffect):
            # An intronic/noncoding DNA call must not hide an observed SV model.
            effect = StructuralVariantEffect(
                variant, transcript, primary_effects=(effect,))
            changed = True
        if hasattr(effect, "with_rna_evidence"):
            refined = effect.with_rna_evidence(observed_tuple)
            replacements.append(refined)
            changed = changed or refined is not effect
        else:
            existing = getattr(effect, "_extra_candidates", ())
            effect._extra_candidates = tuple(existing) + observed_tuple
            replacements.append(effect)
    if changed:
        return _replace_effects(effects, replacements)
    return effects

varcode.make_rna_outcome(effect, *, source: str = 'rna', transcript_model_id: Optional[str] = None, read_count: Optional[int] = None, extra_evidence: Optional[Mapping[str, Any]] = None) -> EffectCandidate

Construct an :class:~varcode.effect_candidates.EffectCandidate carrying RNA-derived provenance.

Convenience factory for the common fields a reads-based or long-read assembly tool wants on each observed outcome — keeps consumers from hand-rolling the evidence dict shape and lets downstream code rely on a small set of well-known keys.

PARAMETER	DESCRIPTION
effect	The effect this RNA-observed outcome represents. TYPE: MutationEffect
source	Producer name; defaults to "rna". Set to a tool-specific string (RNA assembler, long-read caller, etc.) for downstream filtering. Opaque to varcode. TYPE: str DEFAULT: 'rna'
transcript_model_id	Stable ID of the observed transcript model from the producer. Stored under evidence["transcript_model_id"]. TYPE: str or None DEFAULT: None
read_count	Supporting read count. Stored under evidence["read_count"]. TYPE: int or None DEFAULT: None
extra_evidence	Producer-specific extra fields, merged into the evidence dict on top of the well-known keys above. Allows tool-native fields (e.g. "tpm", "junction_id") to ride along without forcing a schema here. TYPE: Mapping or None DEFAULT: None

Source code in varcode/rna_evidence.py

def make_rna_outcome(
        effect,
        *,
        source: str = "rna",
        transcript_model_id: Optional[str] = None,
        read_count: Optional[int] = None,
        extra_evidence: Optional[Mapping[str, Any]] = None) -> EffectCandidate:
    """Construct an :class:`~varcode.effect_candidates.EffectCandidate`
    carrying RNA-derived provenance.

    Convenience factory for the common fields a reads-based or
    long-read assembly tool wants on each observed outcome — keeps
    consumers from hand-rolling the ``evidence`` dict shape and lets
    downstream code rely on a small set of well-known keys.

    Parameters
    ----------
    effect : MutationEffect
        The effect this RNA-observed outcome represents.
    source : str
        Producer name; defaults to ``"rna"``. Set to a tool-specific
        string (RNA assembler, long-read caller, etc.) for downstream
        filtering. Opaque to varcode.
    transcript_model_id : str or None
        Stable ID of the observed transcript model from the producer.
        Stored under ``evidence["transcript_model_id"]``.
    read_count : int or None
        Supporting read count. Stored under ``evidence["read_count"]``.
    extra_evidence : Mapping or None
        Producer-specific extra fields, merged into the evidence dict
        on top of the well-known keys above. Allows tool-native fields
        (e.g. ``"tpm"``, ``"junction_id"``) to ride along without
        forcing a schema here.
    """
    evidence: dict = {}
    if transcript_model_id is not None:
        evidence["transcript_model_id"] = transcript_model_id
    if read_count is not None:
        evidence["read_count"] = read_count
    if extra_evidence:
        evidence.update(extra_evidence)
    return EffectCandidate(
        effect=effect,
        source=source,
        evidence=evidence,
    )