Ensembl, Ensemble, and Your Own Genome: A Practical Annotation Guide

By the end of this you will have your own whole-genome VCF annotated with transcript consequences, population frequencies, and a set of predictor scores, queryable from DuckDB in a few hundred milliseconds, plus a consensus call set built from more than one variant caller so you know which calls are pipeline artifacts. You need: a WGS VCF and the matching BAM or CRAM, the reference FASTA the alignment was done against (this matters more than people expect), Docker or Singularity, about 250 GB of free disk, 16 GB of RAM, and 8 cores if you want annotation to finish in under an hour. Everything here runs offline on a laptop or a single cloud box. Nothing here is a clinical interpretation, and none of it substitutes for a clinical laboratory or a physician.

1. Two different things share the name

People search “ensemble DNA” and mean one of two things.

Ensembl (no “e” at the end) is the genome annotation resource run by EMBL-EBI. It maintains gene, transcript, and regulatory annotation for human and thousands of other species, a browser at ensembl.org, a REST API, a MySQL-accessible set of core databases, and the tool you will use most: the Variant Effect Predictor (VEP)¹. It is free. There is no account, no license fee, and no rate card. The code is open source, the data files are downloadable, and the REST API is open with a cap of roughly 15 requests per second. The taxonomic siblings are bacteria, plants, fungi, protists and metazoa. They live under Ensembl Genomes and use the same schema and the same tooling². If you only care about human, ignore them.

Ensemble methods in bioinformatics means combining several models or pipelines and using their agreement as the signal. Two places this shows up in personal genomics. First, ensemble variant calling: run GATK, DeepVariant, and Strelka2 on the same BAM and reconcile. Second, ensemble predictors: scores like REVEL, CADD, and the dbscSNV splice models are themselves classifiers trained on the outputs of other classifiers³. The same idea has been applied to promoter and transcription start site prediction, where combining models beat every individual model in the comparison⁴.

You want both. Ensembl gives you the coordinate system and the functional annotation. Ensemble methods give you a calibrated sense of which calls and which predictions to trust.

2. Check and normalize the VCF before annotating anything

Annotation on a messy VCF produces confident nonsense. Three things go wrong: the VCF is on a different assembly than you think, multi-allelic sites are unsplit, and indels are not left-aligned.

# What assembly and how many records?
bcftools head sample.vcf.gz | grep -E '^##(reference|contig=<ID=chr1,)'
bcftools stats sample.vcf.gz > sample.stats.txt
grep -E '^SN' sample.stats.txt

For a single human genome at 30x you should see roughly 4.5 to 5 million small variants, of which about 4 million are SNVs and 500,000 to 900,000 are indels. If you see 40,000, you have an exome. You do not have a genome. If you see 8 million, the call set is probably unfiltered or the sample is not the ancestry the joint-calling panel assumed.

Normalize:

bcftools norm -m -any -f GRCh38_full_analysis_set_plus_decoy_hla.fa \
    -c s sample.vcf.gz -Ou \
  | bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

-m -any splits multi-allelics into one record per ALT. -c s sets the REF to the reference and warns you loudly if they disagree, which is how you find out your VCF is GRCh37 in a GRCh38 costume. -d exact drops exact duplicates. Keep the intermediate files. If -c s prints thousands of REF mismatches, stop and fix the assembly before going further.

Contig naming (chr1 vs 1) has to match your cache and your FASTA. VEP will silently skip contigs it does not recognize, which shows up as a suspiciously small output file.

3. Install VEP offline with a cache

Run VEP from the container. Installing Perl dependencies by hand is not a good use of an afternoon.

docker pull ensemblorg/ensembl-vep:release_112.0

mkdir -p $HOME/vep_data && chmod a+rwx $HOME/vep_data

docker run -t -i -v $HOME/vep_data:/data ensemblorg/ensembl-vep:release_112.0 \
  perl INSTALL.pl -a cf -s homo_sapiens -y GRCh38 -c /data --CONVERT

-a cf installs the cache and the FASTA. --CONVERT builds the tabix-indexed cache, which is what makes --offline fast and lets you annotate a region without reading the whole cache. Budget about 20 GB for the human GRCh38 indexed cache and 30 to 60 minutes on a decent connection.

Use --offline for real work. Querying the Ensembl database over the network is fine for ten variants and unusable for five million. Pin the release number and record it next to your results: transcript sets change between releases, and a consequence call of splice_donor_variant in release 108 can become splice_region_variant in 112 because a transcript was retired.

4. Annotate, with the flags that matter

docker run -v $HOME/vep_data:/data ensemblorg/ensembl-vep:release_112.0 vep \
  --input_file /data/sample.norm.vcf.gz \
  --output_file /data/sample.vep.vcf.gz \
  --vcf --compress_output bgzip \
  --cache --offline --dir_cache /data --assembly GRCh38 \
  --fasta /data/homo_sapiens/112_GRCh38/Homo_sapiens.GRCh38.dna.toplevel.fa.gz \
  --everything \
  --mane_select --canonical \
  --pick_allele_gene --flag_pick_allele_gene \
  --buffer_size 50000 --fork 8 \
  --stats_file /data/sample.vep.html

What these do:

--everything is shorthand for about fifteen flags.

• SIFT
• PolyPhen
• HGVS notation
• gene symbols
• exon and intron numbers
• protein domains
• regulatory features
• canonical and MANE flags
• transcript support level
• gnomAD exome and genome allele frequencies
• variant class

It is the right default. It roughly triples runtime compared to bare consequence calling.

--pick_allele_gene with --flag_pick_allele_gene is the choice worth understanding. A variant in a gene with 14 annotated transcripts gets 14 consequence blocks. Unfiltered, a five-million-variant VCF expands to 60 million or more annotation rows. The --flag variant keeps every annotation but marks one per allele-gene pair as the picked one, so you can filter with a boolean instead of losing information. --pick alone throws the rest away, which you will regret the first time you need the non-canonical transcript.

--hgvs (included in --everything) requires the FASTA. Without it you get no HGVSc or HGVSp, and those are the strings you need to search literature or a variant database.

Runtime on 8 cores with a warm cache: 30 to 70 minutes for a 30x WGS VCF. Memory scales with --buffer_size times --fork. If you get OOM-killed, drop buffer size to 10000 first.

5. Add plugins for the scores you care about

Core VEP tells you what a variant does to a transcript. Plugins bring in the precomputed predictions, and this is where the ensemble predictors enter⁵.

Download dbNSFP, CADD whole-genome SNVs plus the indel file, and the SpliceAI precomputed VCFs. The dbNSFP download is one bgzipped, tabixed file, roughly 30 GB for the academic release. Then:

  --plugin dbNSFP,/data/dbNSFP4.7a_grch38.gz,REVEL_score,REVEL_rankscore,\
MetaSVM_pred,MutationTaster_pred,GERP++_RS,phyloP100way_vertebrate,\
ClinPred_score,AlphaMissense_score \
  --plugin CADD,/data/whole_genome_SNVs.tsv.gz,/data/gnomad.genomes.r4.0.indel.tsv.gz \
  --plugin SpliceAI,snv=/data/spliceai_scores.raw.snv.hg38.vcf.gz,\
indel=/data/spliceai_scores.raw.indel.hg38.vcf.gz \
  --plugin dbscSNV,/data/dbscSNV1.1_GRCh38.txt.gz \
  --plugin LoFtool \
  --plugin NMD

Notes from having broken each of these:

• dbNSFP only covers nonsynonymous and splice-site SNVs in coding regions. Asking it for an intronic variant returns nothing, and that is expected. It is not a bug.
• The column names in the plugin call must match the dbNSFP release exactly. Between 4.3a and 4.7a several columns were renamed. Check with zcat dbNSFP4.7a_grch38.gz | head -1 | tr '\t' '\n' | nl.
• SpliceAI precomputed files cover SNVs and small indels within genes. Anything outside that needs the model run live, which needs a GPU and a different pipeline.
• dbscSNV supplies the ada_score and rf_score, two ensemble classifiers trained specifically on splicing-altering variants near splice sites, with the reported advantage that they combine several underlying splice predictors rather than relying on one³. Treat scores near the published thresholds as undecided. Do not treat them as negative.

Do not add every plugin that exists. Each one adds a tabix lookup per variant and a column to a CSQ field that is already 40 fields wide.

6. Build a consensus call set across callers

Annotation quality is capped by call quality. Different pipelines on the same reads disagree more than you would guess, and the disagreement is concentrated in indels, homopolymers, and low-complexity regions. A comparison of whole-genome pipelines found that no single aligner-caller combination recovered everything, and that combining pipelines was what reduced false negatives⁶.

Run at least two callers. We use DeepVariant plus GATK HaplotypeCaller, and add Strelka2 when indels matter.

# DeepVariant
docker run -v "$PWD":/in google/deepvariant:1.6.1 /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/in/GRCh38.fa \
  --reads=/in/sample.cram \
  --output_vcf=/in/dv.vcf.gz \
  --num_shards=16

# Normalize both, then compare
for f in dv gatk; do
  bcftools norm -m -any -f GRCh38.fa -c s $f.vcf.gz -Ou \
    | bcftools norm -d exact -Oz -o $f.norm.vcf.gz
  bcftools index -t $f.norm.vcf.gz
done

bcftools isec -p isec_out -n +1 -c none dv.norm.vcf.gz gatk.norm.vcf.gz

isec_out/sites.txt gives you a per-site membership mask. Expect 95 to 98% concordance on SNVs and 85 to 92% on indels at 30x. If SNV concordance is below 93%, suspect the BAM: duplicate marking, a mismatched reference, or coverage lower than the report claims.

Then pick a policy and write it down:

• Intersection (-n=2) for anything you intend to look at closely. Higher precision, fewer variants.
• Union (-n +1) with a CALLERS tag when you are checking whether a specific known variant exists at all. Higher sensitivity, more junk.

Add the caller membership back as an INFO field so downstream queries can filter on it:

bgzip isec_out/sites.txt && tabix -s1 -b2 -e2 isec_out/sites.txt.gz
bcftools annotate -a isec_out/sites.txt.gz -c CHROM,POS,REF,ALT,INFO/CALLERS \
  -h <(echo '##INFO=<ID=CALLERS,Number=1,Type=String,Description="isec mask">') \
  merged.vcf.gz -Oz -o merged.tagged.vcf.gz

One thing consensus calling will not fix: regions no short-read pipeline covers well. Coverage is uneven, GC-extreme and repetitive regions drop out, and sensitivity in exons is worse from exome capture than from WGS at comparable depth⁷. Systematic biases in read length, error profile, and mappability are properties of the sequencing chemistry. They are not properties of your caller choices⁸. Structural variants, short tandem repeat expansions, and the HLA and pharmacogene regions with high paralogy need dedicated tools (Manta, ExpansionHunter, Hi-C or long reads). A better SNV consensus will not help.

7. Query the result

Flatten the CSQ field into columns and load it into DuckDB. Grepping a compressed VCF is not a query interface.

bcftools +split-vep merged.vep.vcf.gz -d -A tab \
  -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER\t%INFO/CALLERS\t%SYMBOL\t%Gene\t%BIOTYPE\t%Consequence\t%IMPACT\t%MANE_SELECT\t%HGVSc\t%HGVSp\t%gnomADg_AF\t%MAX_AF\t%CADD_PHRED\t%REVEL_score\t%SpliceAI_pred_DS_AG\t%SpliceAI_pred_DS_DL\t%ada_score\n' \
  -s worst \
  | bgzip > variants.tsv.gz

-d deduplicates repeated annotation values, -s worst selects the most severe consequence per variant, -A tab sets the field separator for multi-value fields.

CREATE TABLE v AS SELECT * FROM read_csv_auto('variants.tsv.gz', delim='\t', header=false);

-- rare, high-impact, called by both pipelines, on a MANE transcript
SELECT * FROM v
WHERE column10 = 'HIGH'
  AND (column15 IS NULL OR TRY_CAST(column15 AS DOUBLE) < 0.001)
  AND column11 <> ''
  AND column5 = '11';

A typical genome gives 100 to 200 rare predicted loss-of-function variants that survive that filter. Most are benign or in genes with no established phenotype. Carrier status for recessive conditions is the common finding, and it is a category where the difference between a research-grade annotation and a clinical result is large. Anything you plan to act on needs confirmation in an accredited clinical laboratory and interpretation by a clinical geneticist or genetic counselor. VEP consequence calls are hypotheses about transcript effects. They are not diagnoses.

8. Non-coding variants and regulatory annotation

About 98% of your variants are outside coding exons, and the coding-focused predictors have nothing to say about them. --everything turns on --regulatory, which annotates against Ensembl Regulatory Build features (promoters, enhancers, CTCF binding sites, open chromatin) derived largely from ENCODE and allied projects⁹. You get RegulatoryFeature and MotifFeature consequence types with a feature ID.

Use this to generate hypotheses. Do not use it to draw conclusions. The Regulatory Build is a union across cell types. A promoter-flanking annotation tells you the region is regulatory in some tissue. It does not tell you that it is regulatory in yours. If you have RNA sequencing from the same person, that is where non-coding variants become testable: check whether the gene near the candidate variant shows allele-specific expression or an expression level outside the reference range. Without matched expression data, a regulatory annotation on a rare non-coding variant is roughly as informative as its conservation score.

Common problems

REF mismatch warnings from bcftools norm -c s. Your VCF and FASTA are different assemblies or different patch releases. Do not liftover with a naive coordinate map if you can re-call instead. If you must lift, use bcftools +liftover or CrossMap with the UCSC chain file and then drop everything that fails REF validation.

VEP output has far fewer variants than the input. Contig naming mismatch. The cache uses Ensembl-style names (1, MT). Most alignment pipelines use UCSC-style (chr1, chrM). Add --chr handling or rename with bcftools annotate --rename-chrs.

Empty plugin columns everywhere. Either the plugin file is not tabix-indexed, or the plugin file assembly does not match --assembly, or you asked dbNSFP for a column name from a different release. Test on 100 variants (bcftools view -r chr1:1000000-2000000) before running the whole genome.

Annotation runs for six hours. You left --pick* off entirely and are writing every transcript for every variant, or --fork is 1. Check that the cache is the --CONVERT indexed version.

gnomAD frequency is missing for a variant you know is common. --af_gnomade and --af_gnomadg in the cache are frozen at the release’s gnomAD version and only cover the sites gnomAD reports. For anything current, annotate from the gnomAD VCFs directly with bcftools annotate -a gnomad.genomes.v4.1.sites.chr1.vcf.bgz -c INFO/AF,INFO/AF_nfe.

Two callers agree on a variant and it is still wrong. Correlated errors. Both callers use the same alignments, so a mismapping in a segmental duplication fools both. Load the BAM in IGV and look at the reads: mapping quality zero, soft-clipped ends, or a pileup of identical mismatches are the tells. Visual inspection of the reads is still the fastest check on any single variant you care about.

Cost questions. Ensembl and VEP cost nothing. Sequencing does, and the “$1000 genome” was a target set in the mid-2000s in terms of reagent and instrument cost, not a delivered price with analysis, storage, and interpretation attached¹⁰. Consumer sequencing prices today range widely depending on depth, whether you get the CRAM or only a VCF, and whether anything is clinically reported. Insist on the raw alignment file. A VCF alone forecloses re-calling, consensus work, and every step in section 6.

Footnotes

1. William McLaren, Laurent Gil, Sarah E. Hunt, et al. The Ensembl Variant Effect Predictor. Genome Biology, 2016. https://doi.org/10.1186/s13059-016-0974-4 ↩

2. P. J. Kersey, D. Lawson, E. Birney, et al. Ensembl Genomes: Extending Ensembl across the taxonomic space. Nucleic Acids Research, 2009. https://doi.org/10.1093/nar/gkp871 ↩

3. Xueqiu Jian, Eric Boerwinkle, Xiaoming Liu. In silico prediction of splice-altering single nucleotide variants in the human genome. Nucleic Acids Research, 2014. https://doi.org/10.1093/nar/gku1206 ↩ ↩²

4. Hong-Hee Won, Min-Ji Kim, Seonwoo Kim, et al. EnsemPro: An ensemble approach to predicting transcription start sites in human genomic DNA sequences. Genomics, 2008. https://doi.org/10.1016/j.ygeno.2007.11.001 ↩

5. Michael Yourshaw, S. Paige Taylor, Aliz R. Rao, et al. Rich annotation of DNA sequencing variants by leveraging the Ensembl Variant Effect Predictor with plugins. Briefings in Bioinformatics, 2014. https://doi.org/10.1093/bib/bbu008 ↩

6. Kyu-Baek Hwang, In-Hee Lee, Honglan Li, et al. Comparative analysis of whole-genome sequencing pipelines to minimize false negative findings. Scientific Reports, 2019. https://doi.org/10.1038/s41598-019-39108-2 ↩

7. Alison M Meynert, Morad Ansari, David R FitzPatrick, et al. Variant detection sensitivity and biases in whole genome and exome sequencing. BMC Bioinformatics, 2014. https://doi.org/10.1186/1471-2105-15-247 ↩

8. Martin Kircher, Janet Kelso. High‐throughput DNA sequencing – concepts and limitations. BioEssays, 2010. https://doi.org/10.1002/bies.200900181 ↩

9. The ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature, 2012. https://doi.org/10.1038/nature11247 ↩

10. Simon T Bennett, Colin Barnes, Anthony Cox, et al. Toward the $1000 Human Genome. Pharmacogenomics, 2005. https://doi.org/10.1517/14622416.6.4.373 ↩