Plotting Alignments

0.1 Working with Reads

Google Genomics implements the GA4GH reads API and this R package can retrieve data from that implementation. For more detail, see https://cloud.google.com/genomics/v1beta2/reference/reads

library(GoogleGenomics)
# This vignette is authenticated on package load from the env variable GOOGLE_API_KEY.
# When running interactively, call the authenticate method.
# ?authenticate

By default, this function retrieves reads for a small genomic region for one sample in 1,000 Genomes.

reads <- getReads()

## Fetching reads page.

## Reads are now available.

length(reads)

## [1] 18

We can see that 18 individual reads were returned and that the JSON response was deserialized into an R list object:

class(reads)

## [1] "list"

mode(reads)

## [1] "list"

The top level field names are:

names(reads[[1]])

##  [1] "id"                        "readGroupId"              
##  [3] "readGroupSetId"            "fragmentName"             
##  [5] "properPlacement"           "duplicateFragment"        
##  [7] "fragmentLength"            "readNumber"               
##  [9] "numberReads"               "failedVendorQualityChecks"
## [11] "alignment"                 "secondaryAlignment"       
## [13] "supplementaryAlignment"    "alignedSequence"          
## [15] "alignedQuality"            "nextMatePosition"         
## [17] "info"

And examining the alignment we see:

reads[[1]]$alignedSequence

## [1] "AACAAAAAACTGTCTAACAAGATTTTATGGTTTATAGACCATGATTCCCCGGACACATTAGATAGAAATCTGGGCAAGAGAAGAAAAAAAGGTCAGAGTT"

reads[[1]]$alignment$position$referenceName

## [1] "22"

reads[[1]]$alignment$position$position

## [1] "16051308"

This is good, but this data becomes much more useful when it is converted to Bioconductor data types. For example, we can convert reads in this list representation to GAlignments:

readsToGAlignments(reads)

## GAlignments object with 18 alignments and 1 metadata column:
##                      seqnames strand       cigar    qwidth     start
##                         <Rle>  <Rle> <character> <integer> <integer>
##   ERR251039.44923356    chr22      -        100M       100  16051309
##   ERR251039.28556355    chr22      +        100M       100  16051323
##   ERR016214.27010769    chr22      +     75M1D6M        81  16051330
##   ERR016213.15718767    chr22      +         68M        68  16051400
##   ERR016213.29749886    chr22      -         68M        68  16051403
##                  ...      ...    ...         ...       ...       ...
##    ERR251040.1475552    chr22      +        100M       100  16051454
##    ERR251040.3762117    chr22      +        100M       100  16051469
##   ERR251040.11853979    chr22      +        100M       100  16051478
##   ERR251040.34469189    chr22      +        100M       100  16051486
##   ERR251040.38772006    chr22      -        100M       100  16051498
##                            end     width     njunc |      flag
##                      <integer> <integer> <integer> | <numeric>
##   ERR251039.44923356  16051408       100         0 |       147
##   ERR251039.28556355  16051422       100         0 |       163
##   ERR016214.27010769  16051411        82         0 |        99
##   ERR016213.15718767  16051467        68         0 |       163
##   ERR016213.29749886  16051470        68         0 |       153
##                  ...       ...       ...       ... .       ...
##    ERR251040.1475552  16051553       100         0 |       163
##    ERR251040.3762117  16051568       100         0 |       163
##   ERR251040.11853979  16051577       100         0 |       163
##   ERR251040.34469189  16051585       100         0 |        99
##   ERR251040.38772006  16051597       100         0 |        83
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

0.2 Plotting Alignments

Let’s take a look at the reads that overlap rs9536314 for sample NA12893 within the Illumina Platinum Genomes dataset. This SNP resides on chromosome 13 at position 33628137 in 0-based coordinates.

# Change the values of 'chromosome', 'start', or 'end' here if you wish to plot 
# alignments from a different portion of the genome.
alignments <- getReads(readGroupSetId="CMvnhpKTFhD3he72j4KZuyc",
                       chromosome="chr13",
                       start=33628130,
                       end=33628145,
                       converter=readsToGAlignments)

## Fetching reads page.

## Reads are now available.

alignments

## GAlignments object with 67 alignments and 1 metadata column:
##                                             seqnames strand       cigar
##                                                <Rle>  <Rle> <character>
##    HSQ1008:141:D0CC8ACXX:4:2201:3978:155671    chr13      +        101M
##    HSQ1008:141:D0CC8ACXX:3:1206:7108:129302    chr13      -        101M
##    HSQ1008:141:D0CC8ACXX:2:1307:18541:80804    chr13      +        101M
##    HSQ1008:141:D0CC8ACXX:4:2108:12961:30710    chr13      -        101M
##     HSQ1008:141:D0CC8ACXX:1:2308:18484:6721    chr13      -        101M
##                                         ...      ...    ...         ...
##   HSQ1008:141:D0CC8ACXX:2:1308:19101:126736    chr13      +        101M
##    HSQ1008:141:D0CC8ACXX:2:1302:13881:79089    chr13      -        101M
##     HSQ1008:141:D0CC8ACXX:2:2108:1700:58641    chr13      -        101M
##     HSQ1008:141:D0CC8ACXX:1:2104:4131:31083    chr13      +        101M
##     HSQ1008:141:D0CC8ACXX:4:2302:3851:18256    chr13      +        101M
##                                                qwidth     start       end
##                                             <integer> <integer> <integer>
##    HSQ1008:141:D0CC8ACXX:4:2201:3978:155671       101  33628031  33628131
##    HSQ1008:141:D0CC8ACXX:3:1206:7108:129302       101  33628033  33628133
##    HSQ1008:141:D0CC8ACXX:2:1307:18541:80804       101  33628039  33628139
##    HSQ1008:141:D0CC8ACXX:4:2108:12961:30710       101  33628042  33628142
##     HSQ1008:141:D0CC8ACXX:1:2308:18484:6721       101  33628043  33628143
##                                         ...       ...       ...       ...
##   HSQ1008:141:D0CC8ACXX:2:1308:19101:126736       101  33628133  33628233
##    HSQ1008:141:D0CC8ACXX:2:1302:13881:79089       101  33628136  33628236
##     HSQ1008:141:D0CC8ACXX:2:2108:1700:58641       101  33628136  33628236
##     HSQ1008:141:D0CC8ACXX:1:2104:4131:31083       101  33628139  33628239
##     HSQ1008:141:D0CC8ACXX:4:2302:3851:18256       101  33628141  33628241
##                                                 width     njunc |
##                                             <integer> <integer> |
##    HSQ1008:141:D0CC8ACXX:4:2201:3978:155671       101         0 |
##    HSQ1008:141:D0CC8ACXX:3:1206:7108:129302       101         0 |
##    HSQ1008:141:D0CC8ACXX:2:1307:18541:80804       101         0 |
##    HSQ1008:141:D0CC8ACXX:4:2108:12961:30710       101         0 |
##     HSQ1008:141:D0CC8ACXX:1:2308:18484:6721       101         0 |
##                                         ...       ...       ... .
##   HSQ1008:141:D0CC8ACXX:2:1308:19101:126736       101         0 |
##    HSQ1008:141:D0CC8ACXX:2:1302:13881:79089       101         0 |
##     HSQ1008:141:D0CC8ACXX:2:2108:1700:58641       101         0 |
##     HSQ1008:141:D0CC8ACXX:1:2104:4131:31083       101         0 |
##     HSQ1008:141:D0CC8ACXX:4:2302:3851:18256       101         0 |
##                                                  flag
##                                             <numeric>
##    HSQ1008:141:D0CC8ACXX:4:2201:3978:155671       163
##    HSQ1008:141:D0CC8ACXX:3:1206:7108:129302        83
##    HSQ1008:141:D0CC8ACXX:2:1307:18541:80804        99
##    HSQ1008:141:D0CC8ACXX:4:2108:12961:30710        83
##     HSQ1008:141:D0CC8ACXX:1:2308:18484:6721       147
##                                         ...       ...
##   HSQ1008:141:D0CC8ACXX:2:1308:19101:126736        99
##    HSQ1008:141:D0CC8ACXX:2:1302:13881:79089        83
##     HSQ1008:141:D0CC8ACXX:2:2108:1700:58641       147
##     HSQ1008:141:D0CC8ACXX:1:2104:4131:31083        99
##     HSQ1008:141:D0CC8ACXX:4:2302:3851:18256        99
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Notice that we passed the converter to the getReads method so that we’re immediately working with GAlignments which means that we can start taking advantage of other Bioconductor functionality. Also keep in mind that the parameters start and end are expressed in 0-based coordinates per the GA4GH specification but the Bioconductor data type converters in GoogleGenomics, by default, transform the returned data to 1-based coordinates.

library(ggbio)

We can display the basic alignments and coverage data:

alignmentPlot <- autoplot(alignments, aes(color=strand, fill=strand))

## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

## Scale for 'fill' is already present. Adding another scale for 'fill',
## which will replace the existing scale.

coveragePlot <- ggplot(as(alignments, "GRanges")) +
                stat_coverage(color="gray40", fill="skyblue")
tracks(alignmentPlot, coveragePlot,
       xlab="Reads overlapping rs9536314 for NA12893")

And also display the ideogram for the corresponding location on the chromosome:

ideogramPlot <- plotIdeogram(genome="hg19", subchr="chr13")

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

ideogramPlot + xlim(as(alignments, "GRanges"))

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

0.3 Provenance

sessionInfo()

## R version 3.4.1 (2017-06-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
## 
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.5-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.5-bioc/R/lib/libRlapack.so
## 
## 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] ggbio_1.24.1                           
##  [2] ggplot2_2.2.1                          
##  [3] org.Hs.eg.db_3.4.1                     
##  [4] BSgenome.Hsapiens.UCSC.hg19_1.4.0      
##  [5] BSgenome_1.44.1                        
##  [6] rtracklayer_1.36.4                     
##  [7] TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
##  [8] GenomicFeatures_1.28.5                 
##  [9] AnnotationDbi_1.38.2                   
## [10] GoogleGenomics_1.8.1                   
## [11] VariantAnnotation_1.22.3               
## [12] GenomicAlignments_1.12.2               
## [13] Rsamtools_1.28.0                       
## [14] Biostrings_2.44.2                      
## [15] XVector_0.16.0                         
## [16] SummarizedExperiment_1.6.4             
## [17] DelayedArray_0.2.7                     
## [18] matrixStats_0.52.2                     
## [19] Biobase_2.36.2                         
## [20] GenomicRanges_1.28.5                   
## [21] GenomeInfoDb_1.12.2                    
## [22] IRanges_2.10.3                         
## [23] S4Vectors_0.14.4                       
## [24] BiocGenerics_0.22.0                    
## [25] knitr_1.17                             
## [26] BiocStyle_2.4.1                        
## 
## loaded via a namespace (and not attached):
##  [1] ProtGenerics_1.8.0            bitops_1.0-6                 
##  [3] bit64_0.9-7                   RColorBrewer_1.1-2           
##  [5] httr_1.3.1                    rprojroot_1.2                
##  [7] tools_3.4.1                   backports_1.1.0              
##  [9] R6_2.2.2                      rpart_4.1-11                 
## [11] Hmisc_4.0-3                   DBI_0.7                      
## [13] lazyeval_0.2.0                colorspace_1.3-2             
## [15] nnet_7.3-12                   gridExtra_2.3                
## [17] GGally_1.3.2                  bit_1.1-12                   
## [19] curl_2.8.1                    compiler_3.4.1               
## [21] graph_1.54.0                  htmlTable_1.9                
## [23] labeling_0.3                  scales_0.5.0                 
## [25] checkmate_1.8.3               RBGL_1.52.0                  
## [27] stringr_1.2.0                 digest_0.6.12                
## [29] foreign_0.8-69                rmarkdown_1.6                
## [31] base64enc_0.1-3               dichromat_2.0-0              
## [33] pkgconfig_2.0.1               htmltools_0.3.6              
## [35] ensembldb_2.0.4               htmlwidgets_0.9              
## [37] rlang_0.1.2                   RSQLite_2.0                  
## [39] BiocInstaller_1.26.1          shiny_1.0.5                  
## [41] jsonlite_1.5                  BiocParallel_1.10.1          
## [43] acepack_1.4.1                 RCurl_1.95-4.8               
## [45] magrittr_1.5                  GenomeInfoDbData_0.99.0      
## [47] Formula_1.2-2                 Matrix_1.2-11                
## [49] Rcpp_0.12.12                  munsell_0.4.3                
## [51] stringi_1.1.5                 yaml_2.1.14                  
## [53] zlibbioc_1.22.0               plyr_1.8.4                   
## [55] AnnotationHub_2.8.2           grid_3.4.1                   
## [57] blob_1.1.0                    lattice_0.20-35              
## [59] splines_3.4.1                 rjson_0.2.15                 
## [61] reshape2_1.4.2                biomaRt_2.32.1               
## [63] XML_3.98-1.9                  evaluate_0.10.1              
## [65] biovizBase_1.24.0             latticeExtra_0.6-28          
## [67] data.table_1.10.4             httpuv_1.3.5                 
## [69] gtable_0.2.0                  reshape_0.8.7                
## [71] mime_0.5                      xtable_1.8-2                 
## [73] AnnotationFilter_1.0.0        survival_2.41-3              
## [75] OrganismDbi_1.18.0            tibble_1.3.4                 
## [77] memoise_1.1.0                 cluster_2.0.6                
## [79] interactiveDisplayBase_1.14.0

Plotting Alignments

Contents

0.1 Working with Reads

0.2 Plotting Alignments

0.3 Provenance