```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Exploring VCFs

We read the `.vcf` file generated by using `freebayes` to call germline SNVs in chromosome 20 of the Hapmap sample *NA12878*.

- `rowRanges` extracts `GRanges` information from the object created by `readVcf`

```{r}
library(VariantAnnotation)
NA12878.calls <- readVcf("NA12878.chr20.freebayes.vcf","hg19")
NA12878.calls.ranges <- rowRanges(NA12878.calls)
NA12878.calls.ranges
```


The `Org.hs.eg.db` package can be used to convert between different identifer types. We want the *Entrez ID* for the gene with symbol *PRND*. Types of identifier that can be used to make queries are given by `keytypes`

```{r}
library(org.Hs.eg.db)

keytypes(org.Hs.eg.db)
columns(org.Hs.eg.db)

select(org.Hs.eg.db, keys="PRND",
       keytype = "SYMBOL",
       columns = "ENTREZID")
```

`Txdb` packages in Bioconductor provide gene models for a given organism and genome build. 

```{r}
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
```

We can quickly generate an object that contains the gene models for *all* genes. The object is a *list*, with each entry being indexed by Entrez Gene ID. 

```{r}
all.exons <- exonsBy(txdb,"gene")
all.exons
```

In R, lists are indexed by using the `[[]]` notation. We have already obtained the Entrez Gene ID for our gene of interest (*PRND*) from a previous step.

```{r}
my.gene <- all.exons[["23627"]]
my.gene
```

The chromosome naming conventions of our variant calls (`NA12878.calls.ranges`) and gene are not compatible. A convenient shortcut to re-name the chromosome is `seqlevelsStyle`.

```{r}
seqlevelsStyle(NA12878.calls.ranges)
seqlevelsStyle(my.gene)

seqlevelsStyle(NA12878.calls.ranges) <- "UCSC"
NA12878.calls.ranges <- keepSeqlevels(NA12878.calls.ranges, "chr20")
```


The `%over%` function will perform an overlap between two sets of ranges, a returns `TRUE` or `FALSE` for each interval in the *first* ranges object if it occurs in the *second*. The logical vector can be used to used the `GRanges` object and return the variants that overlap with the gene.

```{r}
in.gene <- NA12878.calls.ranges %over% my.gene
NA12878.calls.ranges[in.gene]
```

If we just want counts, the `countOverlaps` function can be used. Here, the gene has two exons but all variants occur in the second exon

```{r}
countOverlaps(my.gene, NA12878.calls.ranges)
```

## Exercise

The file `regions.of.interest.bed` contains further regions that we want to count the number of variants for. It can be imported using the `rtracklayer` package



```{r}
library(rtracklayer)
regions <- import("regions.of.interest.bed")
varCounts <- countOverlaps(regions, NA12878.calls.ranges)
varCounts
```

How can we produce a `.bed` file with the coordinates of all exons on chromosome 20? Firstly, remember `all.exons` has all exons in the genome, but is a `list` object.

```{r}
all.exons <- unlist(all.exons)
chr20.exons <- all.exons[seqnames(all.exons) == "chr20"]
chr20.exons
export(chr20.exons, con="chr20.exons.bed")
```

Repeating the counting as above, except for introns on chromosome 20. Firstly, `txdb` packages have similar functions to `exonsBy` to obtain other types of genomic feature.

```{r}
all.introns <- intronsByTranscript(txdb)
all.introns <- unlist(all.introns)
intronCounts <- countOverlaps(all.introns, NA12878.calls.ranges)

```

