Analysis of variant data resulting from genome or exome sequencing is fundamental for progress in biology, from basic research to translational genomics in the clinic. It is key for investigating function and for progressing from a system of medical care based on standardized treatment to one targeted to the individual patient.

For sufferers of common or rare disease, the potential benefits of variant analysis include improving patient care, surveillance, and treatment outcomes. In cancer, there have already been numerous successes using data from genetic tests. For example, patients testing positive for the inheritance of BRCA mutations have the option of selective preventative surgery; lung cancer patients showing EGFR gene mutations or triple negative breast cancer patients can have their drug prescriptions tailored to improve success [1, 2].

Rare diseases can individually be difficult to diagnose due to the low incidence and the incomplete penetrance of implicated alleles. However, variant analysis of whole-genome sequencing (WGS) or whole-exome sequencing data can lead to the discovery of underlying genetic mutations [3]. Identifying an associated mutation is advantageous for researching treatment options and for future drug discovery. Meanwhile, even the immediate benefit of diagnosis may result in a more accurate prognosis and remove the burden of additional medical investigations.

The most common non-infectious diseases worldwide are cardiovascular disease, cancer, and diabetes [4]. Despite many array-based genome-wide association studies (GWAS) searching for risk loci, only a relatively small heritable component in these conditions has been elucidated [5]. WGS in large numbers of samples is required to yield enough statistical power to detect rare variants with potential phenotypic or disease associations [6, 7]. WGS studies will also detect variants in regulatory and non-coding regions of the genome, which are thought to comprise the majority of trait-associated variants [8] and play a role in cancer [9].

The potential of large-scale sequencing and variant analysis is revolutionary. Recognizing this value, major population sequencing initiatives have been launched in Iceland [10], the UK [11], and the USA [12]. In other species, efforts such as Genome 10 K [13], the 1001 Arabidopsis genomes [14], and 1000 bull genome project [15] have similar goals but operate under different funding models, often with less support than the Homo sapiens-focused projects.

Ongoing improvements in DNA sequencing technology, and a current cost around $1000 per human genome, have resulted in high volumes of genome, exome, and subsequent variant data requiring interpretation. Meanwhile, the cost of the analysis to determine functional consequences remains substantially higher due to the difficulty of variant interpretation. For example, a typical diploid human genome has around 3.5 million single nucleotide variants (SNVs) and 1000 copy number variants [16] with respect to the genome reference sequence. Around 20,000–25,000 of these variants are protein coding, of which 10,000 change an amino acid but only 50–100 are protein truncating or loss of function variants [16]. Manual review of large numbers of variants is impractical and costly and there are additional difficulties, such as a lack of functional annotation or the interpretation of multiple variants within a haplotype.

Variant interpretation often considers the impact of a variant on a transcript or protein. It is dependent, therefore, on transcript annotation and localizing variants to protein-coding or non-coding regions. There are two major sources of Homo sapiens annotation: GENCODE [17] and Reference Sequence (RefSeq) [18] at the National Center for Biotechnology Information (NCBI). Both sets of transcript annotation are subject to version changes and updates that can modify variant reporting and interpretation. For data reproducibility, transcript isoforms and transcript versions must be rigorously tracked, although in some cases even including the version is not sufficient to avoid all potential misinterpretations [19]. There are differences in how the transcript sets are produced: GENCODE annotation is genome-based while RefSeq transcripts are independent of the reference genome. Although RefSeq transcripts may correct for errors in the reference assembly and provide transcripts with improved biological representation (such as for the genes ABO, ACTN3, and ALMS1 in the GRCh37 reference), differences between a genome and a transcript set can cause confusion and errors when reporting variants at the cDNA and genomic levels (e.g., these descriptions refer to the same variant: NM_000059.3:c.7397C>T, NC_000013.11:g.32355250T=). GENCODE’s aim is to create a comprehensive transcript set to represent expression of each isoform across any tissue and stage of development and, as a result, there are, on average, nearly four transcript isoforms per protein-coding gene. Most genes, therefore, have several annotations for a given variant due to multiple transcript isoforms (the G protein-coupled receptor 56 gene (GPR56) in Ensembl release 79 has 61 transcripts). This number will increase as more experimental data accumulate. Choosing the correct transcript isoform and version for consistent variant annotation is challenging. Finally, in loci where the reference genome has several alternative haplotype representations (“ALTs”), variants may have different interpretations with respect to different ALTs. For example, rs150580082 has mappings to multiple ALTs but introduces a stop codon in only some of these. In this case, considering the primary assembly mapping alone will give misleading results.

Variant reporting using Human Genome Variation Society (HGVS) nomenclature is also based on transcripts or proteins. Therefore, the difficulties with transcript annotation described above may cause confusion and ambiguities when using HGVS nomenclature. Many possible annotations exist for variants in genes with multiple transcript isoforms. For example, rs121908462 is a pathogenic variant associated with polymicrogyria that falls in ADGRG1, an adhesion G protein-coupled receptor G1. This variant has 126 HGVS descriptions in Ensembl [20] (and even more valid HGVS descriptions exist), as it overlaps 75 transcripts, and another 103 different descriptions in dbSNP. Multiple transcripts per locus result in greater numbers of annotations. These require filtering in a consistent manner, which increases the instability and complexity of variant interpretation.

The VEP has been used for analysis of traits in farm animals [24, 25], for patient diagnosis in the clinic and for research on GWAS [26–30]. It has been used for analysis in numerous large-scale projects, including the 1000 Genomes [31] and Exome Aggregation Consortium (ExAC) [32]. VEP’s annotations are used as input to tools for deep exploration of variant annotation such as GEMINI [33]. It is a flexible tool of value to any project requiring detailed annotation of sequence variants.

The VEP annotates two broad categories of genomic variant: (1) sequence variants with specific and well-defined changes (including SNVs, insertions, deletions, multiple base pair substitutions, microsatellites, and tandem repeats); and (2) larger structural variants (greater than 50 nucleotides in length), including those with changes in copy number or insertions and deletions of DNA. For all input variants, the VEP returns detailed annotation for effects on transcripts, proteins, and regulatory regions. For known or overlapping variants, allele frequencies and disease or phenotype information is included.

The VEP can be used to analyze data from any species with an assembled genome sequence and an annotated gene set. The data files necessary for annotation in 80 vertebrate species and many invertebrates are distributed by Ensembl and Ensembl Genomes [34], respectively. These are updated regularly, ensuring analysis can be performed using contemporary biological knowledge. The VEP also supports both the latest GRCh38 and previous GRCh37 human assemblies. Importantly, all results are fully reproducible using Ensembl archived versions. Finally, researchers may use their own transcript data for analysis, e.g., in species not yet in Ensembl or for novel or private annotations. A script is included in the VEP script package to create an annotation set from a general feature format (GFF) and FASTA file pair.

Each version of the VEP is tied to a specific release of Ensembl. This explicit versioning ensures all results are stable across a release, which is critical for provenance and reproducibility. To avoid misinterpretation of a variant based on a previous transcript or protein version, the output includes the identifier and version in HGVS coding descriptions. The VEP is open source, free to use, and actively maintained and developed. A mailing list [35] provides responsive support and the benefits of a shared community. The wide usage helps ensure bugs are found and corrected rapidly and enables suggestions to be gathered from a broad range of project teams.

The nature of the VEP results are described below along with input and output formats, the different interfaces, and details on performance.

VEP interfaces

The VEP is platform independent and available as (1) an online tool, (2) an easily installed Perl script, or (3) via the Ensembl Representational State Transfer (REST) application program interface (API) [72]. Each interface is optimized to support different quantities of data and levels of bioinformatics experience. All three use the same core codebase to ensure results are consistent across each interface. A comprehensive test suite backs all code, with continuous integration performed by Travis CI [73], ensuring high quality code, which must pass stringent quality tests before release.

VEP Web

VEP Web [21] offers a simple point-and-click interface. This is ideal for exploring annotation in an interactive manner. The portal is most suited to first-time use or small-scale analysis. The maximum compressed uploaded data file size currently supported is 50 megabytes, large enough for around two million typical lines of VCF data.

For single variant analysis, the web interface incorporates ‘Instant VEP’ functionality. Pasting or typing a single variant such as a variant in HGVS notation from a manuscript will rapidly return basic consequence prediction data. To submit a request for more than one variant, data can be uploaded, pasted or given via URL and options selected using a simple online form. A limited set of the VEP’s most commonly used plugins is available to use via the web interface. Requests are processed by a resource management system on the Ensembl web servers to distribute the request load.

The output web page (see example in Fig. 1) shows summary statistics and charts to provide an overview of the results. It also has a table with a preview of the detailed results, with a simple interface to configure filtering of the output. Via a series of drop-down menus, multiple filters (see examples in Table 7) can be combined using basic logical relationships, thereby allowing the creation of complex customized queries. This is designed to aid prioritization of smaller numbers of variants. Results can be stored by logging into an Ensembl account.

VEP script

The downloadable Perl script [74] is the most powerful and flexible way to use the VEP. It supports more options than the other interfaces, has no limit on input file size, and includes extensive input, output, filtering, and analysis options.

To install the script, simply download the VEP package and run the installer script, which automatically downloads the necessary API and annotation files (or ‘cache’ files). Updates with the latest data are available for each Ensembl release. The full source code is freely available on the Ensembl GitHub repository.

To process large volumes of data, the VEP script works most efficiently in “offline” mode using a local cache of transcript annotations rather than online public databases. As well as optimizing runtime, this ensures data privacy for clinically or commercially sensitive data. Furthermore, the VEP input can be configured to query overlaps with local, potentially private, variant and phenotype data or other custom data sets in a manner similar to vcfanno [75]. In this way annotation in formats including BED, GFF, GTF, VCF, and bigWig can be incorporated into the VEP output.

Advanced filtering options are available for a smaller result set, either during runtime or as a post-run process (Table 7). Filtering can be performed as a post-run process by an accompanying script that uses a simple field-operator-value language. Filtered results can be fed back to the VEP for further analysis or exported.

With some familiarity of Perl, the VEP can truly be customized, extended, and integrated with other systems. As almost all of the algorithmic content of the VEP is contained within the Ensembl API, the features of the VEP can be accessed using API calls. It is trivial, therefore, to extend the VEP results and perform secondary analyses, such as retrieving all OMIM IDs associated with the genes in the VEP results or calculating known variants in linkage disequilibrium with a subset of variants. Alternatively, the VEP is also customizable via its plugin architecture, which was developed to provide greater scope for expansion. This architecture supports the use of VEP as the backbone of a customized analysis pipeline by writing additional code to extend the VEP’s functionality for specific use cases. Example uses include filtering output, adding annotation from local or remote sources, executing external programs, or rendering graphical representations of the output. Ensembl provides a number of VEP plugins, hosted on GitHub [76], and a variety are published [51, 77] (Table 4).

VEP REST API

Ensembl’s language-independent REST API provides robust computational access in any programming language and returns basic variant annotation and consequence data. Individually or in batches of up to 1000, variants can be submitted to the API server in a single request. Results return in JSON, simple for parsing in most modern programming languages (see Fig. 2 for an example of JSON output). Using this interface, dynamic VEP queries can be integrated into custom-built software for on-demand results, as used, e.g., in the Decipher Genome Browser [78]. For documentation see [79].

As with the web interface, a limited set of the VEP’s most commonly used plugins is configured for use via the REST API.

The Ensembl Variant Effect Predictor software provides tools and methods for a systematic approach to annotate and prioritize variants in both large-scale sequencing projects and smaller analysis studies. By automating annotation in a standard manner and reducing the time required for manual review, it helps manage many of the common challenges associated with analysis of SNVs, short insertions–deletions, copy number variants, and structural variants. The VEP annotates variants using a wide range of reference data, including transcripts, regulatory regions, frequencies from previously observed variants, citations, clinical significance information, and predictions of biophysical consequences of variants.

The quality, quantity, and stability of variant annotation obtained depends on the choice of transcript set used [82]. As such, the VEP allows flexibility of transcript choice. To effectively manage large numbers of variant annotations and transcript isoforms, the VEP provides several methods to prioritize results and reduce the number of variants needing manual review. A selection of these filters is available and VEP also supports building of custom filters. Uniquely, the VEP algorithm can be expanded to perform additional calculations via plugins [77] and can analyze custom, potentially private, data.

Interpreting all variants in a genome remains an unsolved challenge. An increasing number of large-scale WGS will detect rare variants in both coding and non-coding regions of the genome and further possible identification of loci associated with disease. Having these variants available in public repositories such as dbSNP and the European Variant Archive or discoverable using federated resources will be of significant benefit for analysis. Emerging efforts such as the Global Alliance for Genomic Health (GA4GH) Beacon project [83] are currently developing possible distributed solutions.

Improved functional annotation is especially critical for variants in non-coding regions. Many fall in loci that regulate gene expression in specific tissues. Characterizing associations between transcripts and tissues will facilitate a subset of tissue-specific transcript isoforms to be selected for variant annotation, tailoring results. Moreover, upon providing the link from regulatory region to regulated gene, the potential molecular mechanism underlying disease could be explained. Data from large scale efforts such as the Genotype-Tissue Expression project, which aims to systematically characterize the effects of regulatory variants in different tissues [84], will be integrated into the VEP reference data in order to have the most current data available to the VEP for analysis.

As discussed above, standardized SO terms are used for describing variant consequences and VEP results can be output in VCF format. Work is ongoing to develop a comprehensive variant annotation data exchange format within the GA4GH. Furthermore, the GA4GH is defining standards for representation of associations between variants and phenotypes, traits, and diseases. The VEP will support such formats when they are mature.

Current annotation tools, including the VEP, annotate each input variant independently, without considering the potential compound effects of combining alternate alleles across multiple variant loci. This limitation means that having two or more variants affecting the same codon, or a shift in reading frame being corrected by a downstream variant, will not be taken into consideration. In future, given genotype data phased into haplotypes, the VEP will accurately annotate such events.

The VEP is also regularly extended and improved (see release notes at [85]) with new features added to both the core VEP code and the plugin library. Although these developments are frequently driven by new annotations or datasets available for H. sapiens, they are all designed to be compatible with any species. Once additional annotation and sequencing data are available in other species, the VEP extensions can be fully exploited for these too (e.g., 1000 Bulls project, the 1000 Chicken project, the 1001 Arabidopsis project, and the Functional Annotation of ANimal Genomes (FAANG) consortium). To improve genome-wide analysis, the VEP will leverage data from future sequencing projects, implement new algorithms and adopt data exchange standards and, therefore, bring continual benefit to variant interpretation.

The VEP algorithms and code are part of the freely available Ensembl API, coded in the Perl programming language. Time-critical components are written in C and incorporated into the API using the XS framework [86]. Installation of the VEP script triggers automated installation of the Ensembl API, along with the BioPerl API [87] upon which the Ensembl API depends. All interfaces to the VEP use the same underlying API calls, ensuring consistency across the different VEP access platforms when version control is observed.

To process the input data, sequential contiguous blocks of variants (default block size 5000) are read into an input memory buffer. Each variant is converted into an Ensembl VariationFeature object that represents a genomic location and a set of alleles. Variants in tab-delimited and Pileup formats are converted directly to objects; those in HGVS notation are resolved to their genomic coordinates by extracting the relevant reference feature (transcript, protein, or chromosome) using the Ensembl API. VCF input undergoes pre-processing to account for differences in how VCF and Ensembl represent unbalanced substitutions and indels. When using VEP’s forking functionality, the input buffer is divided amongst a number of sub-processes. Each sub-process carries out the analysis described hence and then the results are rejoined and sorted back into input order before being written to output.

Normalization of insertions and deletions in repetitive sequence and decomposition of complex variants are recommended as part of a robust pipeline to ensure consistency of annotation across datasets. Optionally, in a process analogous to that described in [88], VEP’s parser can be forced to decompose alternative alleles in complex variant descriptions to their most minimal representation by stripping identical bases from the 5′ and 3′ ends of the reference and alternative allele. This is not done by default as it may change the input position and allele string provided. Similarly, although it is a recommendation of the VCF format, the VEP does not left-normalize insertion or deletion variants in repetitive sequence. Enforcing this by default would cause discrepancies in input and output coordinates and also for HGVS nomenclature, whose coordinates must be right-normalized with respect to the transcript sequence. Tools such as vt [88] can be used to pre-process VCF input before use in VEP.

Input variants pass through a configurable quality-control process that checks for irregularities and inconsistencies. Variants that fail are reported via standard error output and/or in a warnings file. Checks include, for example, that allele lengths match input coordinates and the input reference allele matches that recorded in the reference genome.

The genomic loci overlapped by the variants in the input buffer are resolved to distinct megabase-sized regions. Each region corresponds to a single file on disk in the VEP cache, which contains objects serialized using Perl’s Storable framework [89]. For each region, the transcripts, regulatory features, and known variants are loaded from disk, deserialized into objects, and cached in memory. This avoids rereading from disk when the same region is overlapped by variants in consecutive input buffers. The publicly available Ensembl databases can be used in place of the cache files to avoid downloading the data in advance, though doing so incurs a performance penalty due to network transfer rates.

Transcripts have a configurable flank (default 5000 base pairs) to allow the VEP to assign upstream and downstream status to variants within the region flanking a transcript. A hash-based tree structure is used to search for overlaps between input variants and genomic features. For each overlap, a VariationFeatureOverlap object is created, with specific sub-classes for each genomic feature type: TranscriptVariation, RegulatoryFeatureVariation, MotifFeatureVariation. Each VariationFeatureOverlap object has two or more child VariationFeatureOverlapAllele objects representing each allele of the input variant—one representing the reference allele and one or more representing each of the alternative or mutant alleles. These objects are also sub-classed, with, for example, a TranscriptVariationAllele representing one allele of a variant overlapping a Transcript object.

For each TranscriptVariationAllele object, the API evaluates consequence types using a set of predicate functions. These assess whether, for example, a variant is predicted to cause a change in protein coding sequence (e.g., missense_variant). Prior to this, a series of pre-predicate checks are performed to improve runtime; for example, a variant does not need to be assessed for change to the protein sequence if it falls entirely within the intron of a transcript. These pre-predicate checks are also cached at each object “level”; for example, the location of a variant relative to the transcript structure is fixed at the TranscriptVariation level but the allele type can be different for each TranscriptVariationAllele. The pre-predicate checks improve runtime by a factor of around two on a typical resequencing-based input file. Without them, runtime is proportional to nfp, where n is the number of input variant alleles, f is the number of overlapped features, and p is the number of predicates; depending on a number of factors this can become as low as nfp/2 with pre-predicate checks enabled.

Predicates also make extensive use of caching: UTR, coding, and translated sequences are all cached on the Transcript object with intron structure and other frequently accessed data. Established components of the Ensembl API handle tasks such as splicing exons and re-translating mutated sequences. Alternative codon tables are used as appropriate for mitochondrial sequences and selenocysteines. If a predicate is true for a given TranscriptVariationAllele, an OverlapConsequence object is assigned representing the consequence type; this object contains the appropriate SO term along with synonyms and ranking information. Each OverlapConsequence object type corresponds to one predicate. Hierarchy in the predicate system preserves the tree structure of the SO such that only the most specific term that applies under any given parent term is assigned; this same tree structure allows for ontological-style querying and filtering of the results. Multiple OverlapConsequence objects may be added to a single VariationFeatureOverlapAllele or TranscriptVariationAllele object to allow for complex cases, such as a variant that falls in a splice-relevant region that also affects the coding sequence of the transcript.

HGVS notations are also derived from TranscriptVariationAlleles, though they undergo significant additional processing to conform to the nomenclature definition [90]. For example, insertions or deletions with respect to the transcript sequence must be reported at the most 3′ position possible when they fall in repetitive sequence.

VariationFeatureOverlapAllele objects are then converted for writing to output, a process that involves several extra stages. VariationFeatureOverlapAlleles can be filtered in various ways which can be configured, for example: reporting only one VariationFeatureOverlapAllele per input variant; removing intergenic VariationFeatureOverlapAlleles (i.e., those produced from variants that don’t overlap a genomic feature); filtering based on allele frequency of a co-located known variant. Additional data fields are retrieved at this stage from relevant objects, for example: external identifiers for transcripts (UniProt, CCDS); exon and intron numbers; clinical significance for co-located variants. It is also at this stage that any configured plugins are executed. They are passed the VariationFeatureOverlapAllele object, which has accessor methods for other objects, e.g., the Transcript, VariationFeature, or genomic Slice. As plugin modules are executed after the VEP consequence calculation, they have access to the VEP and Ensembl API objects before output data are written and return a data structure that is incorporated alongside the VEP’s main output data structure. The output data structure is then written to disk as one of several formats (tab-delimited, VCF, GVF, JSON), with the fields for each data format configurable at runtime. Output files contain headers describing the format and content of data fields, as well as version information for resources used.