Design and implementation contract

This document describes the implemented DuckVEP architecture. Code and tests define executable behavior; open compatibility scope and evidence are linked in the relevant sections.

Purpose and compatibility

DuckVEP is the deterministic variant-consequence kernel inside DuckHTS. It targets Ensembl VEP 116 (57ea5c52340acc1f156267f810ad162e26597082) with Ensembl core c0cf13daa961d80584bad797b2eb0ff3a7500ef3 and variation 2fb834b987ede3824e200197a838ce11e91aeb4b. The behavioral target is not human-only: the kernel accepts every NCBI codon table supported by VEP 116’s BioPerl translator. The declared full-model parity matrix covers human GRCh37 and GRCh38; non-human parity covers Plasmodium falciparum from Ensembl Genomes 63 paired with VEP 116, including codon tables 1, 4, and 11. Evidence for one species or assembly is never treated as evidence for another.

The kernel emits structured transcript consequences. VEP-compatible CSQ text and HGVS are edge projections, not the internal representation. Prepared-event and CDS-edit facts are the shared consequence/HGVS authorities; an HGVS-facing transcript-edit carrier adds VEP’s clipped transcript-slice state without changing those semantic edits. Allocation-free c./n./p. fact/render kernels, a worker-owned indexed-FASTA provider, and the query(duckvep_annotate_sql(..., struct_pack(hgvs := true))) relation path are implemented for independent literal small variants. Strict executable-VEP evidence covers the fixed position-one/right-anchor cases and 56,998 chromosome-21 ClinVar transcript pairs with zero HGVSc/HGVSp differences. Genomic HGVS, transcript models whose sequence differs from the genomic exon sequence, structural/BND HGVS, broader cross-species HGVS distributions, and compound/phased HGVS remain open. VEP 116 is the behavioral authority; pure-C properties and bcftools csq supply independent checks for mechanics and phased edit state.

DuckVEP is alpha: its own interfaces and intermediate representations have no backward-compatibility requirement. Replace or delete them when a simpler ownership contract or shared semantic authority makes them unnecessary; do not retain forwarding APIs or alternate implementations solely for older DuckVEP callers. Preserve the pinned oracle, physical input records, comparison denominators and failure controls; document deliberate result changes and test missing evidence explicitly instead of retaining a second implementation.

Mental model

DuckVEP does two different jobs at two different times.

First, it compiles Ensembl core and funcgen relations plus matching reference sequence into an immutable execution model. DuckDB keeps the rich source relations, stable identifiers, attributes, and hashes. The compiler checks them, assembles transcript-oriented sequence, and publishes compact C arrays containing only facts used repeatedly during annotation.

Second, it runs a sorted stream of variant alleles through that model. Each allele is interpreted once, matched to candidate transcripts, projected through exon and CDS coordinates, and reduced to consequence facts. One generated VEP-116 rule program maps those facts to SO terms. Numeric rows leave the kernel first; identifiers, HGVS, supplementary annotations, and strings are joined or rendered later.

model construction, once per declared source

Ensembl core + funcgen relations + matching FASTA
                -> checked prepared relations + deterministic receipt
                -> immutable transcript/exon/sequence + regulation/motif arrays

annotation, repeatedly

sorted ALT alleles
    -> raw VCF span + VEP feature span + minimized sequence edit
    -> continuing transcript and regulation/motif sweeps
    -> topology / splice / projection / sequence / interval-feature facts
    -> one VEP-116 fact-to-SO program
    -> compact consequence rows keyed to transcripts or core interval features
    -> SQL joins, identifiers, HGVS, CSQ/JSON, clinical evidence

The data placement follows the cost and lifetime of each fact:

Data Where it lives and why
Ensembl identifiers, attributes, source tables, and receipts DuckDB relations, because they are provenance and joinable data rather than hot-loop state.
Transcript spans, exon maps, prepared CDS/flanks, sparse sequence edits, and compact regulatory/motif intervals One immutable named C model, because every variant reuses them.
Variant alleles Sorted batches, because genome-scale inputs stream and coordinate order makes candidate discovery nearly linear.
dbSNP, ClinVar, gnomAD, scores, and interval tracks DuckDB/Parquet streams or bounded tiles, because these sources can dwarf the transcript model and may be selected independently.
Exon cursors, candidate lists, sequence scratch, and output buffers One workspace borrowed exclusively by each active scalar callback, because mutable state must not be shared across concurrent callbacks or named models.

Five rules keep the design understandable:

  1. There is one biological authority for independent consequences: fact producers followed by the generated VEP-116 rule program. A fast path must prove equality with it.
  2. The uploaded VCF span, VEP’s feature span, and the minimized sequence edit are retained separately because VEP predicates genuinely consume different coordinates.
  3. Coordinate order is part of the execution contract. cgranges seeds a run; a continuing sweep, not a fresh transcript search per allele, handles the common sorted workload.
  4. Missing source facts stay explicit. A partial model, unsupported sequence edit, invalid projection, or unknown reference is not silently converted into a supported consequence.
  5. Compatibility and speed are measured separately. VEP differentials compare the union of all emitted variant/object pairs; throughput reports count input alleles, candidate work, emitted rows, bytes, model, threads, and materialization.

The continuing sweep, concretely

A per-allele interval lookup would repeatedly search the same chromosome-sized transcript index. DuckVEP instead exploits the public input contract: model rows and ALT alleles are both nondecreasing by model-local contig and coordinate. At the first allele of one native batch, cgranges supplies every transcript whose span plus directional flank can contain the event. The workspace then retains three monotone frontiers:

  1. the next transcript start not yet admitted;
  2. the active transcript indices not yet expired behind the current event; and
  3. the exon rank last reached for each active transcript on the point and general-feature projection paths.

For each later allele, the worker advances the admission frontier while transcript starts are reachable, compacts expired transcripts out of the active set, and projects only the survivors. A wide allele may additionally admit transcripts whose starts fall inside its full event span; those candidates are unioned with the continuing point active set. The regulation/motif relation uses the same general interval-candidate mechanics but a separate active set and exact event overlap, because it has no upstream/downstream transcript flank.

sorted model spans:       [T1---------] [T2---]       [T3-----------]
sorted event starts:          e1  e2        e3  e4             e5
                              ^ seed once
                                 admit newly reachable starts ->
                                 retire ends left of the event ->
                                 classify the remaining active set

Ignoring emitted pairs, one monotone run is O(transcripts admitted + events) rather than O(events * log(transcripts)) independent searches. Dense regions still cost real work: if an event overlaps many transcripts or core features, DuckVEP must emit and classify those candidates. There is no reward for selecting sparse windows, and the performance suite therefore includes transcript-, exon-, non-coding-RNA-, regulation-, motif-, and observed-variant-dense tiles.

The upstream_distance and downstream_distance values change admission geometry only. They are not allocation sizes and do not clip the event. Statistical sweep scenes exercise zero, 1, 50, 100, 4,999, 5,000, 5,001, 10,000, 50,000, and 65,535-base distances, including alleles wider than the requested flank. The current scalar adapter has no guaranteed state across DuckDB vector edges, so it seeds each native vector and continues inside that vector. A stateful table-function adapter can carry the same frontiers across vectors only when its host lifecycle exposes explicit carry state; it must not create a second candidate-selection authority.

The resident model is immutable and shared among DuckDB workers. Each checked-out workspace owns its active/candidate arrays, exon cursors, result builder, reference window, and faidx handle. Partitioning therefore multiplies bounded workspace scratch rather than the whole transcript/sequence model, and no mutable iterator or reference cache crosses threads.

Phased haplotypes and structural variants reuse projection, sequence editing, and the SO rules, but they are not disguised as independent small variants. Haplotypes group several edits before translation; breakends retain two loci. Supplementary annotation and ACMG/AMP reasoning remain relational work above the deterministic consequence kernel.

Responsibility split

DuckDB owns:

The C kernel owns:

htslib owns VCF/BCF transport and genotype decoding. HGVS consumes the same projected edit facts as consequence prediction; it must not reconstruct biology from rendered SO strings.

Code map

Transcript model build

The model builder is a compiler over relations already staged in DuckDB. It is not a downloader, a MySQL client, or a parser for Ensembl dump archives. The caller first imports the pinned Ensembl tables and supplies matching reference sequence chunks. Building and loading are then separate operations:

Ensembl core tables ─────┐
Ensembl funcgen tables ──┼─> validate and prepare relations ─> receipt
reference FASTA ─────────┘                 │
                                           ├─> region projection
                                           ├─> transcript/exon projection
                                           └─> regulation/motif projection
                                                   │
                                                   v
                                      immutable named C model
                                                   │
                                                   v
                                      sorted consequence queries

This separation is deliberate. DuckDB performs the large joins, sequence assembly, and provenance work once. Annotation workers reuse compact immutable arrays and do not carry a Perl object graph, stable-ID strings, or source-table metadata through the hot loop. Preparation and receipt SQL is emitted by stable-C-API scalar builders and evaluated through query() in the caller’s connection. The C registration file does not iterate transcript rows or implement a second importer.

Inputs

query(duckvep_ensembl_regions_sql(...)) and query(duckvep_ensembl_transcripts_sql(...)) read these Ensembl core relations by name from the supplied schema:

query(duckvep_ensembl_regulation_features_sql(...)) reads the release-matched funcgen regulatory_feature, feature_type, and motif_feature relations. Funcgen uses core sequence-region IDs; the generated SQL joins them to the already prepared region relation, rejects missing feature types or invalid coordinates, and assigns one dense ordinal space across RegulatoryFeature and MotifFeature rows. It also reproduces VEP 116’s source selection by discarding epigenetically_modified_region (EMAR) RegulatoryFeature rows before ordinals are assigned; VEP’s database annotation source removes those rows before constructing overlap objects. Stable IDs, feature types, regulatory-build IDs, binding-matrix IDs, strand, and scores remain cold DuckDB columns. Only region ordinal, start, end, kind, and feature ordinal enter the C model.

Attribute values are cast to their semantic text form while importing these relations. This matters for official dumps where a generic CSV/Parquet staging pass may infer a numeric-only value column as an integer; the model compiler does not require callers to rewrite an otherwise valid staged schema solely to satisfy string predicates.

The reference relation has (chrom, start, end, seq), where start is zero-based and end is half-open. It is normally persisted from tiled fasta_nuc(..., include_seq := true) output. A reduced FASTA deliberately builds a reduced model; the builder never silently fills absent contigs from another source.

GRCh37 is a separate Ensembl source, not a coordinate option applied to the GRCh38 genebuild. Release 116 publishes homo_sapiens_core_116_37 on the GRCh37 archive; it uses the release-116 schema around the frozen release-75/GENCODE-19 annotation. It has GENCODE attributes but no MANE data, because MANE is defined only on GRCh38. A GRCh37 model must therefore use matching GRCh37 core relations and reference sequence and must not invent a MANE mapping. See https://grch37.ensembl.org/index.html and https://www.ensembl.org/info/genome/genebuild/mane.html.

Region preparation

query(duckvep_ensembl_regions_sql(...)):

  1. verifies that every reference chunk is non-null, has valid coordinates, and contains exactly end - start bases;
  2. verifies that chunks cover each supplied contig continuously from zero;
  3. matches each FASTA contig to exactly one same-name, same-length Ensembl region on the requested assembly and species; and
  4. assigns dense model-local region ordinals from zero while retaining the Ensembl source identifiers and coordinate-system fields.

The model contract bounds a region coordinate at UINT32_MAX and the number of regions at 65,536. A mismatch is a query error; no partial region model is published.

Transcript and sequence preparation

query(duckvep_ensembl_transcripts_sql(...)) applies the VEP-116 core-source filter before assigning model ordinals: a transcript must be current and have a non-empty stable ID, and neither the artifact biotype nor a readthrough_tra attribute is admitted. This is why a model built from the full core dump has the same candidate transcript population as VEP’s core cache instead of merely the same coordinate source. For each selected transcript it:

  1. joins its gene, optional translation, ranked exons, and selected attributes;
  2. assigns dense transcript and gene ordinals while retaining Ensembl numeric IDs, stable IDs, versions, and biotypes in the prepared relation;
  3. calculates each exon’s transcript-oriented cDNA span from its rank;
  4. projects the translation start and end from exon-relative coordinates into genomic and cDNA coordinates;
  5. reconstructs exon sequence from the reference chunks, reverse-complements negative- strand exons, joins exons in transcript order, and extracts the CDS;
  6. applies Ensembl start-phase preparation: phase -1 or 0 adds no prefix, while a positive phase adds one or two leading N bases; and
  7. reads the sequence region’s codon_table attribute, defaulting to table 1 exactly as VEP does, and rejects conflicting, malformed, or unsupported table IDs;
  8. extracts the complete transcript-oriented spliced sequence before and after the CDS;
  9. parses supported initial_met, _selenocysteine, amino_acid_sub, and _stop_codon_rt Translation SeqEdits into a sparse reference-peptide relation; and
  10. parses each mature-miRNA cDNA range from the Ensembl miRNA transcript attribute and projects it through the ranked exons into one or more genomic segments.

The compact flag word records translation presence, protein-coding/NMD/miRNA biotype, incomplete CDS ends, selenocysteine, stop readthrough, transcript RNA edits, peptide edits, MANE, GENCODE, CCDS, and upstream-start state. Translation _rna_edit is not one of Translation::get_all_SeqEdits() in VEP 116 and is therefore not treated as a peptide edit. The standalone model ABI still defines a readthrough-transcript flag for explicitly prepared alternate models, but the VEP-compatible core importer filters those transcripts before publication.

For GRCh38 MANE attributes, the prepared DuckDB relation also retains the attribute value as mane_select_refseq or mane_plus_clinical_refseq. That value is the paired versioned RefSeq transcript accession. The builder rejects multiple or empty mappings for a MANE flag. Only the selection bits enter the resident C model; Ensembl/GENCODE and RefSeq identifiers stay in the cold relation for late SQL projection.

Single-position, single-amino-acid Translation SeqEdits are kept with the reference-derived CDS and applied only to VEP’s reference peptide; the alternate peptide remains the raw codon translation. Transcript-level sequence corrections and other Translation SeqEdit shapes withhold sequence with an explicit reason. In particular, the current exon map requires every transcript cDNA exon span to be contiguous and the same length as its genomic exon span; inserted or deleted transcript bases need a richer mapping. The transcript, coordinates, and flags remain in the model. Unsupported reference alphabet has the same fail-closed shape. This is different from malformed topology: bad coordinates, strand, exon phase or rank, translation bounds, or incomplete sequence reconstruction abort the build.

The resident model stores CDS bytes and both non-coding transcript flanks in two packed byte pools with offsets and lengths per transcript. The flank pool has no per-transcript allocation or padding and is cold on ordinary coding rows. It is read only when a VEP predicate rebuilds the 5-prime-UTR-plus-CDS or CDS-plus-3-prime-UTR string. The prepared relation and resident model use the complete flanks; no short-tail projection is stored.

Peptide edits and mature miRNA ranges are prepared once, not remapped for every variant. Peptide edits are packed by transcript and protein position and consulted only for an overlapping reference-peptide window. Ensembl stores mature-miRNA ranges as transcript cDNA intervals, possibly more than one per transcript. The builder splits a range when it crosses an exon boundary and returns the resulting inclusive genomic segments in mature_mirna_regions. The resident model packs all segment starts and ends into flat arrays with one offset per transcript. A non-miRNA transcript pays only the transcript-flag test; a miRNA candidate scans its usually tiny owned slice.

Prepared relations and publication

The builder returns one row per transcript and one row per core regulation/motif feature. The transcript row contains the hot fields, source and stable identifiers, biotypes, the optional prepared CDS, an ordered nested exon list, and nested mature-miRNA genomic segments. The feature row contains a dense ordinal, interval, kind, and cold source metadata. Callers persist the region, transcript, and feature relations and derive the required sorted loader projections plus optional side projections from them:

The core_schema argument names relations, not a transport. It can point at tables loaded from Ensembl’s tab-separated MySQL dumps, or at a read-only MySQL catalog attached through DuckDB’s mysql extension. Downloading, attaching, and staging stay outside the model builder so extension builds remain offline and the same validation runs for either source. The builder does not require a MySQL server once those relations and the matching reference chunks have been persisted.

Versioned model input contract

A DuckVEP model is identified by the tuple (source, source release, VEP behavioral release, species, assembly, transcript-selection policy, reference identity). The source release and behavioral release are separate fields because an adapter can deliberately support a bounded VEP behavioral subset over a particular Ensembl source release; they must never be inferred from a cache filename or silently changed in place.

The release adapter is the only version-specific layer. It obtains the declared Ensembl relations, maps release-specific schema or attribute-policy differences into the canonical input relations below, and records the exact source manifests. It must not reinterpret canonical model fields, alter ordinal ordering, or make a version-specific selection without putting that selection in the receipt and validation evidence.

Builder Required canonical source relations
duckvep_ensembl_regions core_schema.coord_system, core_schema.seq_region, core_schema.seq_region_attrib, core_schema.attrib_type, and a reference-chunk relation with chrom, zero-based start, half-open end, and seq
duckvep_ensembl_transcripts region inputs plus core_schema.seq_region_attrib, transcript, gene, translation, exon_transcript, exon, transcript_attrib, translation_attrib, and attrib_type
duckvep_ensembl_regulation_features funcgen_schema.regulatory_feature, feature_type, and motif_feature, plus the canonical region relation

The extension validates the columns it reads and rejects inconsistent values; the exact source-column projections are executable in the SQL builders and acceptance fixtures. A release adapter that adds source handling must produce the same canonical region, transcript, and regulation relations. The loader accepts an 11-column CDS-only projection or a 13-column complete-flank projection. These express different available sequence evidence; missing transcript flanks remain explicitly unresolved when required. A changed model contract requires an explicit receipt and tests. A receipted artifact retains the meaning of its declared contract.

Promoting a new VEP target therefore requires three independent proofs: a pinned public source/release manifest and reference identity, canonical-model receipt validation, and release-specific differential cases against the pinned upstream VEP target. The evidence must state the supported consequence/HGVS/CSQ subset and every intentional difference; passing VEP 116 evidence does not validate another release.

The prepared region relation contains sequence_length and a Boolean circular sourced from the same Ensembl core release as its seq_region. The builder accepts one circular_seq attribute with value 1 per circular region and rejects duplicate or invalid circular attributes; absence means linear. Model identity is model_sha256, which covers the declared model rows. topology_sha256 separately hashes region name, length and circular flag, and the receipt lists circular regions. A change to hashed rows requires a versioned definition and re-recording every receipt in the same commit. Circular topology is a property of the reference sequence region; codon_table is an independently sourced translation rule. A mitochondrial codon table does not imply that a region is circular, and a circular region does not imply mitochondrial translation. Ordinary human Ensembl-116 MT transcripts in the acceptance fixture do not cross the origin. Circular preparation preserves one-based inclusive coordinates, including start > end spans. origin_crossing and circular are explicit output columns; exon order is determined by transcript rank on either strand. A valid wrapped transcript crosses the origin exactly once across its ranked exons, with the first and last exon anchored to its transcript endpoints. Wrapped single exons use the reference tail followed by its head; regulation and motif features retain their original endpoints. The model loader accepts a four-column region query (seq_region, sequence_length, seq_region_name, circular) and validates wrapped coordinates, exon lengths and rank continuity against topology. Ordinary linear models retain their validated kernel and query behavior.

Circular-coordinate execution (lifted intervals). This is the circular-topology feature; it is unrelated to the mitochondrial codon table, which only selects a translation rule. The consequence kernel stays linear: coordinates are unsigned 32-bit positions and start > end there means an insertion, never a wrap. A circular region that carries at least one wrapped object (a transcript, exon, regulatory feature or motif feature with start > end) is therefore executed on a second, fully linear model built once when the model is loaded (src/kernel/src/duckvep_lift.c):

A region is lifted only when it is circular and holds a wrapped object. A circular region without one, such as human Ensembl-116 MT, runs on the linear kernel and matches the output fingerprint pinned by test/sql/duckvep_circular_mt.test; it reports no upstream or downstream consequence across the origin, exactly as VEP does. Lifting is a property of model contents, not of the circular flag alone. Model loading accepts wrapped models; structural, breakend and phased edit-set entry points reject them with an explicit error because their lifted semantics are not defined. 2B + 3L must stay below 2^31 - 1. HGVS 3’ normalization on a circle shorter than about 2 kb sees the sequence repeat, because the 1000-base shift window on each side wraps onto itself; VEP defines no behavior there.

Validation: make test_properties (and test_properties_sanitized) run native properties over random circular worlds: rotating reference, model and events by random offsets leaves every consequence, projected position, peptide, NMD prediction and escape reason unchanged, one row per event/object holds, and away from the origin the lifted result equals the linear kernel. test/sql/duckvep_circular_lifted.test does the same through SQL with HGVS and projected edits for five rotations, both insertion interbase orientations, CDS starts and ends, exon-intron junctions on both strands at the origin, MNVs, long alleles, four worker threads and several DuckDB vectors, and compares every object, wrapped or not, with an ordinary linear model in which a rotation puts that object mid-sequence. benchmarks/duckvep_circular_origin.py records throughput, output equality and peak memory on an origin-focused workload.

Evidence status. Public origin-crossing transcripts exist and were compared with executable VEP, but VEP is not an oracle for them. Ensembl release 116 has none (only fly and yeast MT carry circular_seq, with no inverted transcript); Ensembl Genomes 63 has 19 inverted transcripts on circular regions, 18 in bacterial and archaeal collections and one trans-spliced plastid rps12, three of whose genomes have a VEP cache. On those genomes (benchmarks/data/circular_vep_differential/) DuckVEP equals VEP on every SO term of the other 152, 937 and 553 transcripts, except flank rows that exist only through the origin, and HGVS 3’ shifts within 1,100 bases of it, where VEP clips its window at the sequence end; comparison with the non-lifted extension baseline is exact except for those origin-reaching rows. VEP models the crossing transcript itself as an interval with reversed bounds (no row for most events inside it, or an intergenic_variant transcript consequence) and aborts --hgvs on some events in its translation. That is a difference in what is modelled, not agreement or refutation. Circular-coordinate support for a crossing object is therefore property-proved and linear-model-proved, not oracle-proved: rotation equivariance, equality with an ordinary linear model per object, and the native properties. The survey, the differential and every disagreement are in benchmarks/data/circular_source_survey.md and ERRATA.md.

query(duckvep_model_receipt_sql(...)) checks dense ordinals, region/transcript agreement, and every regulatory/motif interval against its declared region. It records the declared source, release, assembly, transcript filter, source-manifest hash, reference hash, model counts including CDS, transcript-flank bases, mature-miRNA transcripts and projected segments, peptide edits, regulatory regions, and motif features, and a deterministic hash over every hot model field, including mature-miRNA ranges, reference-peptide edits, and interval-feature geometry. There is no timestamp: identical declared inputs must produce the same receipt.

The checked-in acceptance fixtures under test/data/duckvep/ensembl_core/ are about 120 KiB. The GRCh38 fixture contains complete release-116 MT, KI270395.1, and HG2047_PATCH source rows; it covers mitochondrial codon-table and peptide-edit behavior, an ordinary multi-exon CDS, three real release-116 MotifFeature rows on KI270395.1, and the real ENST00000715685 ↔ NM_032790.4 MANE pair. The GRCh37 fixture contains MT and GL000201.1; it proves sequence-backed coding annotation from the archived GENCODE-19 model and the absence of MANE mappings. The explicit staging script verifies both official core manifests, the GRCh38 funcgen manifest, assembled reference hashes, deterministic model receipts, and exact model counts before writing Parquet. The component manifest hashes are folded into one sorted canonical source-manifest hash. Tests never contact Ensembl.

GRCh37 transcript selection and external MANE mapping

The GRCh37 consequence model uses Ensembl 116 GRCh37 core records without native MANE. For caller selection, scripts/build_mane_grch37.R writes two cold Parquet relations: grch37_transcript_authorities.parquet contains the source gene’s retained canonical transcript state and independent GENCODE-19 Basic flag for every filtered model transcript; mane_grch37_mapping.parquet audits each row of MANE v1.5 against the NCBI GRCh37.p13 annotation. Neither relation changes native consequence masks, ordinals, HGVS, or VEP-parity results. There is no GRCh37 GENCODE Primary or native MANE flag.

Run sh scripts/stage_mane_grch37.sh STAGING_DIR to acquire resumable, checksum-pinned inputs, then Rscript scripts/build_mane_grch37.R STAGING_DIR MODEL.duckdb OUTPUT_DIR. The recipe requires curl, samtools, and R packages DBI, duckdb, data.table, Biostrings, Rsamtools, GenomicRanges, and digest. The script uses R FASTA readers and plain DuckDB SQL; it loads neither DuckHTS nor DuckVEP into its DuckDB connections. The staging directory holds the MANE v1.5 summary, NCBI GCF_000001405.25 assembly report/GFF/RNA/protein FASTAs, and Ensembl GRCh37 primary-assembly FASTA (compressed and indexed uncompressed copies). The script checks pinned SHA-256 digests and the immutable model receipt; it records source URLs and hashes in mane_grch37_receipt.csv. The model’s core tables came from Ensembl’s public MySQL ensembldb.ensembl.org:3337/homo_sapiens_core_116_37, following test/scripts/prepare_duckvep_ensembl_fixture.sql. Its core source-manifest hash is carried from the independently content-verified model receipt because the original source dump manifest cannot be recovered from the compiled model. The FASTA receipt hash covers the compressed Ensembl FASTA, and the build checks the uncompressed indexed FASTA separately. All release inputs and outputs remain external; only the recipe, small policy tests, and receipt ledger are committed.

The mapping resolves exact versioned RefSeq nucleotide accessions in the NCBI GFF before considering the MANE ENST stable root as a candidate. MANE v1.5 contains 19,367 NM_ and 70 NR_ rows; the pinned GRCh37.p13 GFF contains 19,306 NM_ and 61 NR_ accessions. The 61 absent NM_ and nine absent NR_ rows remain rejected. The assembly report links the original NCBI accession to the target FASTA region; no implicit chromosome alias is evidence of transcript identity. For accessions annotated on multiple loci, the model’s contig restricts the candidate locus but cannot itself establish an association; the validation gates still apply. Strict association requires the exon chain in transcript orientation, CDS coordinates and phases, model spliced sequence versus target genomic FASTA, and RefSeq RNA translation versus its versioned protein. exact_model_match retains the full-transcript association. RNA substitutions relative to the target reference are reported explicitly; a RefSeq protein inconsistent with its RNA rejects the strict candidate.

cds_exact_utr_differs is a separate coding-only association. It requires the exact versioned RefSeq accession on the target reference, exactly one candidate with the same contig, strand, CDS segments and phases, a reference-genome CDS that translates to the pinned RefSeq protein, and zero RefSeq RNA/reference mismatches within the CDS. UTR or non-coding exon chains may differ. Evidence columns record true/false when checked and NULL when not evaluable; utr_exon_chain_match distinguishes full exon-chain identity from CDS identity, while cds_reference_difference_bases compares coding positions in the versioned RefSeq RNA against the target reference. If the RNA is absent or its length differs from the target spliced transcript, that comparison is NULL and the coding tier is not admitted. The protein gate uses RefSeq RNA for strict matches and target-reference CDS for coding-only matches. A substitution outside the CDS does not disqualify the coding tier. Ambiguous CDS candidates are not associated. Coding c./p. HGVS and coding consequences are equivalent for this tier; UTR c.-N / c.*N, n. positions, and UTR/non-coding exon consequences are not claimed. Other rows retain their reason, target accession, and candidate where available. The receipt counts statuses over all MANE rows and hashes the canonical ordered TSV representation of the full relation, independently of Parquet metadata. Run offline policy fixtures with make test_mane_grch37; validate an exported full release with make test_mane_grch37 MANE_GRCH37_OUTPUT=OUTPUT_DIR (the receipt test requires the external Parquet directory and the last row of the receipt ledger).

A caller joins the native cold relation by model SHA-256 and transcript_index, then joins mapped MANE rows by the same two fields after consequence expansion. scripts/mane_grch37_caller.sql exposes mane_mapped_to_grch37 for strict matches and mane_coding_region_only separately; the latter is usable only for coding consequences or coding c./p. HGVS. If several MANE rows reference one transcript, aggregate each tier into a list before joining rather than duplicating consequences. Canonical, GENCODE Basic, and mapped MANE are distinct source-attributed facts; selecting a representative is an explicit caller policy, not VEP --pick. If a RefSeq transcript is valid on GRCh37 but has no exact GENCODE-19 match, leave it unmapped or annotate it in a separately receipted RefSeq model.

duckvep_model_load(...) reads committed, non-temporary relations through a private connection, validates and narrows every value, builds independent transcript and regulation/motif seed indexes, and only then publishes the named immutable model. A failed load publishes nothing. Several models may coexist in one database instance, and their numeric ordinals are meaningful only within their model. Stable IDs and provenance remain DuckDB columns. This is also the contract for haplotype-resolved or pangenome paths: each model declares its exact assembly or path set and sequence hashes; contig aliases and mapping confidence remain explicit relations rather than being guessed from chromosome spelling.

Caller preparation may use duckhts_contig_key(...) on both an input contig and the model region name to construct that explicit relation. The key removes one non-empty leading chr prefix and normalizes only mitochondrial M/MT spellings to MT; it does not map numeric sex chromosomes, accessions, patches, or alternate loci. Callers must reject a model-side key collision before joining. The Ensembl model compiler continues to require an exact same-name, same-length reference-region match, so a convenience join key cannot silently substitute sequence from another assembly or region.

The transcript query accepts 11 CDS-only columns or the complete 13-column form ending in pre_cds_sequence and post_cds_sequence. The CDS-only form represents unavailable transcript flanks, not a short-tail approximation. Only the complete form may resolve length-changing edits crossing the CDS start or end. Incomplete sequence inputs return missing_transcript_flank when the required transcript bases are absent.

The optional mature_mirna_query has three columns: transcript ordinal, inclusive genomic start, and inclusive genomic end. Rows must be ordered by transcript and start. The loader proves that each range belongs to a miRNA transcript, stays inside one of its exons, and is ordered before publishing the model. Models without these Ensembl attributes omit the query; the loader does not invent mature regions.

The optional peptide_edit_query has transcript ordinal, one-based protein position, and one uppercase replacement amino acid. Rows must be unique and ordered by transcript and position. The loader packs them into the immutable sequence model and proves that every position lies within its prepared peptide before publication.

The optional interval_feature_query has five columns ordered by region, start, and dense feature ordinal: feature ordinal, region ordinal, inclusive start, inclusive end, and kind (1 RegulatoryFeature, 2 MotifFeature). The loader narrows these to a separate immutable SoA and builds its own cgranges seed index. Cold funcgen metadata is joined later by feature ordinal; it is not copied into each consequence row.

The loader treats transcript coverage as partial by default. Only a model deliberately loaded with transcript_coverage_complete := true may turn “no loaded transcript here” into supported intergenic_variant; a partial model returns an unresolved result instead.

CREATE TABLE reference_chunks AS
SELECT chrom, start, "end", seq
FROM fasta_nuc('GRCh38.primary.fa', bin_width := 1048576, include_seq := true);

CREATE TABLE model_regions AS
SELECT * FROM query(duckvep_ensembl_regions_sql(
  'ensembl_core', 'reference_chunks', 'GRCh38'
));

CREATE TABLE model_transcripts AS
SELECT * FROM query(duckvep_ensembl_transcripts_sql(
  'ensembl_core', 'reference_chunks', 'GRCh38'
));

CREATE TABLE model_regulation AS
SELECT * FROM query(duckvep_ensembl_regulation_features_sql(
  'ensembl_funcgen', 'model_regions'
));

CREATE TABLE model_receipt AS
SELECT * FROM query(duckvep_model_receipt_sql(
  'model_regions', 'model_transcripts',
  'Ensembl', '116', 'GRCh38', source_manifest_sha256,
  reference_sha256, 'VEP 116 core transcript selection',
  {regulation_features_table: 'model_regulation'}
));

The registered duckvep_ensembl116_model artifact has one clean-cache producer in r/duckhtsbench. Its registry rows pin the public Ensembl 116 homo_sapiens_core_116_38 and homo_sapiens_funcgen_116_38 CHECKSUMS, schema, and required table dumps plus the matching primary-assembly FASTA. The producer verifies each dump against its Ensembl manifest, verifies release 116, species homo_sapiens, species ID 1, and assembly GRCh38 from the imported meta and coord_system relations, and preserves the required source table names under ensembl_core and ensembl_funcgen. It then invokes the same public DuckVEP preparation and receipt SQL builders described above. The artifact records the source manifest, reference-sequence, and model hashes without a timestamp; a reused model must reproduce its stored receipt before provider exports are allowed.

What the builder does not implement

The extension’s SQL builder does not itself:

Those richer facts should remain typed DuckDB relations joined by numeric source IDs. They do not belong in every resident C transcript record. Exact VEP-compatible selection and the richer Ensembl relation set can therefore be separate named products built from the same staged release. The conformance report records the declared corpus evidence; throughput accounting remains tracked at https://github.com/RGenomicsETL/duckhts/issues/95. Further species and genetic-code coverage is tracked at https://github.com/RGenomicsETL/duckhts/issues/119.

Ownership

Independent-variant execution

query(duckvep_annotate_sql(events_table, model_name, struct_pack(hgvs := false, upstream_distance := 5000, downstream_distance := 5000))) is the public relation surface. events_table names a narrow, globally coordinate-ordered relation with one row per ALT allele: event identity, model-local region ordinal, one-based position, literal REF/ALT, nullable single-locus structural span and type/copy direction, and nullable mate coordinates. The relation validates that geometry, derives small, exact structural, or paired-breakend family, and dispatches to private native lanes with one fixed compact output schema. It never hides an ORDER BY; invalid input order is rejected by the native sorted-stream checks. Callers retain genotype, confidence interval, raw ALT, orientation, and other wide provenance in the source relation and join selected results back by event_index.

Both direction windows default to VEP’s 5,000 bases and zero disables the corresponding direction. These distances extend candidate admission only beyond transcript endpoints; they do not cap an allele span or clip an event that overlaps a transcript. Literal alleles up to the compact 65,535-byte slice limit retain uploaded, VEP-feature, and minimized-edit geometry independently. An ordinary minimized deletion that contains a complete transcript therefore reaches the same transcript_ablation fact as a symbolic deletion, while an equal-length containing span retains VEP’s endpoint-UTR comparison behavior. The private native adapter copies one DuckDB vector into compact arrays, splits on model/contig/window/order changes, seeds the first candidate set through cgranges, and advances independent sorted transcript and regulation/motif sweeps. The SNV point path keeps a per-transcript exon rank and advances it monotonically; normalized multi-base features keep a separate exon rank. A normalized-coordinate rewind within one vector uses a rewind-capable seek. If that vector was non-monotone, the workspace resets both rank arrays before its next vector: a transcript skipped after an earlier forward jump may still hold an ahead rank even when the next vector begins after the prior vector’s final event. Regulation/motif rows use exact event overlap with no transcript flank. Both sweeps share the generic interval-candidate helper, but own separate active sets because their cardinalities differ. Transcript fast/exhaustive paths and the complete feature sweep are property-checked against independent or brute-force oracles.

The scalar callback has no expression-local state across DuckDB vectors, so it restarts the sweep at each vector. It is a real batch interface, not the stateful whole-stream contract. The latter must expose explicit carry state through a stable host lifecycle rather than assuming DuckDB happens to preserve scalar callback state.

One consequence authority

The intended decision chain is:

  1. event geometry;
  2. transcript topology and projection;
  3. one edit/CDS/peptide context;
  4. consequence facts; and
  5. one static fact-to-SO rule table.

Specialized traversal or SNV codon code is acceptable only when an exhaustive/property oracle proves identical results. It is not a second biological authority.

Production MNVs and indels use the generalized edit/CDS/peptide context exactly once. A projection, sequence, capacity, or unsupported-state failure remains explicit and is never retried through a smaller shape-specific classifier. Narrow direct classifiers remain only as independent pure-C test references; the annotation path cannot select them after a context failure.

NMD prediction

NMD_transcript_variant remains the VEP core consequence for a variant inside a transcript already imported with the nonsense_mediated_decay biotype. It does not say that the current variant creates a new NMD substrate.

Variant-induced NMD is a separate compact result derived from VEP Plugins release/116 NMD.pm (0082591268417af618e03850c5ffdc7c09998a5d). Stop-gained, frameshift, splice-donor, and splice-acceptor consequences are predicted to escape for an intronless transcript, an early-CDS event (cds_end <= 101 in the plugin), an event in the last exon, or an event in the plugin’s inclusive 51-base penultimate-exon-end window. An eligible projected event matching none of those rules is triggering; an eligible event without coding coordinates is unresolved. The SQL result exposes each escape reason as a boolean.

The plugin does not use the minimized sequence edit for those coordinates. Its BaseTranscriptVariation::cds_coords path projects the full VEP VariationFeature and its exon rules read VariationFeature::seq_region_end. For an ordinary span, the first and last mapped coordinates may enclose an internal mapper gap as long as both endpoints map. For a pure insertion, the parent TranscriptVariation instead retains the empty feature as a reversed CDS range such as 102,101; one genomic flank may be sufficient at an exon edge. This is different from the expanded 101..102 allele range rendered in VEP JSON, and NMD.pm reads the parent’s lower cds_end. Immediately before the first coding base, the valid parent range is 1,0: the plugin tests whether both values are defined, not whether they are nonzero, and classifies the zero end as an early-CDS escape.

Equal-length alleles retain unchanged uploaded bases in their feature. DuckVEP must therefore preserve both views: the minimized edit changes the CDS, while the full feature decides NMD position. The early-CDS fact cached by the coding classifier is reusable only when both genomic spans are identical; insertions and wider features use the feature projector. This representation-dependent behavior is pinned by paired executable-plugin witnesses on both strands.

This is the pinned VEP positional policy, not a direct molecular assay or a claim that every transcript follows one universal NMD rule.

Phased edits

The pure C mutation core rebuilds a CDS from several non-overlapping edits in one reverse-coordinate pass, translates once, and partitions interactions while the frame is displaced or the next edit touches the same alternate codon. The model-scoped carrier index groups explicit (transcript, sample, phase_set, haplotype) keys and shares event prefixes. The native literal-event replay stream owns copied alleles and one projection per event/transcript pair in caller-supplied rings. It drains each occupied path once, with complete event provenance and explicit projection or edit-conflict status. Its mutable storage is bounded by the oldest active genomic window, including younger events retained behind a longer-lived transcript; capacity failures latch instead of dropping paths. Projected equal-length edits use the same differing-island decomposition as independent annotation: unchanged internal MNV bases do not mask or conflict with another carried edit. Source-record count and physical-edit count remain separate, with raw alleles preserved. An out-of-CDS contributor is distinct from a failed coding projection. For a coding transcript, the shared topology classifier can prove that its semantic REF span (both flanks for an insertion) has no coding overlap; such a contributor retains its own outside_cds status without suppressing literal CDS replay. A path containing only these events retains the reference CDS with zero edits and complete provenance. Decoded replay treats mixed coding/noncoding spans, invalid sequence slices and other projection errors as failures. Missing/unphased evidence makes the path incomplete. Literal replay does not predict splice alteration; noncoding transcripts still have no CDS. One shared CDS translator serves independent coding contexts and phased replay. It retains every complete codon’s residue and the first-stop position in one pass; the mutation-path protein is a length-delimited prefix through that stop. Later residues stay in worker storage for coding-context consumption, and later source events remain in contributor provenance. Full and stop-truncated translation do not have separate biological implementations. Coding contexts can borrow a completed replay and its complete raw translations. The same opener serves the context builder and native leaf consumers; it checks sequence/translation extents without applying edits or translating again. Displayed first-stop prefixes and curated reference proteins are not interchangeable with these complete raw peptide views. The model and worker buffers remain immutable while a consumer holds the context. The optional protein HGVS consumer builds frame-closed edit spans and merges normalized peptide operations that touch. Physical coding blocks remain intact: a codon-aligned deletion contributes no residue to the following alternate codon, even when the coding-block partition groups that deletion with a later edit. Predicted suffixes contain one complete supported operation set, without accession or source-identity reassignment. HGVS facts borrow a prepared reference protein independently of raw CDS translation and frame predicates. Reference preparation can change residues and length without inventing physical source edits; reference-only raw paths borrow the exact prepared reference. Protein ends compare through the alternate’s first stop. A reference-only stop-marker loss supplies no extension, and insertions require reference flanks, including the prepared terminal stop when present. Interacting terminal operations retain that stop during peptide clipping. Unrepresentable ends, incomplete sequence and unsupported mechanics have NULL HGVS with an explicit status. Caller-selected operation/text capacities are included in the query workspace and cannot grow during execution. This consumer does not establish complete protein or DNA HGVS compatibility. The leaf’s stop_in_displaced_frame fact intersects the first stop’s three rebuilt CDS bases with physical frame excursions. An excursion starts at a frame-changing edit and ends after the restoring edit’s alternate bases, or continues downstream if unrestored. A zero-base excursion cannot intersect a codon. A stop after frame restoration and a sequence with no stop both return false; unavailable sequence or ordered overlapping replacements return NULL. The shared edit geometry does not modify raw frame flags or imply protein rescue, SO classification, or removal of downstream contributors. Translation uses BioPerl’s amino-acid consensus over every A/C/G/T expansion of N for independent coding predicates and phased replay. Uploaded REF/ALT validity remains a separate check. A resolved residue does not clear the input’s unambiguous fact. The translator validates every base, including trailing partial codons and sequence after a stop. Reference protein preparation uses that same translator, then applies Ensembl’s distinct start-methionine, terminal-stop and curated peptide-edit rules before comparison. It retains internal stops and applies single-residue edits to the complete reference. One worker-owned reference peptide is prepared per closing transcript and serves both native replay and SQL difference materialization. A raw sample without a retained exon-overlapping genotype uses this curated reference peptide. Missing calls retain conditional evidence. Any retained exon-overlapping genotype selects mutation translation, including retained REF lanes and shadowed, unmapped or UTR sources with zero physical edits. Entirely intronic source spans retain provenance but do not select mutation translation or alter literal CDS replay. This also applies to short introns classified as frameshift introns for SO; Haplosaurus admits source records through exon overlap before constructing genotypes. CDS equality and edit count do not select the protein route. Strict decoded replay still withholds sequence for missing calls; pure-reference samples remain implicit. Coding and HGVS share codon-rounded peptide windows with distinct reference and alternate offsets. A physical interaction block opens those operands against the complete materialized path: earlier closed blocks can shift the alternate by whole codons, while the local length request uses this block’s change, not the path’s total change. Window access retains reference peptide edits and terminal partial codons without reapplying edits, inventing a single edit, or truncating the complete context at a stop. Substitution-only blocks feed those operands into the same local predicate interpreter as independent events, even when an earlier closed indel shifted the alternate protein. Those local predicates do not decide whether an earlier stop prevents expression: the complete path retains that separate fact. Indel blocks use the same length-change predicate interpreter as independent events, with actual block geometry validated against their physical edit slice. The complete translation retains its first stop; intersection with frame-displaced bases prevents a DNA-restoring block being called in-frame when translation has already stopped. Reference and alternate nucleotide comparisons use their separate codon offsets. Start/terminal-CDS predicates borrow the selected block’s already rebuilt bases between unchanged reference flanks. This isolates local facts from separate blocks without manufacturing an input record or replaying edits. Complete 5-prime sequence and sufficient or explicitly complete 3-prime sequence remain required when the predicate reads them. Single-record genomic insertion-length reach is retained separately; a CDS span cannot infer that distance across introns. Equal-length or identical strings do not erase physical indels or transient frame changes. Restoring indels that recreate reference CDS retain their physical block and can have synonymous local coding facts; an identity substitution is not a coding change. Empty CDS/protein differences do not remove source provenance. Neither local facts nor a net-zero CDS diff constitute a whole-haplotype consequence set. The whole-context compound-indel substitution shortcut remains forbidden. Leaf-specific facts never mutate shared carrier prefixes or immutable model sequence. Carrier-prefix identity includes per-event called/missing/unphased evidence. An uncertain path cannot share the result of a fully known path merely because their called edits agree. The native stream accepts complete decoded calls for a candidate transcript and uses the same phase reducer as SQL. Its host supplies the complete borrowed phase-set domain for each sample/transcript, including sets first encountered later; homozygous/haploid calls and wholly unphased evidence therefore reach future sets without a second native catalogue. Partial-phase uncertainty is confined to its declared set and unresolved slots. Missing or unresolved paths return explicit incomplete-input status and all contributor evidence, with no CDS/protein, including in the VEP-116 slot profile. This conservative incomplete result does not claim parity with upstream conditional sequences for missing genotypes. One stream cannot mix phase policies.

duckvep_phase_call prepares decoded GT/PS calls through a constant-space native reducer. It observes the complete genotype before assigning slots. Strict assignments respect decoded per-allele phase, with unphased slots resolved only when permutation cannot change the called allele. VCF 4.4 genotype fields define each indicator for the following allele. Thus 0|1/2 and /0|1/2 have unresolved first/third slots; |0|1/2 fixes the first two and leaves one possible assignment for the third. A later pipe does not phase the preceding allele. Homozygous calls and haploid calls apply across every phase set; they must not be put into an isolated NULL-PS bucket. Missing alleles and unresolved heterozygous slots remain explicit and must affect whether a completed sequence is known. The named VEP-116 profile ranks called alleles after omitting missing entries, matching the upstream parser’s input to Haplosaurus; the output still retains those missing input slots with no assigned lane. Separators and PS do not affect compatibility assignments. This is decoded-call interpretation, not emulation of VEP’s raw mixed/prefixed-separator parsing. This helper consumes a declared GT/PS phasing source; PSL/PSO or producer-specific phase identities require a separate explicit adapter.

Typed GT is not a lossless representation of raw VCF spelling. For example, HTSlib decodes 0|1 and |0|1 to the same alleles and phase flags, but VEP-116’s raw parser can give them different Haplosaurus sequences. Exact raw-input compatibility therefore requires retained source GT and source-record allele context; reconstructing text from decoded calls cannot recover it. The finite raw-GT audit retains these collisions and missing-call/ploidy disagreements separately from the certified literal-replay cases. Native raw-record replay uses source record IDs plus REF/ALT ordinals. An undefined file slot is an explicit empty-ALT interpretation of the complete source REF span; its sequence is conditional. Missing REF and omitted-call observations retain source evidence with zero physical edits. A full source span crossing coding/noncoding bases has no single Haplosaurus CDS mapping. Raw replay retains it as source_unmapped with conditional evidence and replays the remaining mapped sources. The shared model layout validator, uncached CDS extent and cached-coordinate agreement are checked; every coding-overlap REF segment is verified through the shared REF validator. A completely mapped source whose ALT contains N, U or lowercase bases within the supported ACGTUN/acgtun alphabet retains source_allele_skipped, conditional evidence and zero physical edits. Haplosaurus’s case-sensitive ACGT mutation gate selects this behavior; the remaining sources still replay, and a retained exonic call selects normal CDS translation even when no edit applies. Uppercase ACGT and an undefined-slot empty ALT apply. Unsupported symbols and dashes remain invalid rather than being stripped or coerced. Intronic/UTR reference sequence is not available from a CDS-only pool. Invalid layout, CDS storage, unsupported alleles or coding REF still make sequence unavailable. SQL/R selects this raw-record interface with input_mode := 'source_records' and phase_policy := 'vep_compat'. The default alt_events input is the decoded-call contract; it cannot emulate lexical distinctions absent from those arrays.

duckvep_haplotypes consumes flat event/transcript/sample calls. DuckDB derives phase domains and materializes sorted input; native event ingestion and candidate projection are separate operations, so an entire event’s cohort is never copied into a first-party call matrix. Output pauses retain the transcript drain cursor and reuse worker scratch.

For input_mode := 'alt_events', the SELECT supplies these named columns; DuckDB casts them before execution:

event_index UBIGINT       seq_region UINTEGER        position UBIGINT
reference VARCHAR        alternate VARCHAR          alt_index UINTEGER
transcript_index UINTEGER sample_index UINTEGER      alleles INTEGER[]
phase_before BOOLEAN[]   phase_set BIGINT

event_index uniquely identifies one source ALT; alt_index is its positive source ordinal, and position is one-based. Region/transcript/sample ordinals belong to the selected model and retained cold relations. NULL allele items mean missing, NULL phase flags mean unavailable, and phase_set is nullable. Candidate selection and source/ALT mapping remain explicit relations. Duplicate calls, inconsistent event/ALT geometry, changing sample/transcript ploidy, invalid GTs and out-of-model candidates fail. Evidence bits are 1 called, 2 missing and 4 unphased. Per-call capacities name active pools, leaves, sequences, genotypes and phase domains; workspace_limit sums DuckVEP-owned buffers, excluding DuckDB input/sort/output memory and HTSlib handle/transport storage. Neither result order nor carrier-list order is a SQL ordering guarantee.

For a single ALT contributor, protein HGVS consumes the shared independent-event VEP-116 consequence and genomic-placement facts. Multiple MNV islands still belong to that one source allele. The completed haplotype CDS, protein, differences and provenance are separate outputs: an independent synonymous HGVS label need not describe the contrast against a curated reference protein. Missing genomic context required for placement returns missing_reference, not an unshifted substitute. Compound and reference-only paths retain their completed-path operation builder; neither route certifies complete phased HGVS.

HGVS scratch is query-owned and allocated before execution. Each query opens its own mutable faidx handle. max_hgvs_reference_bytes bounds caller-buffer retrieval, including line-ending scratch; max_sequence_bases, max_leaf_edits and max_allele_bytes bound the sequence, edit and allele workspaces. The reference reader exposes separate VEP shift and complete-upload/duplication lookup views. Capacity failure is an error; a cache refill overwrites fixed storage only after invalidating the previous window. HTSlib transport allocations remain dependency-owned.

For input_mode := 'source_records', the required columns are:

event_index UBIGINT       seq_region UINTEGER        position UBIGINT
reference VARCHAR        alternates VARCHAR[]       gt VARCHAR
transcript_index UINTEGER sample_index UINTEGER

Here event_index identifies one whole source record, not a decomposed ALT. gt is original VCF spelling, and alternates preserves every ALT in source order. Candidate selection remains explicit. Required cells, ALT items and allele strings cannot be NULL; source allele strings are nonempty. An empty ALT list represents a REF-only record. Duplicate record/transcript/sample calls, inconsistent record geometry/ALT lists, or conflicting GT spellings for the same record/sample across candidates fail before replay. The parser rejects invalid GT grammar and out-of-range allele indices, including slots beyond the two used by the file profile. max_ploidy bounds source ploidy and the two file lanes; PS is ignored, and the decoded-only max_phase_sets workspace is not allocated.

Raw replay treats the supplied records as its source universe, including 0|0 records. Preserve original equal-position file order in event_index; filtering records before preparation can change overlapping-edit replay. The native planner uses the loaded model’s sorted transcript spans to form Haplosaurus input buffers. The first overlapping source span closes over connected transcripts; later records enter by start coordinate without extending the buffer with their REF ends. DuckDB retains buffer ordinals, counts all source records and derives the pinned Set::IntervalTree 0.12 preorder before genotype filtering. Its temporary preparation relations are query-owned; native planning uses constant scratch.

DuckDB expands source REF, ALT and undefined-slot interpretations and sorts them by region, position, record ID, allele ordinal, transcript and sample. Each interpretation/candidate is projected once; calls stream without a native cohort matrix. These interpretations share the named active-event/projection/allele limits, including unused interpretations. Original record/GT relations remain the cold provenance authority.

Source-record contributors append nullable alt_index to their struct: 0 means REF, a positive value is the source ALT ordinal, and NULL means the undefined file-slot interpretation. The latter has empty alternate and deletes the complete source REF span. Evidence bit 8 and sequence_status = 'conditional' distinguish missing/undefined-slot replay, validated source-mapping omissions and skipped raw alleles from a known called sequence; bit 2 additionally retains explicit missing-source evidence. Omitted missing observations have no physical edit; retained REF slots participate in ordered replacement. A source_unmapped or source_allele_skipped contributor performs no replacement and retains its interpreted REF/ALT identity. Other projection failures withhold sequence. Blocks, differences and local coding facts describe the displayed conditional sequence, not proven biology. Overlapping raw records replay complete projected REF/ALT spans in descending original CDS start order; equal starts use source-buffer tree order. REF is validated against the model, while replacement acts on the current sequence and clips removal at its current end. A retained REF slot can overwrite an earlier replacement. An omitted missing observation does not execute a REF replacement. Every operation that changes the current sequence retains its source ID, even when a later operation overwrites it or restores the reference. Sources sharing region, position, REF and the complete ordered ALT list use the last retained source in tree order for that transcript. Other calls remain contributors with projection_status = 'shadowed_duplicate' and do not execute replacements. Retention is determined by the native raw-GT parser across the candidate’s samples. This order is not proof that conflicting calls describe a biological haplotype.

Each known leaf exposes coding_blocks in ascending reference CDS order. The same partitioner used by the independent interaction property groups same-alternate-codon edits and keeps a block open while its frame is displaced. Each block retains its source event IDs in physical-edit order and frame flags, and describes a reference span plus a span of the already rebuilt CDS. Retained bases between edits stay in those spans; pure insertions/deletions have one empty span. This is a composite edit representation, not an alignment or HGVS normalization. At most max_leaf_edits blocks are stored in the initialized workspace, and an unknown or failed sequence has NULL blocks. A parallel event-ID array follows physical edits through splitting, heap sorting and reversal; each block borrows its edit slice. An uploaded MNV may contribute several islands within or across blocks, so repeated IDs are retained, not deduplicated. SQL event_indices replaces the count-only block field; its length is the physical edit count. The ID array shares max_leaf_edits and is included in the workspace byte limit. Complete raw contributor provenance remains on the leaf.

For an overlapping raw-record leaf, blocks are disjoint net reference/alternate components rather than differing islands. Their ID lists retain applied full-span operations in ascending reference CDS order, including overwritten operations; edit_count counts those operations. Net-zero components retain their provenance. Block flags describe net component length changes; leaf flags retain nominal source replacement length changes. nominal_length_diff is the signed sum of projected replacement ALT-minus-REF lengths before clipping at the current CDS end. It is zero for reference-only replay and NULL when the leaf has no CDS. Known and conditional paths retain this fact even when the rebuilt CDS length change differs. Local SO is NULL with coding_status='unsupported_ordered_replacements', and stop_in_displaced_frame is NULL because an overlapping operation history does not supply the disjoint physical edits required by those consumers. CDS, protein, aligned differences and after-first-stop positions still describe the literal replay. Disjoint raw records use the differing-island contract.

Each block’s local_consequence_mask evaluates the shared generated SO rules over that block’s physical coding delta in the completed haplotype. It does not OR independent event labels or apply an arbitrary contributor’s uploaded-feature class gates. Decode it with duckvep_so_terms(). coding_status is ok, unsupported, unsupported_ordered_replacements, missing_transcript_tail, missing_transcript_flank, or invalid_argument; only ok has a non-NULL mask. Known zero is distinct from unknown. Substitution predicates consume consensus peptides and retain an X-bearing peptide’s independent coding-unknown flag. Unsupported length-changing contexts remain explicit even when consensus sequence replay is available. after_first_stop means the block’s first alternate codon is strictly after the first translated stop codon; it is false if no stop exists. A block starting before the stop but spanning it is not marked. Later blocks retain their local facts and provenance: this positional fact does not assert biological expression or rescue. The coding context borrows the complete alternate translation and model overlay. One native translation pass per closing transcript prepares two worker-owned reference views: uncurated consensus coding operands and the curated reference used for protein differences. SQL borrows both views. Neither model mutation nor per-leaf replay, translation or allocation is required by this consumer.

cds_differences is a separate alignment view, not a change to physical edit identity. Indel-bearing leaves use the pinned VEP-116 pure-Perl NW score and traceback tie order; substitution-only leaves compare corresponding positions. Differing columns join only when both sides retain the same gap/non-gap type. Each run borrows ungapped reference/alternate spans and names zero-based positions on both sequences and on the alignment. No HGVS normalization or contributor reassignment is inferred from repeat-associated gap placement. Unknown sequences have NULL differences. One worker-local reference view per closing transcript uses replay’s uppercase DNA spelling; model bytes remain immutable and letter case alone cannot introduce differences. A feasible alignment supplies a cost upper bound U under the equivalent nonnegative cost (substitution 4, gap 3); every optimum lies within |i-j| <= floor(U/3). Two worker-owned score rows and a caller-bounded traceback band therefore preserve exact global tie placement without allocating in execution. max_alignment_cells and max_leaf_differences are independent per-call limits within workspace_limit; exceeding either is an error, never approximate output. protein_differences uses the same alignment and span contract in amino-acid coordinates. The reference peptide is prepared once per closing transcript in worker storage of at most max_sequence_bases/3 + 2 bytes: a terminal curated edit can restore a removed residue, and Haplosaurus can append another stop. CDS and protein axes sequentially reuse the same traceback and descriptor arrays; DuckDB copies each list before the next axis resets those arrays. Limits apply separately to each axis. A reference CDS shorter than a complete codon has an unavailable protein comparison (NULL), not a known empty reference. Known equal proteins produce an empty list. The exact raw-suffix stop convention is recorded in the compatibility errata. The registry owns one valid retained query connection: its extension-load database handle must not be retained. Only preparation/materialization uses that connection. A busy slot returns an error, including recursive preparation, while completed scan results and native state have independent ownership. Caller TEMP objects/uncommitted writes are not visible.

This public sequence-mechanics surface is not complete phased annotation. Local block masks and sequence/indel flags are not whole-haplotype SO or protein HGVS, and literal replay does not yet compose typed structural events. Existing executable Haplosaurus comparisons exercise native replay and the public SQL surface for their declared phased-sequence scope; they do not certify compound SO/HGVS. The phased replay benchmark separates the initialized native stream from public SQL sort/materialization and records sparse pool occupancy, native workspace bytes and whole-process RSS. Its shared synthetic cohort measures literal replay, not future combined annotation or real-population performance; those paths require renewed evidence when implemented.

The stream must preserve the original record/ALT identity, decoded allele indexes, ploidy, phasing flag, and PS/PID-like phase-set provenance. The same called local haplotype may arrive as one MNV, several SNVs, or overlapping SNV/indel records; those representations must yield the same alternate CDS and peptide when they encode the same phased edits. Independent annotation of each row is not an acceptable substitute. The existing phased-SNV-versus-MNV property covers one important subset; mixed replacement and indel representations still require generated equivalence and executable csq/Haplosaurus corpora.

Long-read callsets make that contract unavoidable. Clair3 or DeepVariant may emit phased SNVs/MNPs while Sniffles or cuteSV emits a structural record for the same sample and local haplotype; the biological edit set can cross those record classes. HBA1/HBA2-like paralogous loci add a separate ambiguity: an aligner/caller may report one of several near-identical placements. DuckVEP must retain call, alignment, assembly-path, and phase provenance and must not manufacture certainty by merging records solely because their nominal coordinates are close. Haplotype-resolved HPRC assemblies can be loaded as separately receipted models or explicit paths; comparison to GRCh37/38 remains a mapping relation, not a chromosome-name alias.

The annotation-model receipt is not enough to reproduce a long-read result. A callset receipt must also identify the read chemistry, basecaller and model, alignment reference and aligner, small-variant/SV caller versions and options, phasing method, and any callable or confidence masks. Agreeing GT and phase-set fields do not by themselves make overlapping records compatible. The phased executor must prove that their reference and alternate sequences form one consistent local haplotype or return an explicit edit-conflict result while retaining every source record.

An assembled HPRC haplotype is useful truth evidence, but it is not automatically a DuckVEP model. Path-coordinate annotation additionally requires transcripts projected or annotated on that exact path, mapping confidence, and locus-level assembly QC; a nominally haplotype-resolved assembly can still carry a flagged collapse or misassembly. Incremental annotation therefore invalidates work by dependency: an independent new allele can be annotated alone, while a changed phase block, caller interpretation, transcript projection, or model receipt requires recomputing the affected transcript haplotypes.

Sorted input bounds lifetime: a transcript’s phased state can be finalized once the stream passes its end. Reference paths stay implicit; non-reference paths should share compact edit prefixes across samples and translate each distinct leaf once. GT/PS decoding, arbitrary ploidy, transcript-close flushing, and VEP Haplosaurus comparison are tracked at https://github.com/RGenomicsETL/duckhts/issues/92.

DuckDB may decode genotypes, derive explicit sample/phase-set/haplotype columns, and sort the input relation by genomic coordinate before streaming it. That SQL grouping is not the same as biological edit interaction: the C executor still owns transcript lifetime, same-codon interactions, an open displaced frame, frame restoration, and the complete list of contributing variants. Benchmarks must include both an already-sorted stream and the same workload with DuckDB sorting included.

Structural events

Small edits and structural events share overlap, projection, provenance, and output, but not one lossy event shape. Deletion, duplication, inversion, CNV, breakend pairs, inserted sequence, and repeat changes retain their typed geometry. A breakend is two loci, not a wide interval or symbolic point.

The public event relation accepts exact single-locus events plus a typed DEL, DUP, tandem-DUP, tandem-repeat (STR), INV, INS, CNV, or unknown operation and an explicit loss/neutral/gain/unknown copy direction. Span operations use one-based inclusive start/end coordinates. An insertion uses start = end = P for the interbase site after reference base P; preparing symbolic VCF therefore removes the left anchor, maps a span to start = POS + 1, end = INFO/END, and maps an insertion to P = POS. The dispatcher rejects contradictory operation/direction pairs. A BND row supplies both mate coordinates and never pretends that its two loci are one single-locus span.

VEP expands a bounded <CNV:TR> from RN plus RUS/RUC or RB into literal alleles when the result fits its configured structural-size limit. Such an event is then an ordinary VariationFeature and enters DuckVEP’s small-variant path. An oversized or unexpanded repeat remains a structural tandem_repeat. DuckVEP’s STR type preserves that identity for provenance and later HGVS while reproducing VEP’s tandem-duplication gain/insertion predicates. Raw repeat units and counts remain columns in the surrounding relation; the consequence kernel consumes the prepared literal allele or exact structural span, not parser-specific INFO strings.

Typed repeat preparation must distinguish exact sequence from repeat-summary evidence. VCF 4.5 section 5.7 permits a <CNV:TR> summary to omit SNVs and indels present in the corresponding literal phased allele. Equal repeat units, counts, and lengths therefore do not establish sequence identity: (CAG)11 and (CAG)5(CAT)(CAG)5 can share summary metadata but encode different proteins. Sequence replay consumes the explicit literal allele or a caller-qualified exact repeat description, not a summary silently promoted to an exact sequence. VEP-derived expansion is a named compatibility interpretation, not evidence that the sample has a perfect repeat. The preparation contract must retain ordered repeat components, source ALT ordinals, nullable counts and lengths, confidence intervals, and exactness separately; RUC is a floating-point field and must not be silently narrowed to an integer count.

duckvep_repeat_alleles prepares reference and alternate ordered (unit, count) descriptions as one event fact. A required sequence_exact assertion separates exact descriptions from summaries. Missing data and fractional counts withhold both sequences with an explicit status. Complete integral counts expand only after each complete allele fits max_allele_bases; the result includes both lengths, their signed difference and the gain/loss/neutral base-length direction. This is not a copy-number inference. The limit belongs to the SQL call, not to VEP’s per-component parser limit. Empty lists describe empty alleles; a missing list means unavailable sequence. Input component order, case and IUPAC codes are preserved. The caller still supplies reference validation, genomic coordinates and source identity before annotation; summary data cannot enter sequence replay merely because its nominal counts happen to be integers.

The paired-breakend native lane accepts the local and mate regions and raw one-based VCF positions from that same public event row. It queries the resident cgranges transcript and regulatory/motif indexes around both loci, merge and deduplicate each object class, and call shared C evaluators with explicit variant/object pairs. This is intentionally separate from the sorted single-locus sweep: neither a wide span nor two independent endpoint calls can reproduce VEP 116.

The exact VEP state is asymmetric. BaseVCF4::get_start moves the local BND POS to POS + 1; StructuralVariationFeature::_parse_breakends retains the mate coordinate. Ordinary intron, exon, UTR, splice, and coding predicates use only that shifted local feature. Candidate discovery also creates an overlap allele for the mate, but the mate changes transcript consequences only through the exceptional feature_truncation predicate. VEP may therefore emit two internal overlap-allele rows for one BND and transcript. DuckVEP returns their consequence-set union once. A transcript reached only through the mate has a zero region mask and NULL rich region because no local topology exists.

VEP’s RegFeat lane also evaluates both points. A RegulatoryFeature or MotifFeature hit by the shifted local point, the verbatim mate point, or both produces one DuckVEP object row. The result is asymmetric: a shifted-local hit retains regulatory_region_variant or TF_binding_site_variant; a mate-only hit takes VEP’s generic HIGH-impact feature_truncation chromosome-breakpoint branch without requiring deletion or copy-number loss. Once the mate has discovered an object exactly, VEP also attaches a shifted local point on the same contig when it is outside but within the fixed 5000-base structural-feature admission distance. That local allele falls back to intergenic_variant, so the object-level result is feature_truncation&intergenic_variant. A close point does not discover an object by itself. If both points hit one object exactly, the local base term wins. VEP may materialize duplicate identical rows or distinct allele-level rows for one object; DuckVEP’s public contract is their consequence-set union once per (event, object).

The surrounding DuckDB relation retains event identity, bracket orientation, raw ALT, and provenance for HGVS, fusion, and round-trip consumers. Orientation does not change the transcript consequence set, so it is not an ignored kernel argument. The C lane preserves VEP’s fixed 5 kb endpoint admission cap in addition to the configured directional upstream/downstream distances. These are independent controls despite sharing the same default number. Raising the caller distance to 10 kb may widen upstream/downstream transcript terms, but cannot turn an endpoint 5,001 bases from the same transcript or interval feature into a StructuralVariationOverlapAllele. Interval-feature candidate discovery remains exact; the fixed-cap endpoint is attached only after the mate has discovered that same object. This does not clamp ordinary transcript predicates: an overlap allele created by the mate still reads the shifted local feature, so a 10 kb caller window can emit an upstream/downstream term for a local point beyond the fixed allele-admission cap. Pure-C, SQL, and R regressions pin 5,000 versus 5,001 bases under both a 10 kb caller distance and a zero caller distance. In the latter case an admitted local transcript allele has no directional predicate and falls back to intergenic_variant, which is unioned with a mate-derived feature_truncation. Randomized sweep scenes include zero, 1, 50, 100, 4,999, 5,000, 5,001, 10,000, and 65,535-base windows rather than treating 5 kb as a maximum allocation size. A seeded executable differential covering chromosomes 1, 2, 7, 21, and X, all four bracket orientations, same- and cross-chromosome pairs, and transcript/exon/intron/CDS/flank endpoint states matched all 91,428 transcript pairs from 1,004 generated events.

That differential isolates every BND in VEP with buffer_size=1. VEP 116 otherwise puts mate coordinates from several records into one chromosome-blind interval tree, so a neighboring BND can change the executable oracle’s transcript set. One VEP process is retained for cache reuse; only the semantic event buffer is isolated. Chromosomes stay contiguous and positions increase within each chromosome in the generated VCF.

The public SO mask binds all 41 terms registered by VEP 116. Six regulatory-region and transcription-factor-binding-site terms are produced by a separate interval-feature evaluator over the same event geometry. Their five hot columns live in a separate resident SoA, not in fake transcript rows; cold funcgen metadata remains an ordinary DuckDB relation. The small-variant and exact single-locus structural adapters advance the transcript and feature sweeps together and emit typed rows distinguished by overlap_object. Regulatory or motif output therefore does not require a SQL range join, and it preserves the same resumable output cursor as transcript output. sequence_variant, the forty-first registry term, has no VEP overlap predicate and is retained as registry metadata rather than emitted to hide an incomplete model.

Raw BND ALT parsing and repeat-unit/count expansion remain outside the consequence kernel. VEP 116 stores CIPOS/CIEND and structural inserted-sequence payloads, but its registered 41-term consequence predicates use nominal POS/END and its structural inframe_insertion branch explicitly has no inserted-sequence implementation. Those fields therefore remain required provenance and future HGVS/round-trip inputs, not missing consequence facts. The BND statistical differential is part of the executable-VEP harness; broader chromosome, species, and real fusion corpora remain continuing evidence rather than a second implementation. Remaining structural follow-up is tracked at https://github.com/RGenomicsETL/duckhts/issues/98.

Callers such as Sniffles and cuteSV may provide CIPOS/CIEND, mate identity, inserted sequence, copy number, and several records for one event. For strict VEP-116 consequence parity, an upstream relation may annotate the nominal POS/END while retaining every confidence interval and payload beside the result. It must not relabel the nominal span as experimentally exact, discard those fields, or reuse it as HGVS/fusion geometry. Likewise, ambiguous placement between paralogous loci is source evidence for downstream SQL; choosing one locus is not a consequence-kernel inference.

HGVS

HGVS is a consumer of the projected edit, not a formatter over consequence names. The current event keeps the uploaded span, VEP feature span, minimized REF/ALT edit, insertion boundary, and anchor side separately. The transcript model keeps complete spliced pre-CDS/CDS/post-CDS sequence, and the sequence layer can apply one edit or a grouped edit set and compare the resulting peptide. These semantic facts are shared by independent and phased paths.

The implemented internal layer makes that ownership explicit:

query(duckvep_annotate_sql(..., struct_pack(hgvs := true))) exposes those mechanics without changing the public schema: the first 16 fields are the compact consequence row, followed by nullable transcript/protein suffixes, transcript-direction shift, and separate structured status/reason fields. Structural and breakend rows currently leave those HGVS fields NULL. Stable versioned transcript and translation identifiers remain ordinary prepared-model columns and are joined by transcript ordinal rather than copied into every hot result row. A model may bind an existing indexed FASTA through an exact sequence-region ordinal/name/length relation; the loader validates the index without creating or modifying it, while each annotation worker owns its faidx handle and reusable reference-window storage. One bounded fetch supplies distinct borrowed views for shifting and lookup. The adapter performs one candidate sweep and renders into DuckDB-owned output vectors. One worker-owned scratch buffer handles transcript and protein strings in a single render pass; it grows and retries only when a result does not fit.

Consequence predicate flags are reusable evidence, not a closed-world serialization of the later HGVS replay. Positive frameshift evidence can complete a length-changing delta, and absence of frameshift is conclusive for a length-preserving CDS edit. A length-changing splice-overlapping edit without that positive flag must run the complete delta evaluator; otherwise VEP-compatible frameshifts can be misrendered as premature stops or delins. An original partial-codon, retained-stop or available leading-stop peptide operand is explicit negative evidence, captured before HGVS shifting reuses worker scratch. The shared peptide-window guard preserves these exclusions without treating an absent frameshift flag as conclusive.

This implemented surface is not yet a full VEP-HGVS compatibility claim. It covers independent literal small variants and returns explicit unresolved reasons when reference, projection, transcript flank, tail, or protein facts are unavailable. Fixed executable VEP witnesses cover position-one right anchoring, endpoint-clipped deletions/delins, terminal insertion states, both transcript strands, and VEP’s short alternate-CDS trimming bug. A strict chromosome-21 ClinVar run compared 56,998 transcript pairs: 20,782 matched both HGVSc and HGVSp, 24,089 matched HGVSc with HGVSp absent on both sides, and 12,127 had both strings absent, with zero discordant, unresolved, missing, or extra HGVS rows. That is evidence for the exercised GRCh38 independent-event distribution, not for untested species, assemblies, transcript-source-specific sequence/projection states, or other HGVS classes.

The execution split is:

The hot transcript sweep therefore emits or retains numeric projected-edit facts; it does not allocate HGVS strings. Rendering is late, after filtering. The phased executor consumes the same prepared CDS edit set rather than introducing a second trimming or projection authority. External VEP-116 differentials cover independent events on Ensembl/GRCh38; additional evidence must cover transcript catalogs with sequence corrections, other assemblies/species, all shift modes, exact structural events, and compound edits. Apply-then-diff sequence equivalence remains an independent property oracle. Exact structural HGVS can later consume typed exact events. BND HGVS additionally needs the paired relation’s mate and orientation facts; imprecise structural events remain unsupported until their confidence geometry is represented.

The pinned ferro-hgvs v0.9.0 source is an independent HGVS-spec oracle and a useful model for structured fuzzing, large ClinVar corpora, and hermetic reference fixtures. It is not a DuckHTS dependency and does not replace the VEP 116 executable as the compatibility authority: canonical HGVS, Mutalyzer/biocommons behavior, and VEP’s historical output can legitimately disagree. The intended differential therefore has three independent observations: exact VEP-116 output, ferro-hgvs parsing/normalization where its contract applies, and DuckVEP’s apply-then-diff sequence replay.

Supplementary annotations

Exact sources such as dbSNP, ClinVar, and gnomAD should remain sorted numeric-key streams; range sources should remain interval streams. They share prepared variant keys and final output, not a universal physical format. DuckDB/Parquet pruning, parallel reads, memory management, and spill are the first implementation. A custom cache requires a measured workload that those facilities do not meet.

Variant-level exact and positional joins happen before transcript expansion. The compact consequence row already carries transcript and gene ordinals, so affected rows are reduced to distinct gene ordinals before gene resources are joined and then attached back. Ensembl RegulatoryFeature and MotifFeature are core VEP inputs and run in the C consequence kernel. Protein domains and genuinely supplementary interval sources use range joins or sorted interval streams, followed by the smallest source-specific predicate needed for a joined pair. Strings and JSON remain the final projection.

The serving contract is manifest-driven generated SQL, analogous to hts_union_query(...), rather than one C table function per database. A provider receipt must declare at least provider/release identity, assembly, normalization contract, source digest, key kind (variant, position, interval, or gene), partition layout, and hot payload columns. Assembly and normalization are semantic inputs: a source such as a TSV without an assembly declaration may be receipted and inspected, but must not be keyed or joined. Provider-specific licences and cold presentation metadata remain in the receipt or side relations.

For exact human providers, preparation writes two hot lanes. Reversible normalized alleles use only UBIGINT vk plus narrow payload columns. Hashed/nonreversible alleles retain the normalized contig/position/reference/alternate tuple and refine every VariantKey match on that tuple. Non-human providers use a release-scoped numeric contig ordinal plus exact alleles rather than pretending the official human VariantKey contig space is universal. Position, interval, and gene providers use their own typed keys; they do not inherit exact variant semantics merely because all are exposed by one manifest.

FastVEP’s .osa/.osa2, .osi, and .oga files encode exact-allele, interval, and gene provider scopes respectively. DuckHTS closes over those semantics with ordinary typed relations rather than reading or creating another cache format. Sorted/partitioned Parquet holds the hot projections; provider receipts and cold metadata remain relations; tabular or JSON output is rendered after selection. VCF round-trip output remains a separate writer and FORMAT-remapping contract.

The generated annotation query selects touched chromosome or coordinate-tile Parquet shards and emits provider identity/release with every result. DuckDB may cost several narrow exact equality joins in one plan; key count alone is not a reason to serialize them. Every provider must declare whether its hot key is unique, one-to-many, or requires a finer key, because an accidental duplicate multiplies downstream rows even when the hash join itself is fast. Provider matches are reduced or nested before transcript expansion.

At population scale, provider joins run in RSS-bounded groups. The query relation or smaller filtered key set is the hash-build side while large provider relations stream with Parquet min/max and dynamic Bloom filters where the planner can derive them. Several chromosome tasks may run concurrently, but their mutable join state is per task even when immutable Parquet pages are shared. The scheduler therefore bounds concurrency from measured peak RSS. It must not load a dozen full providers or build a whole-dbSNP cgranges index merely to avoid SQL.

cgranges remains the repeated-query choice for compact interval tracks when its measured resident cost is acceptable. The current public registry also owns output coordinates, label validity, labels, and one chromosome string per interval; this is materially larger than cgranges’ 16-byte core interval. Count-only payload elision and chromosome interning are valid future optimizations. Exact dbSNP, TOPMed, REVEL, AlphaMissense, ClinVar, and frequency sources stay sorted/partitioned Parquet streams. The measured storage, serving RSS, interval-index density, real multi-provider plan, join-count stress, and twelve-provider logical-time envelope live in benchmarks/benchmark_variantkey_join_overlap.md.

Supplementary source plumbing does not belong as ignored arguments in the consequence-kernel API.

Validation and performance

The registered VEP-116 consequence engine is closed for its declared contract: supported independent literal small variants, typed exact single-locus structural events, paired breakends, core RegulatoryFeature/MotifFeature objects, and the assemblies/species/model receipts represented by the checked differentials below. “Closed” means new work must preserve this implementation as proven infrastructure and extend it through an explicit new input/output contract. It does not claim phased or compound consequences, raw bounded <CNV:TR> expansion, producer-specific symbolic/BND parsing, supplementary databases, or an untested model receipt. Strict VEP-116 parity here means nominal POS/END consequences while preserving CIPOS/CIEND as relational uncertainty metadata; it does not invent a confidence-aware consequence policy that VEP itself lacks. Held-out fuzzing and executable-VEP comparisons remain permanent regression gates rather than an indefinitely open prerequisite for every release.

Finite-quotient conformance certificate

Statistical differential testing remains a permanent release gate because it reaches rare combinations that fixed witnesses and ordinary corpora do not. It can be strengthened, but not replaced, by a machine-checkable finite quotient of the declared VEP-116 state machine. The formalization starts from the checked test/duckvep/conformance/data/state_machine_transitions.tsv transition relation. Each row names the input fact authorities and output state, every comparison that the transition observes, the implementing source, the pinned VEP authority, and its current proof class. The generated check rejects unproduced inputs, disconnected outputs, missing implementation files, absent fixed tests or randomized property names, a proof class without its required executable campaign receipt, or a claim that an unimplemented transition has evidence. state_machine_campaigns.tsv resolves each campaign identifier to an exact source-revision/artifact/corpus/model/oracle key. The checker parses the selected rows and rejects a nonzero discordance/missing/extra counter, a non-exact all-row, or an HGVS comparison other than match/both-absent; a corpus-name substring cannot certify a failed campaign. Historical campaign evidence remains source-pinned. Before a release, make duckvep-state-current-check additionally requires every named executable campaign to have been run on an ancestor of the checked-out HEAD with no intervening change to the complete extension source/vendor closure, VCF/FASTA readers, DuckVEP implementation, build/catalog inputs and pinned build submodule, property harness, upstream semantic mirrors, or executable conformance inputs. Compiler-like untracked inputs and a dirty build submodule are also fatal. Evidence-only commits do not invalidate the campaign they record. state_machine_outcomes.tsv names fail-closed input/resource errors and unresolved/not-applicable result statuses, with their implementations and fixed witnesses; transitions declare which terminal outcomes they can produce. A + in input_states denotes fact authorities merged by the transition. The paths column names concrete execution lanes; the checker proves that every implemented transition on seven executable lanes has reachable inputs from raw_allele and that each lane reaches its declared terminal state. The combined haplotype-classification lane is reported separately as one structurally connected planned path; its not_implemented transition cannot count as executable reachability. A lane that does not request HGVS or that does not describe a structural/interval event leaves the corresponding optional fact state absent rather than manufacturing a value.

For one annotation object, define the execution state as

S = (uploaded allele, canonical event, candidate object, topology,
     projected edit, coding context, predicate cache, consequence set,
     NMD result, HGVS facts, emission cursor)

and each production stage as a typed partial transition T_i : S_i -> S_(i+1) + Error. An error is an observable terminal state: unsupported, unresolved, missing emission, and extra emission are not erased before comparison. Transcript and regulatory-object paths share event preparation and SO selection but have distinct candidate/topology transitions. The haplotype rows are in the same relation and are explicitly marked mechanics_only or not_implemented; independent-event evidence cannot silently certify them.

This is counterexample-guided abstraction refinement, not an attempt to replace empirical testing with a small hand-picked table. The finite abstraction proposes equivalence classes. Held-out statistical alleles, generated exact SV/BND events, complete ClinVar, and HPRC long/pangenome alleles then search for a counterexample. The SV and HPRC campaigns are therefore state-space exploration inputs: long alleles, graph-derived representations, and rare topology combinations are valuable precisely when they reach states absent from ordinary short-variant corpora. Every counterexample is preserved in pair-level evidence and must either refine an observed dimension, add a satisfiability constraint, expose a wrong transition, or be declared outside the input contract. A successful sampled campaign never proves an unobserved equivalence class. After a sampled counterexample identifies a reachable rare class, the generator must give that class a dedicated stratum and the coverage manifest must require it. Repeatedly hoping that the class appears under an unrelated broad distribution is not an acceptance gate.

Ensembl’s own tests remain a separate source of semantic witnesses. Exact sources are mirrored under test/duckvep/upstream/ using repository/ref-preserving paths, and sources.tsv binds every mirrored file to its repository, exact commit, upstream path, SHA-256, and validation role. The offline generated check rejects drift or an unreceipted test. self_test_receipts.tsv separately binds a successful upstream execution to the exact VEP source commit, checksum-pinned module distribution, Perl version, assertion denominator, and explicit skip counts. The release/116 VEP suite checks the Perl oracle environment and useful parser/regulatory/SV/Haplosaurus paths, while the release/116 ensembl-variation variation_effect.t and hgvs_parser.t files contain the denser consequence and HGVS cases (variation_effect.t contains 189 targeted consequence cases). Extracted DuckVEP cases retain repository commit, file hash, source line/test name, model and expected result. These suites are modest and do not enumerate this document’s whole state relation, so they cannot replace the finite rule proof, named compatibility witnesses, or the statistical and corpus campaigns. The legacy monolithic VEP tests in the ensembl-tools Git history are smaller than the current suite. Ensembl’s CVS history is in Git; VEP’s own “subversion” label denotes the point-release component, not an Apache Subversion repository.

The certificate has three separate proof obligations:

  1. The generated fact-to-SO program is finite. Its generator currently requires one consequence-bearing predicate bit per rule, an empty forbidden mask, a unique binding, and three suppression groups. Empty input, every singleton, every cross-group pair, and the all-bits case therefore form a complete algebraic basis for every possible predicate mask. The property suite must check both ordinary and structural evaluators against the reference interpreter over that basis, and generation must continue to fail if a future rule invalidates the argument.
  2. Geometry and transcript state are quotiented by every comparison that production code observes: event shape; strand; first/internal/penultimate/last exon; relative position at and immediately around transcript, exon, intron, CDS, start/stop, splice, and NMD thresholds; length-change sign and modulo three; phase; coding/biotype state; and exact structural containment/truncation relations. A deterministic relational generator will reject unsatisfiable combinations and retain one minimal concrete VEP witness for every reachable branch signature. Zero-observation classes fail the certificate.
  3. Sequence and HGVS state are quotiented only where an equivalence argument exists: equality/prefix/suffix/internal difference, copied-source availability, stop placement, repeat period, typed peptide edit, and VEP shift distances at 0, 1, 999, 1000, and the first disallowed step. Arbitrary sequence length and repeated-byte loops still require reference-interpreter properties, sanitizer runs, long-allele distributions, or a bounded/symbolic C proof with explicit loop invariants. A short representative list is not evidence for all strings.

The resulting certificate may claim exhaustive agreement only for the declared abstract domain and must publish its dimension tables, satisfiability constraints, reachable class count, witness mapping, and executable-VEP result. Full ClinVar HGVS, held-out generated seeds, long HPRC alleles, structural/BND campaigns, multiple assemblies, and non-human species remain independent gates: they test whether the abstraction omitted a real source state or whether model construction supplied the wrong facts. The fact-to-SO transition already has its finite-basis proof in generated_effect_lookup_matches_rule_interpreter; geometry, sequence, HGVS, SV, and haplotype rows retain their narrower proof status until their quotient and reachability certificate exist.

HGVS and haplotype performance use cumulative lanes over identical prepared input:

Lane Required work
consequence compact candidate discovery, consequence facts, compact rows
consequence + HGVS facts the compact lane plus transcript-edit and typed DNA/protein facts, without strings
consequence + rendered HGVS the preceding lane plus accession lookup and rendered bytes
bcftools local-CSQ projection independent-event consequences plus its declared local output fields
phased edit sets genotype grouping, active transcript state, unique alternate leaves, combined translation, and attributed rows

Every rendered record states input records and ALT alleles, candidate pairs, output rows, rendered bytes, thread count and core pinning, model/corpus identity, and source revision. Sorting-plus-execution is recorded separately from execution over already ordered input. A single independent edit is never relabelled as a haplotype benchmark.

Independent-event parallel execution uses disjoint internally ordered input branches. Each concurrent scalar callback borrows one exclusive workspace from the resident model’s pool; mutable sweep state, result builders, faidx handles, and reference windows are not shared, while the immutable model is. Stable active-set compaction preserves the ordered contents of each variant’s result list, so its per-input annotation_index is stable across DuckDB vector and partition starts. The benchmark’s commutative fingerprints prove equality of the (input_variant_index, annotation_index, public row) multiset, not global emission order. Consumers needing canonical order use ORDER BY input_variant_index, annotation_index after the parallel union. Haplotype execution cannot use arbitrary coordinate cuts because a cut may bisect a transcript/phase-set state.

The generated reports are the numeric authority. At the latest tested ancestor for each corpus, the declared GRCh38 dbSNP, GIAB, ClinVar coding, and ClinVar cross-chromosome samples, GRCh37, and P. falciparum are exact against their indexed caches. The two GRCh38 ClinVar samples contain 604,233 transcript pairs together, with no unresolved, extra, missing, or discordant row. The separate NMD-plugin differential is exact for all 68,554 eligible transcript pairs, including the 29,416 states both implementations leave unresolved.

On the final 644,427-transcript GRCh38 model, the nearest recorded one-core compact baselines preceding source e25c151 are about 841,000 input alleles/s for the transcript-only GIAB topology sample, 801,000/s for its regulation/motif form, 332,000/s for repeated coding SNVs, 112,000/s for repeated coding non-SNVs, 135,000/s for the repeated mixed coding set, and 74,000 semantic BND events/s. They include stable-API list materialization and aggregation but exclude model load and input staging. These recorded rates are nearest baselines, not current-head no-regression evidence. The annotation-dense matrix below is a measurement at its recorded source revision, not a general claim about later builds. The one-million-input-allele target is not met on every final-model workload; output-row rates are a second denominator, not a substitute for the site-rate target. Exact revisions, checksums, row counts, resource measurements, and conditions live in the generated conformance report and throughput report.

The annotation-dense ClinVar matrix is the stricter execution check: at the default 5,000-base transcript distance it expands 517,097 alleles into 26,518,787 rows at 245,535 input alleles/s on one pinned P core and 818,191 input alleles/s on four pinned P cores with four explicit ordered partitions. The full public-row fingerprints match. Those commutative fingerprints establish the same per-input indexed row multiset, not global emission order. Source ownership makes the resident model immutable and shared; GNU time -v observed a 7.03 MiB process peak-RSS increase, not a model-sized step, in this one/four-branch pair. That process measurement is not allocation attribution. At 50,000 bases the same inputs expand to 88,784,213 rows and the observed four-branch RSS premium remains 8.16 MiB.

The conformance report records the closed consequence campaign, while https://github.com/RGenomicsETL/duckhts/issues/95 tracks remaining throughput work. Phased execution remains a distinct semantic vertical at https://github.com/RGenomicsETL/duckhts/issues/92.