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. Matching sense exonic bases become reference-transcript segments; other bases remain observed sequence. No reference bases are appended.

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.
        Matching sense exonic bases become reference-transcript segments;
        other bases remain observed sequence. No reference bases are appended.

    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)
        segments = _reference_segments(bases, variant.genome, sequence)
        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)
            candidate.effect.mutant_transcript = replace(
                candidate.effect.mutant_transcript, reference_segments=segments)
            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)

Reference-completion hypotheses

varcode.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.

PARAMETER	DESCRIPTION
candidate	Observed GeneFusion or TranslocationToIntergenic with materialized cDNA and reference segments. The candidate is never modified. TYPE: EffectCandidate
five_prime_transcripts	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. TYPE: iterable or None DEFAULT: None
three_prime_transcripts	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. TYPE: iterable or None DEFAULT: None
ends	Ends to complete: five_prime, three_prime, or both (default). TYPE: iterable of str DEFAULT: ('five_prime', 'three_prime')

RETURNS

DESCRIPTION

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	DESCRIPTION
ValueError	For unsupported effects, malformed layouts, invalid end evidence, or an already reference-completed input.

Source code in varcode/reference_completion.py

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)

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,
    )