Variant effect prediction as a SQL relation. Load an Ensembl transcript model into DuckDB once, then query duckvep_annotate_sql(...) as a relation to filter, join and aggregate alongside VCF, ClinVar or population-frequency data. A new question is a new query, not another pipeline run.
Variant annotation take time, too much time, variant_myth try to solve this problem using modern and parallel tools.
— rationale of variant_myth
DuckVEP began with that sentence. variant_myth answers it with a parallel Rust annotator built on interval sets and hash tables; DuckVEP pursues the same question inside the database, where annotation is a relation you can query.
DuckVEP returns typed Sequence Ontology consequences and impact, cDNA/CDS/protein positions and amino acids, NMD predictions, and HGVS. Phased calls also support coding whole-haplotype consequences.
Independent-event consequences and HGVS are compared with executable Ensembl VEP 116 and its matching core and variation databases. Discordant and unresolved pairs stay in the counts, deliberate differences are documented, and unavailable results carry a status and reason.
Install · Try it · Function reference · Conformance · Errata · Design · R package
DuckVEP is written in C against the stable DuckDB C API v1, which is the default build. That build is qualified on released DuckDB 1.5.x: all 35 SQL test files pass on DuckDB 1.5.6 (source 069cc9f9b5), and earlier runs cover the DuckDB 1.5.1 CLI and R duckdb 1.5.5. A DuckDB community-extensions entry for this build is in preparation; it is not published, so INSTALL duckvep FROM community does not work yet. Build from source and load the unsigned extension:
git clone --recurse-submodules https://github.com/RGenomicsETL/DuckVEP
cd DuckVEP
make configure release test_release
duckdb -unsigned -c "LOAD 'build/release/duckvep.duckdb_extension'"The extension links its own htslib (for indexed reference FASTA) and cgranges. LOAD registers functions and fetches nothing; transcript models are built from Ensembl data you supply, or restored from snapshot files you saved. From R, Rduckvep builds the same sources offline and adds connection, model and haplotype helpers.
A separate C API v2 host is an optional preview, validated against a pinned DuckDB 2.0 development snapshot. It is neither the default build nor the community artifact; changing the default waits for released DuckDB 2.0.0 (#48).
DuckHTS releases through 1.5.2 register legacy DuckVEP functions under the same SQL names. Loading one alongside DuckVEP can split model loading and annotation across separate model registries, producing unknown model name. Use DuckHTS newer than 1.5.2 with DuckVEP.
A model is loaded from ordinary relations: regions, transcripts with their CDS and flanking sequence, and exons. In production these come from an Ensembl core database through query(duckvep_ensembl_regions_sql(...)) and query(duckvep_ensembl_transcripts_sql(...)). Here a one-transcript fixture stands in:
SELECT loaded FROM duckvep_model_load('demo',
'SELECT * FROM readme_regions ORDER BY seq_region',
'SELECT * FROM readme_transcripts ORDER BY seq_region, transcript_start',
'SELECT * FROM readme_exons ORDER BY transcript_index, exon_start');
#> ┌─────────┐
#> │ loaded │
#> │ boolean │
#> ├─────────┤
#> │ true │
#> └─────────┘The named model stays resident in the DuckDB session. Annotate seven variants across the transcript, from the 5′ UTR to the intron, and join the results back to their input in one query:
CREATE TABLE demo_events AS
SELECT row_number() OVER (ORDER BY position, alternate)::UBIGINT AS event_index,
1::UINTEGER AS seq_region, position::UBIGINT AS position, reference, alternate,
NULL::UBIGINT AS end_position, NULL::VARCHAR AS structural_type,
NULL::VARCHAR AS copy_change, NULL::UINTEGER AS mate_seq_region,
NULL::UBIGINT AS mate_position
FROM (VALUES (110, 'C', 'T'), (124, 'T', 'C'), (125, 'A', 'G'), (130, 'C', 'CA'),
(134, 'C', 'A'), (151, 'G', 'A'), (175, 'A', 'G')) v(position, reference, alternate);
SELECT e.position, e.reference || '>' || e.alternate AS change,
a.consequence, a.impact, a.transcript_hgvs, a.protein_hgvs, a.nmd_prediction
FROM query(duckvep_annotate_sql('demo_events', 'demo', struct_pack(rich := true, hgvs := true))) a
JOIN demo_events e USING (event_index)
ORDER BY e.position;
#>
#> ┌──────────┬─────────┬──────────────────────┬──────────┬─────────────────┬──────────────┬────────────────┐
#> │ position │ change │ consequence │ impact │ transcript_hgvs │ protein_hgvs │ nmd_prediction │
#> │ uint64 │ varchar │ varchar │ varchar │ varchar │ varchar │ varchar │
#> ├──────────┼─────────┼──────────────────────┼──────────┼─────────────────┼──────────────┼────────────────┤
#> │ 110 │ C>T │ 5_prime_UTR_variant │ MODIFIER │ NULL │ NULL │ NULL │
#> │ 124 │ T>C │ missense_variant │ MODERATE │ c.5T>C │ p.Val2Ala │ NULL │
#> │ 125 │ A>G │ synonymous_variant │ LOW │ c.6A>G │ p.Val2= │ NULL │
#> │ 130 │ C>CA │ frameshift_variant │ HIGH │ NULL │ NULL │ escaping │
#> │ 134 │ C>A │ stop_gained │ HIGH │ c.15C>A │ p.Tyr5Ter │ escaping │
#> │ 151 │ G>A │ splice_donor_variant │ HIGH │ NULL │ NULL │ unresolved │
#> │ 175 │ A>G │ intron_variant │ MODIFIER │ NULL │ NULL │ NULL │
#> └──────────┴─────────┴──────────────────────┴──────────┴─────────────────┴──────────────┴────────────────┘The frameshift and the stop gain escape nonsense-mediated decay because they fall in the first CDS positions (nmd_escape_early_cds).
Some HGVS fields are empty, and that isn’t a parsing failure. This model was loaded without a reference FASTA, and those notations need flanking sequence to shift and render. DuckVEP says so rather than guess, and asking why is just another query:
SELECT e.position, a.transcript_hgvs_status, a.transcript_hgvs_reason
FROM query(duckvep_annotate_sql('demo_events', 'demo', struct_pack(hgvs := true))) a
JOIN demo_events e USING (event_index)
WHERE a.transcript_hgvs IS NULL
ORDER BY e.position;
#> ┌──────────┬────────────────────────┬────────────────────────┐
#> │ position │ transcript_hgvs_status │ transcript_hgvs_reason │
#> │ uint64 │ varchar │ varchar │
#> ├──────────┼────────────────────────┼────────────────────────┤
#> │ 110 │ unresolved │ missing_reference │
#> │ 130 │ unresolved │ missing_reference │
#> │ 151 │ unresolved │ missing_reference │
#> │ 175 │ unresolved │ missing_reference │
#> └──────────┴────────────────────────┴────────────────────────┘Phased genotypes become whole-haplotype consequences, again as a relation. The changes at 124 and 125 fall in the same codon (Val2). On their own, one is missense and the other synonymous. Carried on the same haplotype, they are read together as one codon, GTA to GCG, and classified once:
CREATE TABLE demo_calls AS
SELECT event_index, seq_region, position, reference, alternate,
1::UINTEGER alt_index, 0::UINTEGER transcript_index, 0::UBIGINT sample_index,
[1, 0]::INTEGER[] alleles, [false, true]::BOOLEAN[] phase_before, NULL::BIGINT phase_set
FROM demo_events WHERE position IN (124, 125);
SELECT carrier_count, length(contributors) AS contributors, prediction_status,
haplotype_consequences, haplotype_impact
FROM duckvep_haplotypes('SELECT * FROM demo_calls', 'demo');
#>
#> ┌───────────────┬──────────────┬───────────────────┬────────────────────────┬──────────────────┐
#> │ carrier_count │ contributors │ prediction_status │ haplotype_consequences │ haplotype_impact │
#> │ uint32 │ int64 │ varchar │ varchar[] │ varchar │
#> ├───────────────┼──────────────┼───────────────────┼────────────────────────┼──────────────────┤
#> │ 1 │ 2 │ predicted │ [missense_variant] │ MODERATE │
#> └───────────────┴──────────────┴───────────────────┴────────────────────────┴──────────────────┘For VCF or BCF input, duckvep_coding_calls(model, path) builds the calls relation directly; its model’s seq_region_name values must match VCF CHROM. It discovers coding transcripts before decoding GT and PS, then feeds duckvep_haplotypes. If you stage calls yourself, duckvep_coding_transcripts(model, seq_region, position, reference, alternate) returns the transcripts touched by the normalized coding event, so a shared indel anchor alone does not count. R wrappers are Rduckvep::rduckvep_coding_calls() and Rduckvep::rduckvep_coding_transcripts().
Across 9 recorded corpora, executable VEP 116 comparisons contain 1,486,561 exact of 1,486,561 transcript and feature pairs, with 0 unresolved and 0 discordant pairs. The corpora cover GRCh37 and GRCh38, Plasmodium falciparum genetic codes, small variants, exact structural events, paired breakends, regulatory features and motifs. The corpus release criterion is zero discordant, missing, extra or unresolved pairs.
The same evidence set records 357,806 exact HGVSc and 357,806 exact HGVSp transcript pairs against VEP --hgvs; 56,998 pairs come from ClinVar chromosome 21, with 0 discordant. Methods and term-level results are in the conformance report; the consequence and HGVS ledgers hold per-run counts and source revisions. Ledger revisions preceding extraction from DuckHTS resolve through the commit map.
Per-corpus consequence counts
| Corpus | Assembly | Oracle | Pairs | Exact | Unresolved | Different |
|---|---|---|---|---|---|---|
| dbSNP 157 windows | GRCh38 | VEP 116 | 73,620 | 73,620 | 0 | 0 |
| GIAB HG002 small variants | GRCh38 | VEP 116 | 54,905 | 54,905 | 0 | 0 |
| ClinVar coding | GRCh38 | VEP 116 | 287,836 | 287,836 | 0 | 0 |
| ClinVar, all chromosomes | GRCh38 | VEP 116 | 316,397 | 316,397 | 0 | 0 |
| GRCh37 cache corpus | GRCh37 | VEP 116 | 486,464 | 486,464 | 0 | 0 |
| P. falciparum (Ensembl Genomes 63) | GCA_000002765v3 | VEP 116 | 40,732 | 40,732 | 0 | 0 |
| Paired breakends, multi-chromosome | GRCh38 | VEP 116 | 91,428 | 91,428 | 0 | 0 |
| GIAB + regulatory and motif features | GRCh38 | VEP 116 | 14,955 | 14,955 | 0 | 0 |
| Exact structural variants + regulation | GRCh38 | VEP 116 | 120,224 | 120,224 | 0 | 0 |
VEP 116 does not define whole-haplotype consequences. The duckvep-coding contract uses pinned bcftools csq where transcript models and phase semantics are comparable. Independent base-R goldens and hand-derived checks cover supported cases outside csq’s domain: non-diploid calls, alleles over 50 bases, nonstandard genetic codes, and incomplete CDS starts or ends. Unsupported inputs retain explicit statuses and reasons.
Bounded extensions cover compound literal DNA HGVS for non-overlapping edits with observed-cis evidence, explicit diploid phase alternatives, and equal-length exonic noncoding/UTR replay against the final transcript allele. The arrangement relation enumerates admitted unphased calls as labelled hypotheses, separate from strict observed-cis replay. Typed curated translations expose separate, explicitly conditional reference/alternate peptides; they do not assert preserved biological recoding competence or replace raw compatibility proteins. Overlapping edits, general splice and alternative-start biology, and untyped or RNA edits remain unresolved and carry statuses. See the function reference for inputs and limits.
The pinned full-HG002 qualification for issue #34 recorded one-core, cold-process medians against the Ensembl 116 model: 17.25 s for csq, 8.21 s for the fused-reader CLI path (2.10×, including DuckVEP model load), and 8.60 s for the R path (2.01×). The qualification method and receipts identify the input and run contract.
Compatibility is claimed only where it has been measured, and differences are kept, not hidden. The compatibility and errata record separates three verdicts, each with its witness:
- a VEP convention DuckVEP must follow;
- a gap in DuckVEP;
- a likely upstream erratum.
The main cases:
-
Published release annotations are not the oracle. Ensembl’s release VCFs sometimes disagree with executable VEP. In release 116,
X/Y:276322 G>Ais published asintergenic_variant, while the executable emits three5_prime_UTR_variantrows per chromosome. DuckVEP follows the executable, and the PAR witnesses pin both. -
Breakends are evaluated one event at a time. VEP 116’s buffered breakend path uses a chromosome-blind interval tree and can drop valid transcript pairs in multi-chromosome batches. The oracle runs with
--buffer_size 1to isolate that; DuckVEP has no such batching effect (paired-breakend differential). -
Imprecise structural variants use the nominal span, as VEP’s registered predicates do.
CIPOSandCIENDstay on the row as metadata rather than being dropped. -
Structural-span comparison: VEP skips events above
--max_sv_size(5 kb by default). For supported exact structural events, DuckVEP uses the nominal span without that cutoff; conformance comparisons set VEP’s limit to 10 Mb. -
Circular regions use lifted intervals. Origin-crossing transcripts, exons and regulatory or motif features preserve overlap and HGVS behavior across the origin. This is separate from mitochondrial translation; human MT has no wrapped object. VEP cannot oracle origin-crossing models, so they are checked by rotation equivariance and linear-model differential (design, survey).
-
Unknown is a value. When a result can’t be computed, the row carries
duckvep_statusand a reason instead of a guess. Examples are a missing reference sequence, a reference mismatch, or an unsupported symbolic allele. -
Outside the claim:
- species and releases that haven’t been tested;
- compound HGVS beyond the bounded non-overlapping, observed-cis profile above, where the errata record the disagreements that remain;
- genomic and structural HGVS beyond the exact-span, unshifted structural preparation profile.
See the design contract and the roadmap.
Recorded benchmark results use their stated hosts, inputs, output contracts and thread counts; they describe those runs.
- SQL annotation: a one-core compact-output run on GIAB HG002 measured 1,034,768 alleles per second across 4,095,611 model-addressable alleles, with 1,383,580 regulatory and motif features resident. The measurement ran through
FROM query(duckvep_annotate_sql(...)); see the throughput method and receipts. - Comparison with vep-rs: three paired GIAB HG002 runs pinned to one E-core measured median process times of 58.09 s for DuckVEP and 73.26 s for vep-rs, including disk-backed outputs on a shared host. Their common consequence tuples agree on 34,146,531 of 34,148,222 comparisons; executable VEP agrees with DuckVEP on 1,689 of the 1,691 disagreements. These are distinct output formats and measured conditions, not a general speed ranking. See the plots, receipts and reproducible runner.
- Model restore:
duckvep_model_savewrites a loaded model snapshot andduckvep_model_restoremaps it read-only. In an Ensembl 116 GRCh38 measurement, restore took about half a second, and processes restoring the same file share its mapped pages. - Memory limits: native allocations are charged against an enforced budget; exceeding it raises a capacity error rather than truncating results. The scale runner records three concurrent 5-million-allele gnomAD jobs on one 20-thread host.
Rduckvep runs the same native annotation builders on your DuckDB connection. Builder functions return SQL, and rduckvep_annotate() executes it and returns a data frame. Annotation queries can read relations on that connection, including TEMP tables and uncommitted rows.
con <- Rduckvep::rduckvep_connect()
# After loading a model named "demo" and creating "demo_events":
sql <- Rduckvep::rduckvep_annotate_sql(con, "demo_events", "demo", hgvs = TRUE)
results <- Rduckvep::rduckvep_annotate(con, "demo_events", "demo", hgvs = TRUE)
DBI::dbDisconnect(con, shutdown = TRUE)DuckVEP needs no other extension, but it is meant to sit in a query next to them. With DuckHTS reading the VCF and DuckClinVarbitration supplying arbitrated ClinVar decisions, annotation and clinical evidence meet in one query. This one isn’t evaluated here, because it needs a full Ensembl model and both extensions:
-- One row per ALT allele, with the model's region index (DuckHTS reads the VCF)
CREATE TABLE alleles AS
SELECT row_number() OVER (ORDER BY r.seq_region, v.POS, alt_allele)::UBIGINT AS event_index,
r.seq_region, v.POS::UBIGINT AS position, v.REF AS reference, alt_allele AS alternate,
NULL::UBIGINT AS end_position, NULL::VARCHAR AS structural_type,
NULL::VARCHAR AS copy_change, NULL::UINTEGER AS mate_seq_region,
NULL::UBIGINT AS mate_position
FROM read_bcf('cohort.vcf.gz') v CROSS JOIN unnest(v.ALT) AS t(alt_allele)
JOIN grch38_regions r -- Ensembl names contigs 1..22, X, Y, MT; UCSC-style VCFs say chr1, chrM
ON r.seq_region_name = CASE WHEN v.CHROM IN ('chrM', 'M') THEN 'MT'
ELSE regexp_replace(v.CHROM, '^chr', '') END;
-- High-impact consequences next to their arbitrated ClinVar decision
SELECT r.seq_region_name AS contig, e.position, e.reference, e.alternate,
a.consequence, a.protein_hgvs, d.policy_classification, d.gold_stars
FROM query(duckvep_annotate_sql('alleles', 'human_116_grch38', struct_pack(rich := true, hgvs := true))) a
JOIN alleles e USING (event_index)
JOIN grch38_regions r USING (seq_region)
LEFT JOIN clinvar_vcf v -- DuckClinVarbitration
ON v.assembly = 'GRCh38' AND v.contig = r.seq_region_name AND v.position = e.position
AND v.reference = e.reference AND v.alternate = e.alternate
LEFT JOIN clinvar_policy_allele_decisions d USING (allele_id)
WHERE a.impact = 'HIGH';- Reference: SQL functions, with signatures, options, columns and checked examples; R package; host and release profiles
- Contracts: design and implementation (model build, the sweep, NMD, phased edits, structural events, HGVS, supplementary annotation); haplotype contract
duckvep-coding; structural identity and HGVS; corpus workflow - Evidence: conformance and its harness; compatibility and errata, every known difference from VEP classified with witnesses; haplotypes and their scale qualification
- Performance: throughput, vep-rs comparison, scale runner
The project site publishes these pages.
The question DuckVEP set out to answer is the one variant_myth asks, quoted at the top of this page; thank you to its author for asking it. Correctness here rests on two other projects, executable Ensembl VEP and, for whole haplotypes, pinned bcftools csq; both make a runnable oracle possible.
DuckVEP was developed within DuckHTS and extracted with its Git history.
GPL-2.0-or-later; see LICENSE. Vendored htslib and cgranges keep their own licences. Ensembl data and VEP are © EMBL-EBI under their own terms.