fix(hgvs): order the anchor candidates by strand, so a periodic reference resolves to the span the variant touches (#96)

[kuanlinhuang]

Closes #96, taking option 1 from the issue: pass the strand, because the coordinate convention already determines which end of the span protein_start names.

The residue #93 left

#93 gave hgvsp_inframe_indel a second anchor to try when the peptide does not corroborate the first. That fixed every case where the first end is not a match at all - 99.7% of the declines. It left the case where both ends corroborate: find() takes the first, and for a shrinking change on the reverse strand the first is the wrong one.

The wrong answer is not another spelling of the right one. With peptide MEGEGEA, EGE sits at 2-4 and again at 4-6; deleting 2-4 leaves MGEA and deleting 4-6 leaves MEGA, so no later 3'-shift can reconcile them. The description names residues the variant does not touch, which is the family #91 and #58 were both about.

What decides it now

protein_start comes from the genomic left edge. The residues come from the lower of the two CDS coordinates for a shrinking change and from cds_start otherwise (predict_coding_consequence), so it is exactly the shrinking change on the reverse strand whose span arrives anchored at its end. anchored_at_span_end states that as a predicate and anchor_candidates orders the pair by it, instead of leaving the choice to iteration order.

The other end stays as a fallback rather than being dropped. Over 37,122 in-frame ClinVar rows the residues sat at protein_start in 71.1% of cases, one residue earlier in 11.9%, and at neither in 17.0% - so ordering is all the strand can honestly buy, and the peptide still decides. A reverse-strand caller whose residues really do sit at protein_start still normalises; there is a test for exactly that.

Both routing copies pass tr.strand: the CLI pipeline, and the one in fastvep-annotate that serves the web path. Nothing structurally keeps those two in step, so each has its own end-to-end test.

Measured end to end

The real-data gate is in this comment: over the whole 673,660-variant ClinVar 2-star corpus, 15,684,064 consequence rows, 6 HGVSp rows change and nothing else does - Protein_position and every other CSQ field are byte-identical. All six are ClinVar 421481 (PCDH15, reverse strand, PAP/-), verified against a peptide recomputed independently in Python.

Before that corpus was located, the same question was answered synthetically, and the synthetic run is still the broader of the two: 8,000 synthetic in-frame deletions, 4,000 per strand, each on its own single-exon gene with a real translated peptide, driven through fastvep annotate --hgvs on both binaries. Scoring applies each emitted description back to the reference peptide and compares the result to the protein the variant actually produces - the invariant a well-formed-but-wrong span breaks.

	correct	wrong protein
fd982f4 (master)	7,985	15
this branch	8,000	0

15 descriptions changed, every one on the reverse strand, every one from wrong to right. No forward-strand description changes at all. Examples, as peptide del<span>:

MGAGAGAE   del3-5   p.Ala5_Ala7del  ->  p.Ala3_Ala5del
MKAKAKGK   del2-4   p.Lys4_Lys6del  ->  p.Lys2_Lys4del
MEEEGEEGEAEAAG del3-6  p.Glu6_Glu9del -> p.Glu4_Glu7del
MEAGAGAAGAA del3-5  p.Gly9_Ala11del ->  p.Ala3_Ala5del

The last one is worth noting: the old answer had been 3'-shifted off the wrong anchor onto residues 9-11, which really are GAA in that peptide. Nothing about it looks wrong.

The 0.4% rate on reverse-strand rows is inflated by design - the fixtures use a four-residue alphabet so periodic references are common, which is what made the synthetic run useful before the real corpus was to hand. On real ClinVar data the rate is 6 rows in 15.7M: 2,453 in-frame rows satisfied the necessary condition for the ambiguity and 2,447 of those either failed to corroborate both anchors or converged after the 3'-shift.

A differential run over vep_example_GRCh38.vcf, clinvar_inframe_deletions.vcf, chr22_1kgp.vcf and tests/test.vcf against the in-repo fixture annotation is byte-identical too, but that one is weak evidence and worth labelling as such: the fixture FASTA has no chr17 sequence, so those runs carry no protein descriptions at all. The ClinVar run above is the one that counts.

Tests

10 new, each confirmed to fail with the ordering removed ([Some(protein_start), from_end] unconditionally):

a_periodic_reference_on_the_reverse_strand_names_the_span_the_caller_meant
a_periodic_reference_on_the_forward_strand_is_read_from_the_start        (guard)
a_two_residue_homopolymer_reads_the_same_from_either_end                 (guard)
only_a_shrinking_change_on_the_reverse_strand_is_anchored_at_its_end
anchor_candidates_puts_the_determined_end_first
every_description_reconstructs_the_protein_the_variant_produces
hgvsp_reads_a_periodic_reference_from_the_end_the_strand_determines      (fastvep-annotate)
hgvsp_reads_the_same_periodic_reference_from_the_start_on_the_forward_strand
cli_reverse_strand_deletion_names_the_residues_it_deletes                (run_annotate)
cli_forward_strand_deletion_of_the_same_residues_is_unchanged

The five that assert the new behaviour fail without it; the four guards pass either way, which is their job - they exist so the fix cannot trade one wrong end for the other.

every_description_reconstructs_the_protein_the_variant_produces is the one worth keeping an eye on: 4,000 shrinking changes from a seeded LCG, each applied back to its peptide. It found a shape the issue did not mention - shrinking delins, not just pure deletions:

P:p.Lys6_Lys7delinsGly does not reconstruct EGKAGKGEGAKE from EGKAKKKGEGAKE
  (KK/G at 5-6, anchor 6 on Reverse)

The existing real-data tables (test_hgvsp_inframe_indel_never_emits_substitution_shape and the terminal-insertion case) now run under both strands, which pins that insertions are anchored at protein_start either way.

941 tests pass, up from 931. clippy --all-targets -D warnings clean, rustfmt clean.

Not in scope

The E2E fixtures deliberately reproduce the bug through run_annotate rather than by handing protein_start to the unit under test: the anchor convention is a property of the coordinates the pipeline produces, and a unit test that passes the anchor by hand asserts the author's reading of the convention rather than the pipeline's.

🤖 Generated with Claude Code

[kuanlinhuang]

The ClinVar gate, run

The corpus was on this machine after all - data/benchmark/clinvar_2star.vcf, the same 47,013 indels #91 and #93 measured against, annotated with test_data/organisms/human/Homo_sapiens.GRCh38.115.gff3 + the GRCh38 primary assembly. Both binaries, same V3 transcript cache built once with sequences, --hgvs.

Whole corpus: 673,660 variants, 15,684,064 consequence rows

	rows differing
HGVSp	6
Protein_position	0
every other CSQ field	0

Six rows. All six are the same variant - ClinVar 421481, 10:53822439 CAGGAGCAGG>C, a likely-benign in-frame deletion in PCDH15 - seen through its six protein-coding transcripts:

ENST00000320301  pp=1760-1762  PAP/-   p.Pro1762_Pro1764del -> p.Pro1760_Pro1762del
ENST00000361849  pp=1762-1764  PAP/-   p.Pro1764_Pro1766del -> p.Pro1762_Pro1764del
ENST00000373957  pp=1767-1769  PAP/-   p.Pro1769_Pro1771del -> p.Pro1767_Pro1769del
ENST00000395430  pp=1757-1759  PAP/-   p.Pro1759_Pro1761del -> p.Pro1757_Pro1759del
ENST00000395433  pp=1737-1739  PAP/-   p.Pro1739_Pro1741del -> p.Pro1737_Pro1739del
ENST00000437009  pp=1691-1693  PAP/-   p.Pro1693_Pro1695del -> p.Pro1691_Pro1693del

Protein_position byte-identical across all 15.7M rows, so the gate the last three of these had passes. Every other field too - the only thing this branch moves is the HGVSp of those six rows, and each moves onto the span Protein_position already named.

That variant, verified without our code

PCDH15 is reverse-strand and the reference is PAP, which is periodic with period 2, so both anchors are corroborated and the old resolution was a coincidence. Recomputing the peptide independently in Python from the GFF3 CDS blocks and the FASTA (32 CDS blocks, reverse strand, 5,868 bases, 1,956 residues):

residues 1752-1772:  PISPPSPPPAPAPLAPPPDIS
                              ^^^ ^^^        PAP at 1760-1762 and again at 1762-1764
deleted CDS 5278-5286: CCT GCT CCT  =  P A P   (= residues 1760, 1761, 1762)

del 1760-1762 -> PISPPSPPAPLAPPPDISP     the protein this variant makes
del 1762-1764 -> PISPPSPPPALAPPPDISP     the protein the old description named

Different proteins, so the 3'-rule cannot reconcile them: only one of the two descriptions belongs to this variant, and the branch now emits it. ClinVar's own record agrees on the event (CLNHGVS=NC_000010.11:g.53822445_53822453del, MC=SO:0001822|inframe_deletion).

An Ensembl VEP cross-check is not included: VEP refuses --gff --fasta --offline without a species cache directory, and pulling the ~25 GB homo_sapiens 115 cache for one row was not worth it. The independent peptide recomputation above settles the same question without depending on VEP's shift behaviour, which we already know diverges near the terminus (#94).

How exposed the class is, on real data

From the 47,013-indel subset (1,203,053 consequence rows, 80,681 in-frame):

in-frame rows                                          80,681
  on the reverse strand                                35,579
    shrinking                                          26,599
      reference length >= 3                             8,191
        first residue == last (ambiguity possible)      2,453
          resolved differently before this branch           6

The last two lines are the answer to the question the issue asked first. 2,453 rows could have hit it - reverse strand, shrinking, n >= 3, and the one-residue overlap between the two candidate spans satisfied. In 2,447 of them the two anchors either did not both corroborate or converged after the 3'-shift. Six did not.

So the issue's own estimate was close: "the answer might be twice". It is six rows on one variant, and #96's option 3 - document it and move on - would have been a defensible call at that rate.

Two things still argue for the fix as written. It removes the class rather than bounding it, using information already at the call site: the strand was always the determinant, and the old code was reading a coincidence instead. And the property sweep found the same fault in shrinking delins, which the issue's analysis did not cover - the exposed set there is not this 2,453 and has not been measured, because it needs the peptide-carrying corroboration check to enumerate rather than the necessary condition above.

941 tests pass, CI green, and the five behaviour tests still fail with the ordering removed.