This document describes the implemented DuckVEP architecture. Code and tests define executable behavior; open compatibility scope and evidence are linked in the relevant sections.
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.
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:
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:
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.
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.
docs/functions.md is the public SQL contract.src/duckvep_ensembl.c registers the relation-to-model SQL described below.src/duckvep_sql.c registers the public event-relation dispatcher and generated
SO metadata relation.src/duckvep_model.c validates, owns, publishes, pins, and drops named models.src/duckvep_annotate.c adapts DuckDB vectors to the pure C batch interface
and materializes results.src/core/ holds the row logic both hosts share; host_v2/ is the C API v2 host.src/kernel/include/duckvep_kernel.h is the host-neutral kernel ABI; the sibling
kernel sources own traversal, topology, projection, coding edits, phased edits,
structural geometry, and fact-to-SO evaluation.test/sql/duckvep_ensembl.test exercises model preparation and publication;
test/sql/duckvep_annotate.test exercises the SQL adapter; test/duckvep/property/ and
test/duckvep/conformance/ own pure-C properties and VEP differentials.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.
query(duckvep_ensembl_regions_sql(...)) and query(duckvep_ensembl_transcripts_sql(...)) read these Ensembl
core relations by name from the supplied schema:
coord_system and seq_region identify the requested assembly and its regions;gene, transcript, exon, and exon_transcript define transcript topology;translation defines coding start and end within ranked exons; andattrib_type, seq_region_attrib, transcript_attrib, and translation_attrib supply
codon tables, consequence-relevant flags, and exceptional sequence edits.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.
query(duckvep_ensembl_regions_sql(...)):
end - start bases;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.
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 or 0 adds no prefix, while a
positive phase adds one or two leading N bases; andcodon_table attribute, defaulting to table 1 exactly as
VEP does, and rejects conflicting, malformed, or unsupported table IDs;initial_met, _selenocysteine, amino_acid_sub, and
_stop_codon_rt Translation SeqEdits into a sparse reference-peptide relation; andmiRNA 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.
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.
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):
p of a lifted region of length L executes at p + B, where B is a
multiple of L no smaller than the widest reference window (the 1000-base HGVS shift limit,
two maximal alleles and slack), so a window around any event lies inside the lifted interval;p + B + k*L for k in {-1, 0, +1},
which cover every relative placement of two spans shorter than L; a wrapped span [s, e]
becomes [s, e + L], and the CDS and mature-miRNA endpoints of a wrapped transcript move by
L exactly when they lie before the transcript start. Exon rank order, cDNA coordinates,
phases, peptide edits and the CDS and flank sequence pools are copied or shared unchanged, so
continuity comes from exon rank on both strands, never from sorting genomic starts;p + B. A reference allele that runs past L wraps onto base 1 (it must be
shorter than L). A lifted reference window [a, b] is fetched as [start, L], whole laps,
and [1, end] into the worker’s existing scratch through its existing faidx_t; no second
handle or buffer exists;duckvep_lift_resolve maps them to the source
ordinals and keeps one row per event/object pair. The image with the smallest gap to
the event wins; ties go to the larger overlap, then to the image lying after the event, then to
the unshifted image. Only relative geometry decides, so the choice does not depend on where
the origin is. Output contains no genomic coordinate other than ordinals, so nothing else
needs to be mapped back to source coordinates.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.
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.
The extension’s SQL builder does not itself:
r/duckhtsbench producer
owns that transport for the registered Ensembl 116 GRCh38 artifact;estgene and otherfeatures/RefSeq;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.
.fai, and optional .gzi. Linux workers reopen those descriptors through
/proc/self/fd; Windows workers use the resolved source while the model retains a
deny-write handle. Other POSIX systems use independent handles on the resolved source and
verify their identities around lazy open, avoiding /dev/fd handles whose seek offsets
may be shared. External replacement or in-place mutation of a loaded source is outside the
immutable-model contract on those systems and fails when the identity change is observed.
Contained reference requests reuse the worker
window, and bounded forward read-ahead amortizes faidx fetches over coordinate-sorted
alleles. Workspaces are pooled only after their call completes.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.
The intended decision chain is:
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_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.
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.
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 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:
duckvep_model_open(...) is the single constructor for prepared transcript, exon,
sequence and optional interval-feature views used by consequence and HGVS. It derives
the first complete reference stop once per coding transcript, or validates a supplied
immutable cache, so per-row HGVSp
does not rescan an unchanged CDS;duckvep_transcript_edit_t is an HGVS-facing carrier built from that prepared allele and
one transcript. Projection first adds VEP’s endpoint-clipped transcript-slice
coordinates without trimming, reinterpreting, or clamping the semantic allele; CDS
projection is attached lazily only when reference validation or protein replay needs it;c. from non-coding n. edits,
preserve insertion and range geometry, and render into caller-owned bounded buffers;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:
g./m. HGVS is computed once per allele and needs a genomic reference-window
provider; VariationFeature::get_all_hgvs_genomic uses
get_3prime_seq_offset;c./n. HGVS is per transcript, but VEP first runs
TranscriptVariationAllele::_genomic_shift over a forward-reference window in that
transcript’s strand direction. Its perform_shift loop has a hard-coded 1000-base
search and retains its own allele-length-dependent loop limit. The shifted genomic
event is then projected back to transcript coordinates. Complete uploaded-REF checking
and hgvs_variant_notation duplication-source comparison are separate consumers of a
wider bounded reference lookup and must not widen that exact shift slice;perform_shift over the corrected transcript sequence.
VEP encounters this in some RefSeq models; it is a separate model capability, not a
property of RefSeq identifiers themselves; andp. HGVS consumes the alternate peptide difference produced by the same edit
set used for phased consequence classification.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.
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.
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.
make test-duckvep-kernel runs fixed and randomized pure-C properties; ASan and UBSan
execute the same suite.make test-duckvep-kernel-statistical raises randomized targets to an explicit seed and
trial count.DUCKVEP_PROPERTY_ARGS='-t hgvs' and '-t haplotype' are diagnostic filters for local
iteration. An official history run leaves this argument empty so no favorable family can
replace the complete state-machine denominator.test/duckvep/conformance/property_history.R records those counters in a long-form
append-only ledger, so a passing run cannot hide that a declared edit shape, strand,
terminal state, shift limit, or haplotype interaction received zero observations.
data/property_coverage_requirements.tsv classifies selected counters as statistically
required or covered by a named deterministic C witness. The history runner fails on a
missing counter, a statistically required zero, or any other zero without such a witness,
and preserves the complete seed-specific failure log.KI270395.1.make test-duckvep-differential compares generated witnesses to pinned VEP 116.make duckvep-corpus-differential records the union of emitted variant/transcript or
variant/core-feature pairs, including mismatches, misses, extras, and unresolved rows.--regulatory is
enabled, in addition to transcript-derived points. Independent seeds can run
concurrently against one read-only attached model.20260719 executed 100,000 trials per randomized
property (175 tests; 204,759 assertions) and passed under the ordinary,
AddressSanitizer, and UndefinedBehaviorSanitizer targets. It found and minimized a rare
terminal-codon oracle error after 93,064 generated frame-changing cases; the pinned VEP
source showed that the production kernel was correct, and the corrected oracle retained
both concrete-local-peptide and endpoint-reconstruction coverage. The independent
executable-VEP run with seed 20260716 compared 100,268 generated/fixed alleles and all
100,268 transcript pairs were exact. Its generated set contained SNVs, MNVs, insertions,
deletions, and delins; duplicate rejection deliberately makes the accepted shape counts
non-uniform rather than resampling a cosmetically balanced result.<CNV:TR> expansion, producer-specific symbolic encodings, or untested species.
The separate source-pinned CIPOS/CIEND contract retains inner/outer coordinates but
uses nominal POS/END for VEP-116 consequence terms. The checked-in 12-record
structural_confidence_grch38.vcf witness pairs nominal and imprecise forms of all six
executable exact single-locus kinds. It matched all 466 VEP transcript pairs, and each
engine produced equal nominal/imprecise consequence multisets for all six pairs while
the oracle VCF retained IMPRECISE, CIPOS, and CIEND.--regulatory comparisons cover a 1,196-site chromosome-21 GIAB sample
(14,955 annotation-object pairs) and 2,700 generated exact structural events
(120,224 pairs). Both are exact with no unresolved, extra, missing, or discordant row.
The structural corpus observes every one of the six regulatory/motif SO terms across
exact, containing, and partial feature geometries. These are declared sampled
distributions, not evidence for confidence-aware policy beyond VEP’s nominal-coordinate
terms or for untested funcgen releases.make bench-duckvep-release-parquet reads the official Ensembl variation consequence
VCF through typed CSQ columns and records complete versus consequence-only Parquet size,
checksum, cardinality, and elapsed time without committing the large artifacts.VE relation rather than the lossy CSQ presentation, and retain differences
as product-lineage evidence. Release 116 differs from cache-mode VEP at X/Y:276322 G>A,
so executable-VEP witnesses and generated corpora remain the semantic compatibility gate.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:
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.