varcode.transforms.combine_cis_snvs — combine in-cis adjacent SNVs into MNVs

[iskandr]

Follow-up from #364 / #367 (the varcode.transforms module shape). One of four planned transforms enumerated in docs/transforms.md and the module docstring.

Scope

Reduce two or more adjacent SNVs that share a phase set (or are otherwise determined cis) into a single MNV. Removes the need for HaplotypeEffect/PhaseAmbiguousEffect wrapping in the common case where two SNVs sit in the same codon on the same haplotype and can be classified as one combined codon substitution.

Signature:

def combine_cis_snvs(vc, phase_resolver) -> VariantCollection:
    """Combine adjacent in-codon SNVs sharing a phase set into MNVs.

    Cardinality: reduces.
    """

Contract (per transforms module conventions)

• Cardinality: reduces.
• Provenance: combined MNV carries source_variants=(snv1, snv2, ...).
• Metadata behavior: GT must agree across all source SNVs (raises on mismatch — same rule as pair_breakends); other FORMAT fields taken from the lex-earlier source.

Pairing rule

Two SNVs are combined when all of:

1. Both share a phase set (the phase_resolver answers in_cis(a, b) is True), OR both are homozygous-alt for the same sample at distinct positions within a single codon.
2. They sit within a transcript codon window (3 bp at the same codon position).
3. They're on the same contig.

Tests

• Two SNVs at codon positions 1+2 of the same codon, in cis -> single MNV; effect prediction emits one Substitution instead of two adjacent ones.
• Two SNVs in cis but spanning a codon boundary -> not combined (out of codon window).
• Two SNVs in trans -> not combined.
• Phase unknown -> not combined; existing PhaseAmbiguousEffect path still handles them.
• Three SNVs in cis within one codon -> one MNV with source_variants=(a, b, c).

Composition

Should compose cleanly after pair_breakends (separate scopes; both reduce; both idempotent).

See also

• Pair BND breakend rows into a single StructuralVariant via MATEID #364 — BND mate pairing (the first transform).
• varcode.transforms + pair_breakends (#364, planning + scaffold) #367 — varcode.transforms module shape PR.
• Phasing: cis/trans-aware effect prediction for nearby variants #269 — phase resolver infrastructure.

[iskandr]

Implementation plan

Drafting this here for visibility before any code lands. Not starting yet — flagged as complicated.

Signature

def combine_cis_snvs(
    vc: VariantCollection,
    *,
    phase_resolver=None,
) -> VariantCollection

phase_resolver is keyword-only so the call site reads combine_cis_snvs(vc, phase_resolver=resolver). When None, only the homozygous-alt path can establish cis (matches the pattern in apply_phase_resolver_to_effects).

Core idea

Two SNVs are combined when:

1. Same codon on at least one shared transcript — transcript.spliced_offset(pos_a) // 3 == transcript.spliced_offset(pos_b) // 3 (strand-handled by spliced_offset).
2. In cis — homozygous-alt for both, OR a shared PS (phase set) tag in a sample's FORMAT field, OR phase_resolver.in_cis(a, b) is True.

If both hold, replace the SNVs with a single combined variant covering the codon span. Strictly an MNV if positions are contiguous; if there's a gap of 1 bp between the SNVs (positions 1+3 of a codon), we fill the middle base from the reference (via reference_range, the helper already added in #373).

Algorithm

1. Bucket all SNVs by (transcript_id, codon_index) using spliced_offset.
   Variants not in any coding region pass through unchanged.

2. For each bucket with >= 2 SNVs:
   - Establish cis across the group: every pair must satisfy _all_cis().
   - If cis: build one combined MNV; mark sources as consumed.
   - If not cis: leave all members unchanged.

3. Output = (untouched variants) + (combined MNVs), sorted.

A variant might land in multiple transcript buckets (overlapping transcripts). The dedup rule: a variant gets combined at most once. If two transcripts each say "this SNV combines with a different partner," we pick the first transcript by lex order of transcript_id and log a debug. The other transcript's view degrades to pass-through SNVs — annotator already handles that.

Helpers

• _codon_buckets(vc) -> dict[(transcript_id, codon_idx), list[Variant]] — single pass over coding SNVs.
• _all_cis(snvs, resolver) -> bool — pairwise check across the bucket; short-circuits.
• _build_combined_mnv(snvs, reference) -> Variant — sorts by position, fills gap bases via reference_range, constructs MNV with merged metadata.

Metadata

Field	Rule
source_variants	tuple of the merged SNVs, sorted by position
id	lex-earliest source ID
qual	min(qual) across sources where present
filter	union of FILTER tokens; PASS drops if any non-PASS
info	from lex-earliest source; adds original_starts: tuple[int, ...] and combined_from_cis_snvs: True
sample_info	per-sample GT must match across all sources (raises with all source IDs + GT values on disagreement, same rule as pair_breakends); other FORMAT fields from lex-earliest

Test matrix (~15 cases)

Pass-throughs

• No SNVs in VC -> unchanged
• Single SNV -> unchanged
• Two SNVs different codons -> unchanged
• Two SNVs same codon, trans (resolver says not in_cis) -> unchanged
• Two SNVs same codon, phase unknown (no homozygous, no PS, no resolver) -> unchanged
• Two SNVs in an intron -> unchanged (no codon bucket)

Real combines

• Codon positions 1+2 in cis via shared PS -> single MNV
• Codon positions 1+3 in cis (middle base from reference) -> MNV with filled gap
• Three SNVs same codon in cis -> one MNV, source_variants has 3 entries
• Both homozygous-alt -> combined (no resolver needed)
• Resolver-determined cis -> combined

Edge cases

• GT disagreement -> raises with both source IDs and GTs
• Combined MNV → effects() emits one Substitution instead of two
• Idempotence: running twice produces same VC (second pass finds no candidate buckets)
• Composition with pair_breakends and left_align_indels: order-independent (SVs aren't SNVs; left-aligned SNVs are unchanged)
• Multi-transcript overlap: same SNV in two transcript buckets with different partners → combines with the lex-earlier transcript's partner; other left as pass-through SNV

Tricky bits

1. Strand handling: codon index from spliced_offset is strand-aware, so a - strand transcript's codon at genomic positions {100, 99, 98} resolves to the same codon_index. The combined MNV's ref/alt are still on the + strand — that's varcode's normalized form.

2. Reference for gap fill: only needed for non-contiguous SNVs within a codon (positions 1+3). Reads via reference_range(genome, contig, start, end) — works under both varcode.Genome(fasta=...) and pyensembl-only setups (cDNA fallback).

3. Codon boundaries: two cis SNVs that happen to sit at the boundary of two codons (positions 3+1 of consecutive codons) do not combine — they're in different codon buckets. This is intentional: the effect is two adjacent substitutions, not a multi-codon event.

4. Same-position SNVs: two A→C and A→G at the same locus — multi-allelic. Different codon-index? No, same. Different alt allele. We skip multi-allelic candidates entirely; that's a separate concern (see varcode.transforms.split_multiallelic — explicit transform for multi-ALT split #370 discussion — closed as not-planned).

5. In-degree analog for buckets: a codon bucket with 4 SNVs is structurally suspicious (only 3 positions per codon). If we see this, log a warning and leave them all unchanged.

Commit plan

1. Scaffold (~80 LOC): function signature, codon-bucket helper, pass-through paths, no actual combining yet. Tests for pass-throughs.
2. Body (~250 LOC): _all_cis, _build_combined_mnv, gap fill via reference_range, metadata merge, multi-transcript dedup.
3. Tests (~400 LOC): full matrix from above.
4. Docs (~60 LOC): docs/transforms.md section, API reference entry, roadmap → empty.

Open questions

• Should the same-position multi-allelic case raise, warn, or silently skip? Leaning warn + skip (matches pair_breakends's "in-degree > 1" handling).
• Should combined_from_cis_snvs be a flag or a count? Flag (boolean) is simpler; the count is recoverable from len(source_variants).
• Worth supporting combine_cis_snvs without a phase_resolver (homozygous-only path) as the default? Yes — keeps the simple case ergonomic.