Design: one pipeline for splice, germline-phase, haplotype and SV effects (varcode 8)

Summary

varcode has 48 effect classes and five separate ways to express uncertainty. Splice outcomes, unknown germline phase, known-cis haplotypes and structural-variant alternatives each wrap the output of the step before them. As #421 shows, the wrappers don't compose: germline phase ignores splicing, germline-aware annotation ignores the SV annotator, and haplotypes containing a splice variant or an SV are silently dropped.

This proposes replacing all of that with one pipeline used for every variant type:

hypotheses(variants, transcript, germline, phase)   # every combination of phase × splice mechanism × SV resolution
  → realize(hypothesis)                              # genomic layout with edits applied, then spliced into mRNA
  → classify(mutant vs. patient baseline)            # the existing protein-diff classifier
  → merge identical local proteins, order by weight  # effect = top candidate; alternatives in .candidates

What changes for users:

• One Variant class. SVs are variants with a symbolic alt and an end or mate. No placeholder N>A alleles (#417).
• 18 effect classes (down from 48), each naming a consequence or a location: Deletion, FrameShift, StartLoss, Intronic, …
  • LargeDeletion, ExonSkipping, SpliceOutcomeSet, PhaseCandidateSet, HaplotypeEffect and friends are gone (#420).
  • How a consequence arises (mechanism: exon_skip, sv_type: DEL, phase: cis) is evidence on the candidate, not a class.
• The top-level effect is an ordinary effect. type(effect), aa_ref and mutant_protein_sequence mean what they say. Alternatives are in effect.candidates, ordered by weight.
• Genomic sequence is optional and raises resolution. Without a FASTA, exon-level outcomes are fully resolved, and intron retention and cryptic sites are Unresolved. With one, everything is realized from sequence, including splice-signal detection on the mutated haplotype. In all 33 examples below, the top effect is the same with or without genomic sequence; the extra tier fills in the alternatives.

Worked examples follow in the comments below. Every number there was computed on CFTR ENST00000003084 (Ensembl 81, GRCh38 FASTA), using a throwaway prototype of this design that reuses varcode's classifier and splice-motif scorer. The prototype is attached at the end so reviewers can rerun or change it.

Related: #412, #416, #417, #418, #420, #421, #422.

The model

1. Variants

One class, Variant(contig, start, ref, alt, end=None, mate=None, info=None):

• alt can be symbolic: a nucleotide string, <DEL> / <DUP> / <INV> / <INS> / <CNV>, or a breakend.
• kind property: snv, mnv, indel, DEL, DUP, INV, INS, CNV or BND.
• Allele predicates: is_snv, is_indel and is_transition are only true for nucleotide alleles.
• SV-only fields: end, mate_contig / mate_start / mate_orientation, ci_start / ci_end and alt_assembly are None on point variants.

Breakend stays as an internal value type. Loaders and effects() dispatch on kind, never on the class.

2. Hypotheses

A hypothesis is one full assignment of every open question for a (variants, transcript) pair. Its weight is the product of the weights of its choices. There are three axes.

Phase. One axis per germline variant in the window:

• Het, phase unknown: cis 0.5 / trans 0.5. This is the no-information prior; varcode already enumerates the same two hypotheses in germline.enumerate_phase_hypotheses.
• Phase known from a resolver (VCF PS, reads, assembly): a single choice, weight 1.
• Homozygous: always in the baseline, weight 1.
• Somatic–somatic phase: never enumerated. Joint effects are built only when a resolver says cis; otherwise variants are annotated one at a time (section 8).

Splicing. Axes are splice sites, keyed by (donor|acceptor, exon, occurrence), so the same site has the same key in the baseline and the mutant even when an SV removes exons between them. Each site's status is computed on the realized baseline and on the realized mutant (section 3). The comparison decides whether the site becomes an axis:

baseline status vs. mutant status	caused by	axis?	baseline splices as
same, strong	germline	yes, shared: baseline and mutant take the same choice	the shared choice (C2)
same, weak	germline	no, treated as intact	canonical (E2)
differ	the somatic variant	yes, mutant only	canonical, or its own most likely choice if the baseline is itself strongly disrupted (C1, C3)

Why each rule:

• Shared choice for strong germline disruptions. Both molecules carry the broken site, so they splice the same way. Classifying the somatic change under the same choice isolates what the somatic variant does.
• Germline weak disruptions are ignored. SNPs in the last three exonic bases or at donor +3..+6 are common, and most don't change splicing. Enumerating them without a splice scorer would multiply every result by four for no information. A tier-2 scorer (section 7) can promote them to axes.
• The baseline keeps its own outcome when the somatic variant causes the disruption. The baseline is a reference frame, what the patient's haplotype makes without the somatic variant. An earlier draft shared the choice here, which classified a somatic donor +1 variant as "no change" whenever a germline SNP also weakened that donor.

Options and weights by the mutant's disruption strength:

disruption	detected when	options (weight)
strong	exon end/start removed, donor +1/+2 not GT, or acceptor −2/−1 not AG	exon skip 0.4, intron retention 0.3, cryptic site 0.2, leaky normal 0.1
weak, intronic	donor +3..+6 or acceptor −3 changed	normal 0.4, exon skip 0.3, intron retention 0.2, cryptic 0.1
weak, exonic	last 3 exonic bases before a donor, or first exonic base after an acceptor, changed	normal 0.4, exon skip 0.3, cryptic 0.2, intron retention 0.1

• Where the order comes from: these are varcode's existing _MECHANISM_ORDER_SPLICE_DONOR / _ACCEPTOR / _EXONIC_SPLICE_SITE / _INTRONIC_SPLICE_SITE_* tuples (splice_outcomes.py:479-520), unchanged.
• The numbers are placeholders. The 0.4/0.3/0.2/0.1 values are only a rank-to-weight mapping. Calibrating them, for example against RNA-seq studies of splice-site variant outcomes, is an open question below.
• The windows are today's position windows (SpliceDonor, SpliceAcceptor, IntronicSpliceSite, ExonicSpliceSite).
• Impossible options are dropped and the remaining weights renormalized. Examples: skipping the first or last exon, or no cryptic site with a positive motif score within ±50 bp (cryptic_exons.score_donor / score_acceptor).

Structural variants. An SV's junctions are part of the realized layout, not a separate axis. Its alternatives come out of the splice axis automatically: a junction that removes a donor is a disrupted junction like any other (examples D4, D5, D8). Breakpoint-uncertainty hypotheses (ci_start / ci_end) could be added later as another axis.

Which germline variants are enumerated. A variant is included if it falls within:

• the protein-changing region of any reference-relative candidate, ±90 nt,
• the splice windows of any disrupted junction, or
• the span of an SV.

Everything else is ignored, which keeps enumeration small; the cap is 64 hypotheses (varcode's max_hypotheses). Merging (section 5) makes the result insensitive to germline variants that don't change the protein near the somatic change (B5).

3. Realization

Each hypothesis is realized twice: once as the patient baseline (germline edits for that phase assignment only) and once as the mutant (baseline + somatic edits).

1. Genomic layout. Start from the reference bases of the transcript's span (± flank) and apply edits in genomic coordinates:

   • SNVs and indels replace bases,
   • DEL removes a span,
   • DUP inserts a tagged copy,
   • INV inserts the reverse complement tagged antisense,
   • BND concatenates two layouts.

   Every base keeps its origin (reference position, inserted, duplicated or inverted).

2. Exon runs. Walk the layout and collect maximal runs of sense-strand bases that fall in the transcript's exons. An antisense (inverted) exon is not an exon run (D7).

3. Junction status. For each consecutive pair of runs, check the donor at the end of the first and the acceptor at the start of the second (table above). With genomic sequence the check reads the realized bases, so a germline or SV that creates or destroys GT/AG is seen. Without it, the check is positional.

4. Splice. Apply the hypothesis's choice at each disrupted junction:

   • skip the run,
   • include the realized intervening sequence,
   • move the boundary to the best cryptic site in the realized sequence,
   • or join as normal.

   Concatenate to get the mRNA.

5. Translate from the reference start codon's base (if it survives in the mRNA) to the first stop.

Which choices the baseline takes, shared or its own, follows the table in section 2.

4. Classification

• Protein changed: classify_from_protein_diff(mutant_protein, baseline_protein, cds_length_delta), the existing classifier, unchanged. It yields one of Substitution, ComplexSubstitution, Insertion, Deletion, FrameShift, PrematureStop, StopLoss, StartLoss or Silent.

• Protein identical: location from the first mRNA difference:

  • mRNA identical → Intronic
  • first difference before the CDS → FivePrimeUTR
  • inside the CDS → Silent
  • after it → ThreePrimeUTR

  An exonic variant whose exon is spliced out on that haplotype therefore comes out Intronic. On that molecule its bases leave with the intron (C2, C3). Whether that deserves its own name is an open question.

• Start codon deleted from the genome: StartLoss for every hypothesis, even without genomic sequence, because no splicing choice can restore it (D3).

• mRNA contains exons from two genes: GeneFusion(in_frame=...), carrying mutant_protein_sequence when translatable.

• Protein needs sequence the caller didn't provide: Unresolved(mechanism=...).

• Unchanged from today: Intergenic, NoncodingTranscript, IncompleteTranscript and Failure.

Every effect also carries splice_signal (e.g. "donor+1"), set when the variant touches a splice window, independent of its class.

5. Merging

Hypotheses with the same class and the same local mutant protein merge into one candidate whose weight is the sum. Local protein means the residues within 30 aa of the change, enough to cover any MHC-I or MHC-II peptide spanning it.

• Why local, not full protein: a germline SNP 40 codons away splits the full protein but no peptide (B5).
• Why not class only: a germline SNP 4 codons away changes the peptides (B4), and a cis germline next to a new exon junction changes the junction peptides (E2).
• What's kept: each candidate lists all of its hypotheses, so no provenance is lost.
• Unresolved candidates never merge across phase hypotheses, since nothing shows their proteins are equal. Without this, an unresolved intron retention that merged across phases could outrank resolved candidates that didn't.

6. Ordering and the top effect

Within one transcript, candidates are ordered by:

1. weight, descending;
2. effect_priority of the class, descending (the existing priority list, minus removed classes; Unresolved sits just above Intronic);
3. fewer cis germline variants (the reference-like reading first);
4. enumeration order.

effect is candidates[0]. It is an ordinary effect of that class, with .candidates, .weight and .highest_priority_candidate (the most severe candidate, for "flag if anything could be a frameshift" filters).

Across transcripts, top_priority_effect is unchanged: it applies the existing filters and tie-breaks (priority, protein-coding gene/transcript, complete transcript, CDS length, transcript length, exon count, name) to the per-transcript top effects.

Every example below states which ordering rule decided the top effect.

7. Resolution tiers

Genomic sequence is optional and each extra input sharpens the same candidates rather than adding classes:

tier	inputs	what is resolved
0	pyensembl only	all exon-level outcomes: point changes, exon skips, whole-exon SV deletions and duplications, junctions between existing exons. Intron retention and cryptic sites are Unresolved (weight kept, protein unknown). Splice-signal detection is positional.
1	+ genome FASTA	intron retention and cryptic donors/acceptors realized from sequence. Impossible cryptic options are dropped. Splice signals are checked on the realized haplotype, so SV junctions and germline variants that create or destroy GT/AG are caught. SV junctions into introns or intergenic sequence are translated. Intronic reference alleles are validated.
2	+ a splice scorer (SpliceAI, Pangolin, …)	replaces the placeholder mechanism weights per junction. No structural change.
3	+ RNA junctions or assembled contigs / alt_assembly	each observed isoform is itself a realized mRNA. It reweights matching candidates by junction signature and adds observed-only candidates (section 9).

8. Multiple somatic variants

• Per-variant effects (unchanged contract): each somatic variant is classified alone against the patient baseline.
• Joint effect for a cis group: when a phase resolver puts somatic variants in cis, one extra effect is built from a hypothesis containing all their edits (G1–G3). It is an ordinary effect class with variants=(…), replacing HaplotypeEffect, and the per-variant effects remain.

9. RNA evidence

A candidate's realized mRNA has a junction signature, the list of non-reference donor^acceptor joins (shown in every example's evidence column).

• Match: an RNA-observed junction or isoform matches the candidates whose signatures contain it.
• Unobserved candidates stay in candidates with rna_support=0 rather than being deleted.
• Observed-only isoforms are realized and classified like any other hypothesis, and added with source="rna".

This replaces the mechanism-class matching in SpliceOutcomeSet.with_rna_evidence, and works the same for SV and haplotype candidates.

Classes after the change

kind	classes
location, no protein change	Intergenic, NoncodingTranscript, IncompleteTranscript, FivePrimeUTR, ThreePrimeUTR, Intronic, Silent
protein change	Substitution, ComplexSubstitution, Insertion, Deletion, FrameShift, PrematureStop, StopLoss, StartLoss, GeneFusion
other	Unresolved, Failure
data	Variant, MutantTranscript, EffectCandidate, EffectCollection, VariantCollection

Removed:

• StructuralVariant (folded into Variant)
• Intragenic (folded into Intergenic)
• the abstract intermediates TranscriptMutationEffect, Exonic, CodingMutation, NonsilentCodingMutation, KnownAminoAcidChange, SpliceSite
• the splice-location classes IntronicSpliceSite, SpliceDonor, SpliceAcceptor, ExonicSpliceSite, which become effect.splice_signal
• the splice-mechanism classes SpliceMechanismEffect, NormalSplicing, ExonSkipping, IntronRetention, CrypticSpliceSiteEffect, CrypticDonor, CrypticAcceptor, CrypticExonCandidate
• the SV classes StructuralVariantEffect, LargeDeletion, LargeDuplication, Inversion, TranslocationToIntergenic
• the wrappers ExonLoss, MultiOutcomeEffect, SpliceOutcomeSet, PhaseCandidateSet, HaplotypeEffect
• AlternateStartCodon (becomes Silent with ref_codon / alt_codon)
• FrameShiftTruncation (becomes FrameShift with an empty shifted sequence)

Downstream code uses very little of this:

• topiary maps class names to labels (predictor.py:315-325) and checks NonsilentCodingMutation (filters.py:52).
• vaxrank checks MultiOutcomeEffect (varcode_effects.py:18) and type(effect).__name__ == "FrameShift" (external_input.py:116).
• isovar uses ExonicSpliceSite and CodingMutation (effect_prediction.py:13).

Each needs a small PR alongside varcode 8.0.

Rollout

1. One Variant class. Also fixes #417.
2. Genomic-layout MutantTranscript, with realization of edits and splicing. Tier 0 first; the existing SV and splice builders become thin callers.
3. Hypotheses, merging, ordering. effects() switches to the new engine behind a flag. The current parity and divergence test suites plus the examples below become the acceptance tests.
4. Remove the old classes. Bump to 8.0; downstream PRs land the same day.
5. Tier 1 / 2 / 3 inputs, as follow-ups.

Open questions for reviewers

1. Weights. The mechanism order is varcode's, but the 0.4 / 0.3 / 0.2 / 0.1 values are a placeholder rank mapping. What should calibrate them? And should a tier-0 run report weights at all, or only an order?
2. Top effect rule. Highest weight first means a canonical donor variant's top effect is exon skipping (Deletion), not the "worst" reading. vaxrank already has a highest-priority selection mode; should the default be weight or priority?
3. "Intronic" for a skipped exon. An exonic variant that is spliced out on its haplotype (C2, C3) classifies as Intronic. Is that the right name, or does it need ExcludedExon or a flag?
4. Intron retention labels. The classifier calls a 3,161 bp retained intron a FrameShift, because the CDS length change isn't a multiple of 3, even though translation reads a few codons into the intron and stops (A4). Should read-through-to-stop get its own label, or is FrameShift / PrematureStop by length enough?
5. Local merge window. ±30 aa is sized for MHC-I and MHC-II peptides. Should it be configurable, and should merging ever look at the full protein?
6. Weak germline splice changes. Ignored without a scorer. Is that too aggressive for germline variants at +3..+6, where some are known pathogenic?
7. Somatic–somatic phase. Unknown somatic phase never enumerates a joint effect (G1–G4 assume cis from a resolver). Should unknown somatic phase enumerate cis/trans like germline does?
8. Hypothesis cap. 64 by default. When it's hit, should varcode drop low-weight branches, or fall back to per-axis marginals?
9. Baseline splicing uncertainty. When the baseline is strongly disrupted but the somatic variant causes a different disruption, the baseline takes its single most likely outcome instead of enumerating (C1, C3). Enumerating it too would make the baseline a distribution; is that ever worth it?

How the examples were computed

Each example below was run three ways:

• Today: varcode main at 1f65898 on the same inputs (Variant.effects, VariantCollection.effects with germline= / phase_resolver=, or StructuralVariant), keeping the top-level effects on ENST00000003084.
• New model with genomic sequence: a prototype of this design over GRCh38 (hg38.fa).
• New model without genomic sequence: the same prototype with intron retention and cryptic sites left Unresolved.

What the prototype does:

1. Builds the genomic layout.
2. Applies the edits.
3. Finds exon runs and splice-site status.
4. Enumerates the axes in the section 2 table.
5. Splices and translates from the reference ATG.
6. Classifies with varcode's own classify_from_protein_diff.
7. Finds cryptic sites with varcode's _best_cryptic_site / score_donor / score_acceptor, applied to the realized mutant sequence.

Checks:

• Reference data: exon-concatenated hg38.fa reproduces pyensembl's 6,132 nt CFTR cDNA and 1,480 aa protein exactly.
• Agreement with varcode where both should agree:
  • p.L159M, p.F143fs, p.R153=
  • exon 4 skip p.E92_K163del
  • intron 4 retention (184 aa)
  • exon 5 skip p.T164_E193del
  • the B1 / B2 germline pairs
• Where they disagree, the difference is either explained in the example or already filed: #422 (varcode's cryptic donor scan reads the reference, off by one) and #412 / #421.

Every reference allele in the inputs was read from the FASTA. Coordinates are 1-based GRCh38; SV spans are the affected bases (VCF POS+1..END).

The prototype only handles forward-strand, single-contig transcripts. That's enough for these examples, but it isn't an implementation plan. Its source is in the last comment.

Example index

Every example is worked in full in the comments below: inputs, today's output, every hypothesis with its realized mRNA and evidence, the merge, the ordered candidates, and why the top effect wins.

#	case	today (top-level)	new top effect	without genomic sequence
A1	Missense SNV, nothing else in play	Substitution	Substitution p.L159M	same
A2	1 bp deletion in an exon	FrameShift	FrameShift p.F143fs, protein 151 aa	same
A3	Synonymous SNV mid-exon	Silent	Silent no protein change	same
A4	Canonical donor +1 SNV	SpliceOutcomeSet	Deletion p.E92_K163del (72 aa)	same
A5	Donor +5 SNV (splice region)	SpliceOutcomeSet	Intronic no protein change	same
A6	Synonymous SNV on the last base of an exon	SpliceOutcomeSet	Silent no protein change	same
A7	Canonical acceptor -2 SNV	SpliceOutcomeSet	Deletion p.T164_E193del (30 aa)	same
A8	1 bp deletion in the last three bases of an exon	SpliceOutcomeSet	FrameShift p.K163fs, protein 164 aa	same
B1	Missense SNV + het germline SNV in the same codon, phase unknown	PhaseCandidateSet	Substitution p.L159M	same
B2	Same as B1, phase known (cis)	Substitution	Substitution p.S159T	same
B3	Same as B1, germline homozygous	—	Substitution p.S159T	same
B4	Missense SNV + het germline missense 4 codons away	Substitution	Substitution p.L159M	same
B5	Missense SNV + het germline missense far away	Substitution	Substitution p.L159M	same
B6	Missense SNV + het germline stop-gain upstream	Substitution	Substitution p.L159M	same
C1	Donor +1 SNV + het germline in the exon that gets skipped	PhaseCandidateSet	Deletion p.E92_K163del (72 aa)	same
C2	Missense SNV + het germline that breaks the donor of the same exon	Substitution	Substitution p.L159M	same
C3	Het germline that breaks exon 5's acceptor + somatic deletion of exon 5	PhaseCandidateSet	Deletion p.T164_E193del (30 aa)	same
D1	Deletion of exon 5 (intronic breakpoints)	LargeDeletion	Deletion p.T164_E193del (30 aa)	same
D2	Deletion of exons 5–6	LargeDeletion	FrameShift p.T164fs, protein 171 aa	same
D3	Deletion from exon 1 (including the ATG) into intron 3	LargeDeletion	StartLoss p.M1?	same
D4	Deletion that removes exon 5's donor (anchored on the last exon base)	LargeDeletion	Deletion p.T164_E193del (30 aa)	same
D5	Deletion starting at exon 5 donor +3	Intronic	Intronic no protein change	same
D6	Tandem duplication of exon 5	LargeDuplication	Insertion p.T164_E193dup (30 aa)	same
D7	Inversion of exon 5 (intronic breakpoints)	Inversion	Deletion p.T164_E193del (30 aa)	same
D8	Inversion with one breakpoint inside exon 5	Inversion	Deletion p.T164_E193del (30 aa)	same
D9	Tandem duplication with one breakpoint inside exon 5	LargeDuplication	Intronic no protein change	same
E1	Exon 5 deletion + het germline inside exon 5	PhaseCandidateSet	Deletion p.T164_E193del (30 aa)	same
E2	Exon 5 deletion + het germline on the first base of exon 6	LargeDeletion	Deletion p.T164_E193del (30 aa)	same
F1	Donor-removing deletion (D4) + het germline inside exon 5	LargeDeletion	Deletion p.T164_E193del (30 aa)	same
G1	Two somatic SNVs in the same codon, known cis	HaplotypeEffect; Substitution; Substitution	Substitution p.L159T	same
G2	Two somatic 1 bp indels 9 bp apart, known cis (compensating frameshifts)	HaplotypeEffect; FrameShift; FrameShift	ComplexSubstitution p.F143_145delinsLAY	same
G3	Donor +1 SNV + missense SNV in the same exon, known cis	Substitution; SpliceOutcomeSet	Deletion p.E92_K163del (72 aa)	same
G4	Exon 5 deletion + missense SNV inside exon 5, known cis	LargeDeletion; Substitution	Deletion p.T164_E193del (30 aa)	same

Start by reading splice_outcomes.py:479-520, germline.enumerate_phase_hypotheses, cryptic_exons.score_donor and score_acceptor, and classify_from_protein_diff. Review the attached prototype and the related issues before attempting the proposed unified pipeline. Done requires resolving the open design questions and validating the behavior against the stated examples.