Contents

0.1 Load Data

Below we compare the results of annotating variants via VariantAnnotation specifically repeating a subset of the steps in vignette Introduction to VariantAnnotation. We compare using data from 1,000 Genomes Phase 1 Variants:

0.1.1 VCF Data

First we read in the data from the VCF file:

library(VariantAnnotation)
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
vcf <- renameSeqlevels(vcf, c("22"="chr22"))
vcf
## 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, C...
## info(header(vcf)):
##              Number Type    Description                                   
##    LDAF      1      Float   MLE Allele Frequency Accounting for LD        
##    AVGPOST   1      Float   Average posterior probability from MaCH/Thu...
##    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 imprecis...
##    CIPOS     2      Integer Confidence interval around POS for imprecis...
##    END       1      Integer End position of the variant described in th...
##    HOMLEN    .      Integer Length of base pair identical micro-homolog...
##    HOMSEQ    .      String  Sequence of base pair identical micro-homol...
##    SVLEN     1      Integer Difference in length between REF and ALT al...
##    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...
##    AF        1      Float   Global Allele Frequency based on AC/AN        
##    AMR_AF    1      Float   Allele Frequency for samples from AMR based...
##    ASN_AF    1      Float   Allele Frequency for samples from ASN based...
##    AFR_AF    1      Float   Allele Frequency for samples from AFR based...
##    EUR_AF    1      Float   Allele Frequency for samples from EUR based...
##    VT        1      String  indicates what type of variant the line rep...
##    SNPSOURCE .      String  indicates if a snp was called when analysin...
## geno(vcf):
##   SimpleList 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 .      Float  Genotype Likelihoods

The file chr22.vcf.gz within package VariantAnnotation holds data for 5 of the 1,092 individuals in 1,000 Genomes, starting at position 50300078 and ending at position 50999964.

HG00096 HG00097 HG00099 HG00100 HG00101

0.1.2 Google Genomics Data

Important data differences to note:

library(GoogleGenomics)
# This vignette is authenticated on package load from the env variable GOOGLE_API_KEY.
# When running interactively, call the authenticate method.
# ?authenticate
# We're just getting the first few variants so that this runs quickly.
# If we wanted to get them all, we sould set end=50999964.
granges <- getVariants(datasetId="10473108253681171589",
                       chromosome="22",
                       start=50300077,
                       end=50303000,
                       converter=variantsToGRanges)
## Fetching variants page.
## Continuing variant query with the nextPageToken: CIWK_hcQ_Nr-8ZGlk78i
## Fetching variants page.
## Continuing variant query with the nextPageToken: CIaM_hcQvM_un46SmcJ0
## Fetching variants page.
## Continuing variant query with the nextPageToken: CJGN_hcQiYna3cSwse_RAQ
## Fetching variants page.
## Continuing variant query with the nextPageToken: CMKQ_hcQr9OUpKu4n_xz
## Fetching variants page.
## Continuing variant query with the nextPageToken: CMKS_hcQn-nIneaph7WJAQ
## Fetching variants page.
## Continuing variant query with the nextPageToken: CKOU_hcQ2eCCwsmBysdd
## Fetching variants page.
## Continuing variant query with the nextPageToken: CLWV_hcQmo-n-tX08uPyAQ
## Fetching variants page.
## Continuing variant query with the nextPageToken: CMqX_hcQ__3H4KKV9tUR
## Fetching variants page.
## Continuing variant query with the nextPageToken: CKuZ_hcQ9pHXwa742bKgAQ
## Fetching variants page.
## Continuing variant query with the nextPageToken: CLqa_hcQ-6zjo_Dc86gy
## Fetching variants page.
## Continuing variant query with the nextPageToken: CLSb_hcQ3vrDo5XJwsB_
## Fetching variants page.
## Continuing variant query with the nextPageToken: CKOd_hcQ6IPgk-PVjskK
## Fetching variants page.
## Continuing variant query with the nextPageToken: CPGe_hcQ-8Saw6Ojg8yXAQ
## Fetching variants page.
## Continuing variant query with the nextPageToken: CPGf_hcQuKDdodiCvLUm
## Fetching variants page.
## Variants are now available.

0.1.3 Compare the Loaded Data

Ensure that the data retrieved by each matches:

vcf <- vcf[1:length(granges)] # Truncate the VCF data so that it is the same 
                              # set as what was retrieved from the API.
library(testthat)
expect_equal(start(granges), start(vcf))
expect_equal(end(granges), end(vcf))
expect_equal(as.character(granges$REF), as.character(ref(vcf)))
expect_equal(as.character(unlist(granges$ALT)), as.character(unlist(alt(vcf))))
expect_equal(granges$QUAL, qual(vcf))
expect_equal(granges$FILTER, filt(vcf))

0.2 Compare the Annotation Results

Now locate the protein coding variants in each:

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

rd <- rowRanges(vcf)
vcf_locations <- locateVariants(rd, txdb, CodingVariants())
## 'select()' returned many:1 mapping between keys and columns
vcf_locations
## GRanges object with 7 ranges and 9 metadata columns:
##     seqnames               ranges strand | LOCATION  LOCSTART    LOCEND
##        <Rle>            <IRanges>  <Rle> | <factor> <integer> <integer>
##   1    chr22 [50301422, 50301422]      - |   coding       939       939
##   2    chr22 [50301476, 50301476]      - |   coding       885       885
##   3    chr22 [50301488, 50301488]      - |   coding       873       873
##   4    chr22 [50301494, 50301494]      - |   coding       867       867
##   5    chr22 [50301584, 50301584]      - |   coding       777       777
##   6    chr22 [50302962, 50302962]      - |   coding       698       698
##   7    chr22 [50302995, 50302995]      - |   coding       665       665
##       QUERYID        TXID         CDSID      GENEID       PRECEDEID
##     <integer> <character> <IntegerList> <character> <CharacterList>
##   1        24       75253        218562       79087                
##   2        25       75253        218562       79087                
##   3        26       75253        218562       79087                
##   4        27       75253        218562       79087                
##   5        28       75253        218562       79087                
##   6        57       75253        218563       79087                
##   7        58       75253        218563       79087                
##            FOLLOWID
##     <CharacterList>
##   1                
##   2                
##   3                
##   4                
##   5                
##   6                
##   7                
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
granges_locations <- locateVariants(granges, txdb, CodingVariants())
## 'select()' returned many:1 mapping between keys and columns
granges_locations
## GRanges object with 7 ranges and 9 metadata columns:
##     seqnames               ranges strand | LOCATION  LOCSTART    LOCEND
##        <Rle>            <IRanges>  <Rle> | <factor> <integer> <integer>
##   1    chr22 [50301422, 50301422]      - |   coding       939       939
##   2    chr22 [50301476, 50301476]      - |   coding       885       885
##   3    chr22 [50301488, 50301488]      - |   coding       873       873
##   4    chr22 [50301494, 50301494]      - |   coding       867       867
##   5    chr22 [50301584, 50301584]      - |   coding       777       777
##   6    chr22 [50302962, 50302962]      - |   coding       698       698
##   7    chr22 [50302995, 50302995]      - |   coding       665       665
##       QUERYID        TXID         CDSID      GENEID       PRECEDEID
##     <integer> <character> <IntegerList> <character> <CharacterList>
##   1        24       75253        218562       79087                
##   2        25       75253        218562       79087                
##   3        26       75253        218562       79087                
##   4        27       75253        218562       79087                
##   5        28       75253        218562       79087                
##   6        57       75253        218563       79087                
##   7        58       75253        218563       79087                
##            FOLLOWID
##     <CharacterList>
##   1                
##   2                
##   3                
##   4                
##   5                
##   6                
##   7                
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
expect_equal(granges_locations, vcf_locations)

And predict the effect of the protein coding variants:

library(BSgenome.Hsapiens.UCSC.hg19)
vcf_coding <- predictCoding(vcf, txdb, seqSource=Hsapiens)
vcf_coding
## GRanges object with 7 ranges and 17 metadata columns:
##                   seqnames               ranges strand | paramRangeID
##                      <Rle>            <IRanges>  <Rle> |     <factor>
##       rs114335781    chr22 [50301422, 50301422]      - |         <NA>
##         rs8135963    chr22 [50301476, 50301476]      - |         <NA>
##   22:50301488_C/T    chr22 [50301488, 50301488]      - |         <NA>
##   22:50301494_G/A    chr22 [50301494, 50301494]      - |         <NA>
##   22:50301584_C/T    chr22 [50301584, 50301584]      - |         <NA>
##       rs114264124    chr22 [50302962, 50302962]      - |         <NA>
##       rs149209714    chr22 [50302995, 50302995]      - |         <NA>
##                              REF                ALT      QUAL      FILTER
##                   <DNAStringSet> <DNAStringSetList> <numeric> <character>
##       rs114335781              G                  A       100        PASS
##         rs8135963              T                  C       100        PASS
##   22:50301488_C/T              C                  T       100        PASS
##   22:50301494_G/A              G                  A       100        PASS
##   22:50301584_C/T              C                  T       100        PASS
##       rs114264124              C                  T       100        PASS
##       rs149209714              C                  G       100        PASS
##                        varAllele     CDSLOC    PROTEINLOC   QUERYID
##                   <DNAStringSet>  <IRanges> <IntegerList> <integer>
##       rs114335781              T [939, 939]           313        24
##         rs8135963              G [885, 885]           295        25
##   22:50301488_C/T              A [873, 873]           291        26
##   22:50301494_G/A              T [867, 867]           289        27
##   22:50301584_C/T              A [777, 777]           259        28
##       rs114264124              A [698, 698]           233        57
##       rs149209714              C [665, 665]           222        58
##                          TXID         CDSID      GENEID   CONSEQUENCE
##                   <character> <IntegerList> <character>      <factor>
##       rs114335781       75253        218562       79087    synonymous
##         rs8135963       75253        218562       79087    synonymous
##   22:50301488_C/T       75253        218562       79087    synonymous
##   22:50301494_G/A       75253        218562       79087    synonymous
##   22:50301584_C/T       75253        218562       79087    synonymous
##       rs114264124       75253        218563       79087 nonsynonymous
##       rs149209714       75253        218563       79087 nonsynonymous
##                         REFCODON       VARCODON         REFAA
##                   <DNAStringSet> <DNAStringSet> <AAStringSet>
##       rs114335781            ATC            ATT             I
##         rs8135963            GCA            GCG             A
##   22:50301488_C/T            CCG            CCA             P
##   22:50301494_G/A            CAC            CAT             H
##   22:50301584_C/T            CCG            CCA             P
##       rs114264124            CGG            CAG             R
##       rs149209714            GGA            GCA             G
##                           VARAA
##                   <AAStringSet>
##       rs114335781             I
##         rs8135963             A
##   22:50301488_C/T             P
##   22:50301494_G/A             H
##   22:50301584_C/T             P
##       rs114264124             Q
##       rs149209714             A
##   -------
##   seqinfo: 1 sequence from hg19 genome; no seqlengths
granges_coding <- predictCoding(rep(granges, elementNROWS(granges$ALT)),
                                txdb,
                                seqSource=Hsapiens,
                                varAllele=unlist(granges$ALT, use.names=FALSE))

granges_coding
## GRanges object with 7 ranges and 16 metadata columns:
##               seqnames               ranges strand |            REF
##                  <Rle>            <IRanges>  <Rle> | <DNAStringSet>
##   rs114335781    chr22 [50301422, 50301422]      - |              G
##     rs8135963    chr22 [50301476, 50301476]      - |              T
##   rs200080075    chr22 [50301488, 50301488]      - |              C
##   rs147801200    chr22 [50301494, 50301494]      - |              G
##   rs138060012    chr22 [50301584, 50301584]      - |              C
##   rs114264124    chr22 [50302962, 50302962]      - |              C
##   rs149209714    chr22 [50302995, 50302995]      - |              C
##                              ALT      QUAL      FILTER      varAllele
##               <DNAStringSetList> <numeric> <character> <DNAStringSet>
##   rs114335781                  A       100        PASS              T
##     rs8135963                  C       100        PASS              G
##   rs200080075                  T       100        PASS              A
##   rs147801200                  A       100        PASS              T
##   rs138060012                  T       100        PASS              A
##   rs114264124                  T       100        PASS              A
##   rs149209714                  G       100        PASS              C
##                   CDSLOC    PROTEINLOC   QUERYID        TXID         CDSID
##                <IRanges> <IntegerList> <integer> <character> <IntegerList>
##   rs114335781 [939, 939]           313        24       75253        218562
##     rs8135963 [885, 885]           295        25       75253        218562
##   rs200080075 [873, 873]           291        26       75253        218562
##   rs147801200 [867, 867]           289        27       75253        218562
##   rs138060012 [777, 777]           259        28       75253        218562
##   rs114264124 [698, 698]           233        57       75253        218563
##   rs149209714 [665, 665]           222        58       75253        218563
##                    GENEID   CONSEQUENCE       REFCODON       VARCODON
##               <character>      <factor> <DNAStringSet> <DNAStringSet>
##   rs114335781       79087    synonymous            ATC            ATT
##     rs8135963       79087    synonymous            GCA            GCG
##   rs200080075       79087    synonymous            CCG            CCA
##   rs147801200       79087    synonymous            CAC            CAT
##   rs138060012       79087    synonymous            CCG            CCA
##   rs114264124       79087 nonsynonymous            CGG            CAG
##   rs149209714       79087 nonsynonymous            GGA            GCA
##                       REFAA         VARAA
##               <AAStringSet> <AAStringSet>
##   rs114335781             I             I
##     rs8135963             A             A
##   rs200080075             P             P
##   rs147801200             H             H
##   rs138060012             P             P
##   rs114264124             R             Q
##   rs149209714             G             A
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths
expect_equal(as.matrix(granges_coding$REFCODON), as.matrix(vcf_coding$REFCODON))
expect_equal(as.matrix(granges_coding$VARCODON), as.matrix(vcf_coding$VARCODON))
expect_equal(granges_coding$GENEID, vcf_coding$GENEID)
expect_equal(granges_coding$CONSEQUENCE, vcf_coding$CONSEQUENCE)

0.3 Provenance

sessionInfo()
## R version 3.3.0 (2016-05-03)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.4 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] testthat_1.0.2                         
##  [2] ggbio_1.20.1                           
##  [3] ggplot2_2.1.0                          
##  [4] org.Hs.eg.db_3.3.0                     
##  [5] BSgenome.Hsapiens.UCSC.hg19_1.4.0      
##  [6] BSgenome_1.40.0                        
##  [7] rtracklayer_1.32.0                     
##  [8] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
##  [9] GenomicFeatures_1.24.2                 
## [10] AnnotationDbi_1.34.2                   
## [11] GoogleGenomics_1.4.2                   
## [12] VariantAnnotation_1.18.1               
## [13] GenomicAlignments_1.8.0                
## [14] Rsamtools_1.24.0                       
## [15] Biostrings_2.40.0                      
## [16] XVector_0.12.0                         
## [17] SummarizedExperiment_1.2.2             
## [18] Biobase_2.32.0                         
## [19] GenomicRanges_1.24.0                   
## [20] GenomeInfoDb_1.8.2                     
## [21] IRanges_2.6.0                          
## [22] S4Vectors_0.10.0                       
## [23] BiocGenerics_0.18.0                    
## [24] knitr_1.13                             
## [25] BiocStyle_2.0.2                        
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.1.0                    jsonlite_0.9.20              
##  [3] AnnotationHub_2.4.2           splines_3.3.0                
##  [5] Formula_1.2-1                 shiny_0.13.2                 
##  [7] interactiveDisplayBase_1.10.3 latticeExtra_0.6-28          
##  [9] RBGL_1.48.0                   yaml_2.1.13                  
## [11] RSQLite_1.0.0                 lattice_0.20-33              
## [13] biovizBase_1.20.0             chron_2.3-47                 
## [15] digest_0.6.9                  RColorBrewer_1.1-2           
## [17] colorspace_1.2-6              htmltools_0.3.5              
## [19] httpuv_1.3.3                  Matrix_1.2-6                 
## [21] plyr_1.8.3                    OrganismDbi_1.14.0           
## [23] XML_3.98-1.4                  biomaRt_2.28.0               
## [25] zlibbioc_1.18.0               xtable_1.8-2                 
## [27] scales_0.4.0                  BiocParallel_1.6.2           
## [29] nnet_7.3-12                   crayon_1.3.1                 
## [31] survival_2.39-4               magrittr_1.5                 
## [33] mime_0.4                      memoise_1.0.0                
## [35] evaluate_0.9                  GGally_1.0.1                 
## [37] foreign_0.8-66                graph_1.50.0                 
## [39] BiocInstaller_1.22.2          tools_3.3.0                  
## [41] data.table_1.9.6              formatR_1.4                  
## [43] stringr_1.0.0                 munsell_0.4.3                
## [45] cluster_2.0.4                 ensembldb_1.4.2              
## [47] grid_3.3.0                    RCurl_1.95-4.8               
## [49] dichromat_2.0-0               rjson_0.2.15                 
## [51] labeling_0.3                  bitops_1.0-6                 
## [53] rmarkdown_0.9.6               gtable_0.2.0                 
## [55] DBI_0.4-1                     reshape_0.8.5                
## [57] curl_0.9.7                    reshape2_1.4.1               
## [59] R6_2.1.2                      gridExtra_2.2.1              
## [61] Hmisc_3.17-4                  stringi_1.0-1                
## [63] Rcpp_0.12.5                   rpart_4.1-10                 
## [65] acepack_1.3-3.3