# ensemblVEP: using the REST API with Bioconductor

# Introduction

Ensembl’s Variant Effect Predictor is described in McLaren et al. (2016).

Prior to Bioconductor 3.19, the ensemblVEP package provided access to Ensembl’s predictions through an interface between Perl and MySQL.

In 3.19 VariantAnnotation supports the use of the VEP component of the REST API at https://rest.ensembl.org.

# Acquire annotation on variants from a VCF file

The function `vep_by_region` will accept a VCF object as defined in VariantAnnotation.

    library(VariantAnnotation)
    fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
    r22 = readVcf(fl)
    r22

    ## class: CollapsedVCF
    ## dim: 10376 5
    ## rowRanges(vcf):
    ##   GRanges with 5 metadata columns: paramRangeID, REF, ALT, QUAL, FILTER
    ## info(vcf):
    ##   DataFrame with 22 columns: LDAF, AVGPOST, RSQ, ERATE, THETA, CIEND, CIPOS,...
    ## info(header(vcf)):
    ##              Number Type    Description
    ##    LDAF      1      Float   MLE Allele Frequency Accounting for LD
    ##    AVGPOST   1      Float   Average posterior probability from MaCH/Thunder
    ##    RSQ       1      Float   Genotype imputation quality from MaCH/Thunder
    ##    ERATE     1      Float   Per-marker Mutation rate from MaCH/Thunder
    ##    THETA     1      Float   Per-marker Transition rate from MaCH/Thunder
    ##    CIEND     2      Integer Confidence interval around END for imprecise var...
    ##    CIPOS     2      Integer Confidence interval around POS for imprecise var...
    ##    END       1      Integer End position of the variant described in this re...
    ##    HOMLEN    .      Integer Length of base pair identical micro-homology at ...
    ##    HOMSEQ    .      String  Sequence of base pair identical micro-homology a...
    ##    SVLEN     1      Integer Difference in length between REF and ALT alleles
    ##    SVTYPE    1      String  Type of structural variant
    ##    AC        .      Integer Alternate Allele Count
    ##    AN        1      Integer Total Allele Count
    ##    AA        1      String  Ancestral Allele, ftp://ftp.1000genomes.ebi.ac.u...
    ##    AF        1      Float   Global Allele Frequency based on AC/AN
    ##    AMR_AF    1      Float   Allele Frequency for samples from AMR based on A...
    ##    ASN_AF    1      Float   Allele Frequency for samples from ASN based on A...
    ##    AFR_AF    1      Float   Allele Frequency for samples from AFR based on A...
    ##    EUR_AF    1      Float   Allele Frequency for samples from EUR based on A...
    ##    VT        1      String  indicates what type of variant the line represents
    ##    SNPSOURCE .      String  indicates if a snp was called when analysing the...
    ## geno(vcf):
    ##   List of length 3: GT, DS, GL
    ## geno(header(vcf)):
    ##       Number Type   Description
    ##    GT 1      String Genotype
    ##    DS 1      Float  Genotype dosage from MaCH/Thunder
    ##    GL G      Float  Genotype Likelihoods

In this example we confine attention to single nucleotide variants.

There is a limit of 200 locations in a request, and 55000 requests per hour. We’ll base our query on 100 positions in the chr22 VCF.

    dr = which(width(rowRanges(r22))!=1)
    r22s = r22[-dr]
    res = vep_by_region(r22[1:100], snv_only=FALSE, chk_max=FALSE)
    jans = toJSON(content(res))

There are various ways to work with the result of this query to the API. We’ll use the rjsoncons JSON processing infrastructure to dig in and understand aspects of the API behavior.

First, the top-level concepts produced for each variant can be retrieved using

    library(rjsoncons)
    names(jsonlite::fromJSON(jmespath(jans, "[*]")))

    ##  [1] "end"                             "allele_string"
    ##  [3] "id"                              "seq_region_name"
    ##  [5] "strand"                          "start"
    ##  [7] "transcript_consequences"         "most_severe_consequence"
    ##  [9] "input"                           "assembly_name"
    ## [11] "colocated_variants"              "motif_feature_consequences"
    ## [13] "regulatory_feature_consequences"

Annotation of the most severe consequence known will typically be of interest:

    table(jsonlite::fromJSON(jmespath(jans, "[*].most_severe_consequence")))

    ##
    ##   5_prime_UTR_variant        intron_variant splice_region_variant
    ##                    22                    76                     2

There is variability in the structure of data returned for each query.

    head(fromJSON(jmespath(jans, "[*].regulatory_feature_consequences")))

    ## [[1]]
    ##        biotype regulatory_feature_id variant_allele consequence_terms   impact
    ## 1 CTCF_bin....          ENSR22_5....              T      regulato.... MODIFIER

Furthermore, the content of the motif feature consequences field seems very peculiar.

    table(unlist(fromJSON(jmespath(jans, "[*].motif_feature_consequences"))))

    ##
    ##                  -0.018                      -1                       0
    ##                       1                       2                       4
    ##                       1                      13                       2
    ##                       3                       1                       1
    ##                       4                       5                       7
    ##                       1                       1                       1
    ##                       A         ENSM00000606562         ENSM00000606563
    ##                       1                       1                       1
    ##         ENSM00000606566         ENSM00000606569              ENSPFM0109
    ##                       2                       1                       1
    ##              ENSPFM0290              ENSPFM0357                       G
    ##                       2                       2                       4
    ##                   HNF4A                MODIFIER                       N
    ##                       2                       5                       3
    ##                   NR2F1                    RXRA                     SP1
    ##                       2                       2                       2
    ##                   TBX21 TF_binding_site_variant                       Y
    ##                       1                       5                       2

# Transforming the API response to GRanges

We’ll consider the following approach to converting the API response to a GenomicRanges GRanges instance. Eventually this may become part of the package.

    library(GenomicRanges)
    .make_GRanges = function( vep_response ) {
      stopifnot(inherits(vep_response, "response"))  # httr
      nested = fromJSON(toJSON(content(vep_response)))
      ini = GRanges(seqnames = unlist(nested$seq_region_name),
        IRanges(start=unlist(nested$start), end=unlist(nested$end)))
      dr = match(c("seq_region_name", "start", "end"), names(nested))
      mcols(ini) = DataFrame(nested[,-dr])
      ini
    }
    tstg = .make_GRanges( res )
    tstg[,1]  # full print is unwieldy

    ## GRanges object with 100 ranges and 1 metadata column:
    ##         seqnames    ranges strand | allele_string
    ##            <Rle> <IRanges>  <Rle> |        <list>
    ##     [1]       22  50300078      * |           A/G
    ##     [2]       22  50300086      * |           C/T
    ##     [3]       22  50300101      * |           G/A
    ##     [4]       22  50300113      * |           C/T
    ##     [5]       22  50300166      * |           C/T
    ##     ...      ...       ...    ... .           ...
    ##    [96]       22  50304748      * |           G/A
    ##    [97]       22  50304805      * |           G/A
    ##    [98]       22  50304935      * |           C/A
    ##    [99]       22  50304943      * |           G/C
    ##   [100]       22  50305084      * |           T/C
    ##   -------
    ##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

    names(mcols(tstg))

    ##  [1] "allele_string"                   "id"
    ##  [3] "strand"                          "transcript_consequences"
    ##  [5] "most_severe_consequence"         "input"
    ##  [7] "assembly_name"                   "colocated_variants"
    ##  [9] "motif_feature_consequences"      "regulatory_feature_consequences"

Now information about variants can be retrieved with range operations. Deep annotation requires nested structure of the metadata columns.

    mcols(tstg)[1, "transcript_consequences"]

    ## [[1]]
    ##    variant_allele   hgnc_id      biotype transcript_id gene_symbol_source
    ## 1               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 2               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 3               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 4               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 5               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 6               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 7               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 8               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 9               G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 10              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 11              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 12              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 13              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 14              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 15              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 16              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 17              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 18              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 19              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 20              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 21              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 22              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 23              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 24              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 25              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 26              G HGNC:9104 protein_....  ENST0000....               HGNC
    ## 27              G HGNC:9104 protein_....  ENST0000....               HGNC
    ##    gene_symbol consequence_terms      gene_id strand   impact      flags
    ## 1       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 2       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER cds_end_NF
    ## 3       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER cds_end_NF
    ## 4       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 5       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 6       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 7       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 8       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 9       PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 10      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 11      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 12      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 13      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 14      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 15      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 16      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 17      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 18      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 19      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 20      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 21      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 22      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 23      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 24      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 25      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 26      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER
    ## 27      PLXNB2      intron_v.... ENSG0000....     -1 MODIFIER

# Further work

An important element of prior work in ensemblVEP supports feeding annotation back into the VCF used to generate the effect prediction query. This seems feasible but concrete use cases are of interest.

# References

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