---
title: "1. Introduction to VariantAnnotation"
author: "Valerie Obenchain"
date: "`r format(Sys.time(), '%B %d, %Y')`"
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{1. Introduction to VariantAnnotation}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document:
    number_sections: yes
    toc: yes
    toc_depth: 4
---

# Introduction

This vignette outlines a work flow for annotating and filtering genetic
variants using the `r Biocpkg("VariantAnnotation")` package. Sample data are
in VariantCall Format (VCF) and are a subset of chromosome 22 from [1000
Genomes](http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20110521/).
VCF text files contain meta-information lines, a header line with column
names, data lines with information about a position in the genome, and
optional genotype information on samples for each position. Samtools
organisation and repositories describes the [VCF
format](http://samtools.github.io/hts-specs/VCFv4.4.pdf) in detail.

Data are read in from a VCF file and variants identified according to region
such as `coding`, `intron`, `intergenic`, `spliceSite` etc. Amino acid coding
changes are computed for the non-synonymous variants and SIFT and PolyPhen
databases provide predictions of how severly the coding changes affect protein
function.

# Variant Call Format (VCF) files

## Data import and exploration

Data are parsed into a `VCF` object with `readVcf`.

```{r readVcF,message=FALSE}
library(VariantAnnotation)
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
vcf
```

### Header information

Header information can be extracted from the VCF with `header()`. We see there
are 5 samples, 1 piece of meta information, 22 info fields and 3 geno fields.

```{r readVcf_showheader}
header(vcf)
```

Data can be further extracted using the named accessors.

```{r headeraccessors}
samples(header(vcf))
geno(header(vcf))
```

### Genomic positions

`rowRanges` contains information from the CHROM, POS, and ID fields of the VCF
file, represented as a `GRanges`. The `paramRangeID` column is meaningful when
reading subsets of data and is discussed further below.

```{r readVcf_rowRanges}
head(rowRanges(vcf), 3)
```

Individual fields can be pulled out with named accessors. Here we see `REF` is
stored as a `DNAStringSet` and `qual` is a numeric vector. 

```{r readVcf_fixed}
ref(vcf)[1:5]
qual(vcf)[1:5]
```

`ALT` is a `DNAStringSetList` (allows for multiple alternate alleles per
variant) or a `DNAStringSet`. When structural variants are present it will be a
`CharacterList`.

```{r readVcf_ALT}
alt(vcf)[1:5]
```

### Genotype data

Genotype data described in the `FORMAT` fields are parsed into the geno slot.
The data are unique to each sample and each sample may have multiple values
variable. Because of this, the data are parsed into matrices or arrays where the
rows represent the variants and the columns the samples.  Multidimentional
arrays indicate multiple values per sample. In this file all variables are
matrices.

```{r geno_hdr}
geno(vcf)
sapply(geno(vcf), class)
```

Let's take a closer look at the genotype dosage (DS) variable. The header
provides the variable definition and type.

```{r explore_geno}
geno(header(vcf))["DS",]
```

These data are stored as a 10376 x 5 matrix. Each of the five samples (columns)
has a single value per variant location (row).

```{r dim_geno}
DS <-geno(vcf)$DS
dim(DS)
DS[1:3,]
```

DS is also known as 'posterior mean genotypes' and range in value from [0, 2].
To get a sense of variable distribution, we compute a five number summary of the
minimum, lower-hinge (first quartile), median, upper-hinge (third quartile) and
maximum.

```{r fivenum}
fivenum(DS)
```

The majority of these values (86%) are zero.

```{r DS_zero}
length(which(DS==0))/length(DS)
```

View the distribution of the non-zero values.

```{r DS_hist, fig=TRUE}
hist(DS[DS != 0], breaks=seq(0, 2, by=0.05),
    main="DS non-zero values", xlab="DS")
```

### Info data 

In contrast to the genotype data, the info data are unique to the variant and
the same across samples. All info variables are represented in a single
`DataFrame`.

```{r info}
info(vcf)[1:4, 1:5]
```

We will use the info data to compare quality measures between novel (i.e., not
in dbSNP) and known (i.e., in dbSNP) variants and the variant type present in
the file. Variants with membership in dbSNP can be identified by using the
appropriate SNPlocs package for the hg19 genome (GRCh37).

```{r examine_dbSNP, message=FALSE, warning=FALSE}
library(SNPlocs.Hsapiens.dbSNP144.GRCh37)
vcf_rsids <- names(rowRanges(vcf))
chr22snps <- snpsBySeqname(SNPlocs.Hsapiens.dbSNP144.GRCh37, "22")
chr22_rsids <- mcols(chr22snps)$RefSNP_id
in_dbSNP <- vcf_rsids %in% chr22_rsids
table(in_dbSNP) 
```

Info variables of interest are 'VT', 'LDAF' and 'RSQ'. The header offers
more details on these variables.

```{r header_info}
info(header(vcf))[c("VT", "LDAF", "RSQ"),]
```

Create a data frame of quality measures of interest ...

```{r examine_quality}
metrics <- data.frame(QUAL=qual(vcf), in_dbSNP=in_dbSNP,
    VT=info(vcf)$VT, LDAF=info(vcf)$LDAF, RSQ=info(vcf)$RSQ)
```

and visualize the distribution of qualities using `r CRANpkg("ggplot2")`.  For
instance, genotype imputation quality is higher for the known variants in dbSNP.

```{r examine_ggplot2, message=FALSE, warning=FALSE, fig=TRUE}
library(ggplot2)
ggplot(metrics, aes(x=RSQ, fill=in_dbSNP)) +
    geom_density(alpha=0.5) +
    scale_x_continuous(name="MaCH / Thunder Imputation Quality") +
    scale_y_continuous(name="Density") +
    theme(legend.position="top")
```

## Import data subsets

When working with large VCF files it may be more efficient to read in subsets of
the data. This can be accomplished by selecting genomic coordinates (ranges) or
by specific fields from the VCF file.

### Select genomic coordinates

To read in a portion of chromosome 22, create a `GRanges` with the regions of
interest.

```{r subset_ranges}
rng <- GRanges(seqnames="22", ranges=IRanges(
           start=c(50301422, 50989541), 
           end=c(50312106, 51001328),
           names=c("gene_79087", "gene_644186")))
```

When ranges are specified, the VCF file must have an accompanying Tabix index
file. See `indexTabix` for help creating an index.

```{r subset_TabixFile}
tab <- TabixFile(fl)
vcf_rng <- readVcf(tab, "hg19", param=rng)
```

The `paramRangesID` column distinguishes which records came from which param
range.

```{r}
head(rowRanges(vcf_rng), 3)
```

### Select VCF fields 

Data import can also be defined by the `fixed`, `info` and `geno` fields. Fields
available for import are described in the header information. To view the header
before reading in the data, use `ScanVcfHeader`.

```{r subset_scanVcfHeader}
hdr <- scanVcfHeader(fl)
## e.g., INFO and GENO fields
head(info(hdr), 3)
head(geno(hdr), 3)
```

To subset on "LDAF" and "GT" we specify them as `character` vectors in the
`info` and `geno` arguments to `ScanVcfParam`. This creates a `ScanVcfParam`
object which is used as the `param` argument to `readVcf`.

```{r subset_ScanVcfParam}
## Return all 'fixed' fields, "LAF" from 'info' and "GT" from 'geno'
svp <- ScanVcfParam(info="LDAF", geno="GT")
vcf1 <- readVcf(fl, "hg19", svp)
names(geno(vcf1))
```

To subset on both genomic coordinates and fields the `ScanVcfParam` object must
contain both.

```{r subset_ScanVcfParam_new}
svp_all <- ScanVcfParam(info="LDAF", geno="GT", which=rng)
svp_all
```

# Locating variants in and around genes

Variant location with respect to genes can be identified with the
`locateVariants` function. Regions are specified in the `region` argument and
can be one of the following constructors: CodingVariants, IntronVariants,
FiveUTRVariants, ThreeUTRVariants, IntergenicVariants, SpliceSiteVariants or
PromoterVariants. Location definitions are shown in Table \@ref(tab:table).

| Location   | Details                                                    |
|------------|------------------------------------------------------------|
| coding     | falls *within* a coding region                             |
| fiveUTR    | falls *within* a 5' untranslated region                    |
| threeUTR   | falls *within* a 3' untranslated region                    |
| intron     | falls *within* an intron region                            |
| intergenic | does not fall *within* a transcript associated with a gene |
| spliceSite | overlaps any portion of the first 2 or last 2              |
| promoter   | falls *within* a promoter region of a transcript           |

: (\#tab:table) Variant locations

For overlap methods to work properly the chromosome names (seqlevels) must be
compatible in the objects being compared. The VCF data chromosome names are
represented by number, i.e., '22', but the TxDb chromosome names are preceded
with 'chr'. Seqlevels in the VCF can be modified with the `seqlevels` function.

```{r locate_rename_seqlevels, message=FALSE, warning=FALSE}
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
seqlevels(vcf) <- "chr22"
rd <- rowRanges(vcf)
loc <- locateVariants(rd, txdb, CodingVariants())
head(loc, 3)
```

Locate variants in all regions with the `AllVariants()` constructor,

```{r AllVariants, eval=FALSE}
allvar <- locateVariants(rd, txdb, AllVariants())
```

To answer gene-centric questions data can be summarized by gene reguardless
of transcript.

```{r locate_gene_centric}
## Did any coding variants match more than one gene?
splt <- split(mcols(loc)$GENEID, mcols(loc)$QUERYID) 
table(sapply(splt, function(x) length(unique(x)) > 1))

## Summarize the number of coding variants by gene ID.
splt <- split(mcols(loc)$QUERYID, mcols(loc)$GENEID)
head(sapply(splt, function(x) length(unique(x))), 3)
```

# Amino acid coding changes

`predictCoding` computes amino acid coding changes for non-synonymous variants.
Only ranges in query that overlap with a coding region in the `subject` are
considered. Reference sequences are retrieved from either a `BSgenome` or fasta
file specified in `seqSource`.  Variant sequences are constructed by
substituting, inserting or deleting values in the `varAllele` column into the
reference sequence. Amino acid codes are computed for the variant codon sequence
when the length is a multiple of 3.

The query argument to `predictCoding` can be a `GRanges` or `VCF`. When a
`GRanges` is supplied the `varAllele` argument must be specified. In the case of
a `VCF`, the alternate alleles are taken from `alt(<VCF>)` and the `varAllele`
argument is not specified.

The result is a modified `query` containing only variants that fall within
coding regions. Each row represents a variant-transcript match so more than one
row per original variant is possible.

```{r predictCoding, warning=FALSE}
library(BSgenome.Hsapiens.UCSC.hg19)
coding <- predictCoding(vcf, txdb, seqSource=Hsapiens)
coding[5:7]
```

Using variant rs114264124 as an example, we see varAllele `A` has been
substituted into the `refCodon` `CGG` to produce `varCodon` `CAG.` The
`refCodon` is the sequence of codons necessary to make the variant allele
substitution and therefore often includes more nucleotides than indicated in the
range (i.e. the range is 50302962, 50302962, width of 1). Notice it is the
second position in the `refCodon` that has been substituted. This position in
the codon, the position of substitution, corresponds to genomic position
50302962. This genomic position maps to position 698 in coding region-based
coordinates and to triplet 233 in the protein. This is a non-synonymous coding
variant where the amino acid has changed from `R` (Arg) to `Q` (Gln).

When the resulting `varCodon` is not a multiple of 3 it cannot be translated.
The consequence is considered a `frameshift` and `varAA` will be missing.

```{r predictCoding_frameshift}
## CONSEQUENCE is 'frameshift' where translation is not possible
coding[mcols(coding)$CONSEQUENCE == "frameshift"]
```

# SIFT and PolyPhen Databases

From `predictCoding` we identified the amino acid coding changes for the
non-synonymous variants. For this subset we can retrieve predictions of how
damaging these coding changes may be. SIFT (Sorting Intolerant From Tolerant)
and PolyPhen (Polymorphism Phenotyping) are methods that predict the impact of
amino acid substitution on a human protein. The SIFT method uses sequence
homology and the physical properties of amino acids to make predictions about
protein function. PolyPhen uses sequence-based features and structural
information characterizing the substitution to make predictions about the
structure and function of the protein.

Collated predictions for specific dbSNP builds are available as downloads from
the SIFT and PolyPhen web sites. These results have been packaged into
`r Biocpkg("SIFT.Hsapiens.dbSNP132")` and
`r Biocpkg("PolyPhen.Hsapiens.dbSNP131")` and are designed to be searched by
rsid. Variants that are in dbSNP can be searched with these database packages.
When working with novel variants, SIFT and PolyPhen must be called directly.
See references for home pages.

Identify the non-synonymous variants and obtain the rsids.

```{r nonsynonymous}
nms <- names(coding) 
idx <- mcols(coding)$CONSEQUENCE == "nonsynonymous"
nonsyn <- coding[idx]
names(nonsyn) <- nms[idx]
rsids <- unique(names(nonsyn)[grep("rs", names(nonsyn), fixed=TRUE)])
```

Detailed descriptions of the database columns can be found with `?SIFTDbColumns`
and `?PolyPhenDbColumns`. Variants in these databases often contain more than
one row per variant. The variant may have been reported by multiple sources and
therefore the source will differ as well as some of the other variables. 

It is important to keep in mind the pre-computed predictions in the SIFT and
PolyPhen packages are based on specific gene models. SIFT is based on Ensembl
and PolyPhen on UCSC Known Gene. The `TxDb` we used to identify the coding snps
was based on UCSC Known Gene so we will use PolyPhen for predictions. PolyPhen
provides predictions using two different training datasets and has considerable
information about 3D protein structure. See `?PolyPhenDbColumns` or the PolyPhen
web site listed in the references for more details.

Query the PolyPhen database,

```{r polyphen, message=FALSE, warning=FALSE}
library(PolyPhen.Hsapiens.dbSNP131)
pp <- select(PolyPhen.Hsapiens.dbSNP131, keys=rsids,
          cols=c("TRAININGSET", "PREDICTION", "PPH2PROB"))
head(pp[!is.na(pp$PREDICTION), ]) 
```

# Other operations

## Create a SnpMatrix

The 'GT' element in the `FORMAT` field of the VCF represents the genotype. These
data can be converted into a `SnpMatrix` object which can then be used with the
functions offered in `r Biocpkg("snpStats")` and other packages making use of the `SnpMatrix` class.

The `genotypeToSnpMatrix` function converts the genotype calls in `geno` to a
`SnpMatrix`. No dbSNP package is used in this computation. The return value is a
named list where 'genotypes' is a `SnpMatrix` and 'map' is a `DataFrame` with
SNP names and alleles at each loci.  The `ignore` column in 'map' indicates
which variants were set to NA (missing) because they met one or more of the
following criteria,

- variants with >1 ALT allele are set to NA
- only single nucleotide variants are included; others are set to NA
- only diploid calls are included; others are set to NA

See ?`genotypeToSnpMatrix` for more details.

```{r snpMatrix, message=FALSE}
res <- genotypeToSnpMatrix(vcf) 
res
```

In the map DataFrame, allele.1 represents the reference allele and allele.2 is
the alternate allele.

```{r snpMatrix_ALT}
allele2 <- res$map[["allele.2"]]
## number of alternate alleles per variant
unique(elementNROWS(allele2))
```

In addition to the called genotypes, genotype likelihoods or probabilities can
also be converted to a `SnpMatrix`, using the `r Biocpkg("snpStats")` encoding
of posterior probabilities as byte values.  To use the values in the 'GL' or
'GP' `FORMAT` field instead of the called genotypes, use the `uncertain=TRUE`
option in `genotypeToSnpMatrix`.

```{r message=FALSE}
fl.gl <- system.file("extdata", "gl_chr1.vcf", package="VariantAnnotation")
vcf.gl <- readVcf(fl.gl, "hg19")
geno(vcf.gl)

## Convert the "GL" FORMAT field to a SnpMatrix
res <- genotypeToSnpMatrix(vcf.gl, uncertain=TRUE)
res
t(as(res$genotype, "character"))[c(1,3,7), 1:5]

## Compare to a SnpMatrix created from the "GT" field
res.gt <- genotypeToSnpMatrix(vcf.gl, uncertain=FALSE)
t(as(res.gt$genotype, "character"))[c(1,3,7), 1:5]

## What are the original likelihoods for rs58108140?
geno(vcf.gl)$GL["rs58108140", 1:5]
```

For variant rs58108140 in sample NA06989, the \"A/B\" genotype is much more
likely than the others, so the `SnpMatrix` object displays the called genotype.

## Write out VCF files

A VCF file can be written out from data stored in a `VCF` class.

```{r writeVcf, message=FALSE, warning=FALSE}
fl <- system.file("extdata", "ex2.vcf", package="VariantAnnotation")
out1.vcf <- tempfile()
out2.vcf <- tempfile()
in1 <- readVcf(fl, "hg19")
writeVcf(in1, out1.vcf)
in2 <- readVcf(out1.vcf, "hg19")
writeVcf(in2, out2.vcf)
in3 <- readVcf(out2.vcf, "hg19")
identical(rowRanges(in1), rowRanges(in3))
identical(geno(in1), geno(in2))
```

# Performance

Targeted queries can greatly improve the speed of data input. When all data from
the file are needed define a `yieldSize` in the `TabixFile` to iterate through
the file in chunks.

```{r eval=FALSE}
readVcf(TabixFile(fl, yieldSize=10000))
```

`readVcf` can be used with a to select any combination of INFO and GENO fields,
samples or genomic positions.

```{r eval=FALSE}
readVcf(TabixFile(fl), param=ScanVcfParam(info='DP', geno='GT'))
```

While `readvcf` offers the flexibility to define combinations of INFO, GENO and
samples in the `ScanVcfParam`, sometimes only a single field is needed. In this
case the lightweight `read` functions (`readGT`, `readInfo` and `readGeno`) can
be used. These functions return the single field as a matrix instead of a `VCF`
object.

```{r eval=FALSE}
readGT(fl)
```

The table below highlights the speed differences of targeted queries vs reading
in all data. The test file is from 1000 Genomes and has 494328 variants, 1092
samples, 22 INFO, and 3 GENO fields and is located at
http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20101123/. `yieldSize` is
used to define chunks of 100, 1000, 10000 and 100000 variants. For each chunk
size three function calls are compared: `readGT` reading only GT, `readVcf`
reading both `GT` and `ALT` and finally `readVcf` reading in all the data.

```{r eval=FALSE}
library(microbenchmark)
fl <- "ALL.chr22.phase1_release_v3.20101123.snps_indels_svs.genotypes.vcf.gz"
ys <- c(100, 1000, 10000, 100000)

## readGT() input only 'GT':
fun <- function(fl, yieldSize) readGT(TabixFile(fl, yieldSize))
lapply(ys, function(i) microbenchmark(fun(fl, i), times=5)

## readVcf() input only 'GT' and 'ALT':
fun <- function(fl, yieldSize, param) 
           readVcf(TabixFile(fl, yieldSize), "hg19", param=param)
param <- ScanVcfParam(info=NA, geno="GT", fixed="ALT")
lapply(ys, function(i) microbenchmark(fun(fl, i, param), times=5)

## readVcf() input all variables:
fun <- function(fl, yieldSize) readVcf(TabixFile(fl, yieldSize), "hg19")
lapply(ys, function(i) microbenchmark(fun(fl, i), times=5))
```

n records   readGT   readVcf (GT and ALT)   readVcf (all)
----------- -------- ---------------------- ---------------
100         0.082    0.128                  0.501
1000        0.609    0.508                  5.878
10000       5.972    6.164                  68.378
100000      78.593   81.156                 693.654

: (\#tab:performance) Targeted queries (time in seconds)

# References

Wang K, Li M, Hakonarson H, (2010), ANNOVAR: functional annotation of
genetic variants from high-throughput sequencing data. Nucleic Acids
Research, Vol 38, No. 16, e164.

McLaren W, Pritchard B, RiosD, et. al., (2010), Deriving the
consequences of genomic variants with the Ensembl API and SNP Effect
Predictor. Bioinformatics, Vol. 26, No. 16, 2069-2070.

SIFT home page: http://sift.bii.a-star.edu.sg/

PolyPhen home page: http://genetics.bwh.harvard.edu/pph2/

# Session Information

```{r sessionInfo, echo=FALSE}
sessionInfo()
```