About This Document »

Package name liftOver
Built with Bioconductor (R) 3.6 (3.4.1)
Last Built Fri, 06 Oct 2017 11:12:44 -0700
Last Modified Fri, 06 Oct 2017 10:31:43 -0700 (r132129)
Source Package liftOver_1.0.2.tar.gz
Windows Binary liftOver_1.0.2.zip
R Script liftOver.R

To install this workflow under Bioconductor 3.6, start R and enter:

source("http://bioconductor.org/workflows.R")
workflowInstall("liftOver")

Changing genomic coordinate systems with rtracklayer::liftOver

Contents

The liftOver facilities developed in conjunction with the UCSC browser track infrastructure are available for transforming data in GRanges formats. This is illustrated here with an image of the EBI/NHGRI GWAS catalog that is, as of May 10 2017, distributed with coordinates defined by NCBI build hg38.

Setup: The NHGRI GWAS catalog as an hg38-based GRanges

## Loading required package: Homo.sapiens

## Loading required package: AnnotationDbi

## Loading required package: stats4

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, cbind, colMeans, colnames,
##     colSums, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, lengths, Map, mapply, match,
##     mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rowMeans, rownames, rowSums, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which, which.max, which.min

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## Loading required package: IRanges

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: OrganismDbi

## Loading required package: GenomicFeatures

## Loading required package: GenomeInfoDb

## Loading required package: GenomicRanges

## No methods found in "RSQLite" for requests: dbGetQuery

## Loading required package: GO.db

## 

## Loading required package: org.Hs.eg.db

## 

## Loading required package: TxDb.Hsapiens.UCSC.hg19.knownGene

## gwascat loaded.  Use data(ebicat38) for hg38 coordinates;

##  data(ebicat37) for hg19 coordinates.

library(gwascat)
cur = makeCurrentGwascat()  # result varies by day

data(cur)
cur

## gwasloc instance with 36761 records and 37 attributes per record.
## Extracted:  2017-05-10 
## Genome:  GRCh38 
## Excerpt:
## GRanges object with 5 ranges and 3 metadata columns:
##       seqnames                 ranges strand | DISEASE/TRAIT        SNPS
##          <Rle>              <IRanges>  <Rle> |   <character> <character>
##   [1]        1 [203186754, 203186754]      * | YKL-40 levels   rs4950928
##   [2]       13 [ 39776775,  39776775]      * |     Psoriasis   rs7993214
##   [3]       15 [ 78513681,  78513681]      * |   Lung cancer   rs8034191
##   [4]        1 [159711078, 159711078]      * |   Lung cancer   rs2808630
##   [5]        3 [190632672, 190632672]      * |   Lung cancer   rs7626795
##         P-VALUE
##       <numeric>
##   [1]     1e-13
##   [2]     2e-06
##   [3]     3e-18
##   [4]     7e-06
##   [5]     8e-06
##   -------
##   seqinfo: 23 sequences from GRCh38 genome

Resource: The chain file for hg38 to hg19 transformation

The transformation to hg19 coordinates is defined by a chain file provided by UCSC. rtracklayer::import.chain will bring the data into R.

library(rtracklayer)
path = system.file(package="liftOver", "extdata", "hg38ToHg19.over.chain")
ch = import.chain(path)
ch

## Chain of length 25
## names(25): chr22 chr21 chr19 chr20 chrY chr18 ... chr6 chr5 chr4 chr3 chr2 chr1

str(ch[[1]])

## Formal class 'ChainBlock' [package "rtracklayer"] with 6 slots
##   ..@ ranges  :Formal class 'IRanges' [package "IRanges"] with 6 slots
##   .. .. ..@ start          : int [1:6842] 16367189 16386933 16386970 16387001 16387128 16395491 16395528 16395841 16395860 16395956 ...
##   .. .. ..@ width          : int [1:6842] 19744 36 31 112 8362 36 312 18 95 33 ...
##   .. .. ..@ NAMES          : NULL
##   .. .. ..@ elementType    : chr "integer"
##   .. .. ..@ elementMetadata: NULL
##   .. .. ..@ metadata       : list()
##   ..@ offset  : int [1:6842] -480662 -480702 -480702 -480726 -480726 -480726 -480726 -480726 -480726 -480726 ...
##   ..@ score   : int [1:1168] -1063867308 68830488 21156147 20814926 7358950 3927744 2928210 991419 880681 802146 ...
##   ..@ space   : chr [1:1168] "chr22" "chr14" "chr22" "chr21" ...
##   ..@ reversed: logi [1:1168] FALSE FALSE FALSE FALSE FALSE FALSE ...
##   ..@ length  : int [1:1168] 1124 1280 173 465 398 110 43 173 342 84 ...

Some more details about the chain data structure are available in the import.chain man page

   A chain file essentially details many local alignments, so it is
   possible for the "from" ranges to map to overlapping regions in
   the other sequence. The "from" ranges are guaranteed to be
   disjoint (but do not necessarily cover the entire "from"
   sequence).

Action: liftOver

The liftOver function will create a GRangesList.

seqlevelsStyle(cur) = "UCSC"  # necessary
cur19 = liftOver(cur, ch)
class(cur19)

## [1] "GRangesList"
## attr(,"package")
## [1] "GenomicRanges"

We unlist and coerce to the gwaswloc class, a convenient form for the GWAS catalog with its many mcols fields.

cur19 = unlist(cur19)
genome(cur19) = "hg19"
cur19 = new("gwaswloc", cur19)
cur19

## gwasloc instance with 36740 records and 37 attributes per record.
## Extracted:   
## Genome:  hg19 
## Excerpt:
## GRanges object with 5 ranges and 3 metadata columns:
##       seqnames                 ranges strand | DISEASE/TRAIT        SNPS
##          <Rle>              <IRanges>  <Rle> |   <character> <character>
##   [1]     chr1 [203155882, 203155882]      * | YKL-40 levels   rs4950928
##   [2]    chr13 [ 40350912,  40350912]      * |     Psoriasis   rs7993214
##   [3]    chr15 [ 78806023,  78806023]      * |   Lung cancer   rs8034191
##   [4]     chr1 [159680868, 159680868]      * |   Lung cancer   rs2808630
##   [5]     chr3 [190350461, 190350461]      * |   Lung cancer   rs7626795
##         P-VALUE
##       <numeric>
##   [1]     1e-13
##   [2]     2e-06
##   [3]     3e-18
##   [4]     7e-06
##   [5]     8e-06
##   -------
##   seqinfo: 23 sequences from hg19 genome; no seqlengths

We see that the translation leads to a loss of some loci.

length(cur)-length(cur19)

## [1] 21

setdiff(cur$SNPs, cur19$SNPs)

## NULL

It may be interesting to follow up some of the losses.