# Chunks options

```{r }
library(knitr)
knitr::opts_chunk$set(fig.width=8, 
    fig.height=6, echo=T, 
    warning=FALSE, message=FALSE,
    prompt=T,tidy=T,tidy.opts=list(width.cutoff=50),
    include=TRUE,cache=TRUE)
# To output a pdf or html document simply substitute github_document with pdf_document or html_document in the header of the file

# Working directory - To be changed once the repository is cloned
dir <- "."
```

# Install packages

```{r message = FALSE, eval = F}
# Bioconductor packages
source("https://bioconductor.org/biocLite.R")
biocLite(c("GenomicRanges","rtracklayer","Rsamtools","GenomicAlignments","TxDb.Hsapiens.UCSC.hg19.knownGene","AnnotationHub"))

# CRAN packages 
install.packages("ggplot2")
install.packages("knitr")
```

# Load packages

```{r}
# Load packages
library("ggplot2")
library("GenomicRanges")
library("Rsamtools")
library("GenomicAlignments")
library("rtracklayer")
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
library("AnnotationHub")
library("knitr")
```


# Advanced `GenomicRanges`

Imagine you have done an experiment on *Sample 1* and your output is a set of genomic ranges. Let's start by manually creating a very simple `GRanges` objects with ranges on *chr1* and *chr2*. 

```{r}
# Create and example: Sample 1
gr_S1 <- GRanges(
seqnames = rep(c("chr1", "chr2"), times = c(3,2)),   
ranges = IRanges(start = c(5,8,20,8,18), 
                 end = c(11,15,26,16,21), 
                 names = c(paste("Range_",1:5,sep = ""))), 
strand = rep(strand(c("*")), times = 5))

gr_S1
```

A very common analysis to perform is to evaluate to what extent and where your ranges overlap with some **features**, such as genes, exons, etc...
As an example we create a simple fake gene annotation to use with the ranges created above.

```{r}
# Gene annotation
genes <- GRanges(
seqnames = rep(c("chr1", "chr2"), times = c(2, 2)),  
ranges = IRanges(start = c(7,17,7,23), 
                 end = c(15,23,14,26), 
                 names = c(paste("Gene_",1:4,sep = ""))), 
strand = rep(strand(c("+")), times = 4))

genes
```

Figure 1 is a simple way of plotting ranges stored into `Granges` object using  `geom_rect()` from `ggplot2`. More advanced plotting options are available but this simple strategy is sufficient for now and also every excuse is good to use `ggplot`!

```{r  fig.cap = "Graphic representation of the GRanges objects created above."}
# Plot ranges with ggplot2
gr <- as.data.frame(rbind(as.data.frame(gr_S1), as.data.frame(genes)))
gr$rangeID <- c(names(gr_S1),names(genes))
gr$Sample <- c(rep("Peaks",length(gr_S1)), rep("Genes",length(genes)))

ggplot(data = gr, aes(xmin = start, xmax = end, ymin = 0, ymax = 1)) + 
  geom_rect( aes(fill = rangeID),alpha = 0.4) + facet_wrap(~ Sample + seqnames) + theme_bw() +
  labs(x = "Genomic position", y = "Ranges") +theme(axis.ticks.y=element_blank(),axis.text.y = element_blank())
```

\clearpage

## Overlaps between two `GRanges` objects

### `findOverlaps()` 

By default looks for overlaps of a minimum of 1bp between a `query` and a `subject`.

```{r}
### `findOverlaps()` 
?findOverlaps
overlaps <- findOverlaps(query = gr_S1, subject = genes)
overlaps
# Query the output
queryHits(overlaps)
subjectHits(overlaps)
subjectLength(overlaps)
```

```{r}
# You can allow for a gap to be ignored 
overlaps <- findOverlaps(gr_S1,genes, maxgap = 5)
overlaps
# You can specify a minimum overlap
overlaps <- findOverlaps(gr_S1,genes, minoverlap = 5)
overlaps
# You can specify a minimum overlap
overlaps <- findOverlaps(gr_S1,genes, type = "start")
overlaps
overlaps <- findOverlaps(gr_S1,genes, type = "end")
overlaps
overlaps <- findOverlaps(gr_S1,genes, type = "within")
overlaps
```

### `countOverlaps()`

```{r}
## `countOverlaps()`
?countOverlaps
```


```{r}
N_overlaps <- countOverlaps(gr_S1,genes)
N_overlaps
# You can play around with the same options as findOverlaps()
```

## Nearest-methods in `GenomicRanges`

### `nearest()`

```{r}
## Nearest-methods in `GenomicRanges`
### `nearest()`
?nearest
```

It returns a vector of indeces referring to the nearest neighbour range in `subject` for every range in `x`.
By default if one range overlaps with multiple genes then one overlap will be chosen at random:

```{r}
nearest(x = gr_S1, subject = genes)
nearest(x = gr_S1,subject = genes, select = "all")
nearest(gr_S1)
nearest(gr_S1,gr_S1)
```

### `distance()`

```{r}
### `distance()`
?distance
```

```{r}
distance(x = gr_S1[1],y = genes)
distance(x = gr_S1[1:4],y = genes)
distance(x = gr_S1, y = genes)
```

`distance()` is a symmetric function which means that it requires x and y to have to have the same lenght and if one is shorter than the other one it will be recycled to match the length of the longest. Also, the distance between two consecutive blocks is 0 not 1 which affects the notion of overlaps. If `distance(x, y) == 0` then x and y can be either adjacent or overlapping ranges.

### `distanceToNearest()`

```{r}
### `distanceToNearest()`
?distanceToNearest
```

For every range in `x` it will return the index and the distance to its nearest neighbour in `subject`.

```{r}
distanceToNearest(x = gr_S1, subject = genes)
```


\clearpage

## `GRangesList`

`GRangesList` are lists of `GRanges` objects. 

```{r}
# Create a GRangesList 
# Sample 1
gr_S1 <- GRanges(
seqnames = rep(c("chr1", "chr2"), times = c(3, 2)),   
ranges = IRanges(start = c(5,8,20,8,18), 
                 end = c(11,15,26,16,21), 
                 names = c(paste("Region_",1:5,sep = ""))), 
strand = rep("*", times = 5) )

gr_S1

# Sample 2
gr_S2 <- GRanges(
seqnames = rep(c("chr2", "chr3"), times = c(3, 5)),   
ranges = IRanges(start = c(1:8), 
                 width = 10, 
                 names = c(paste("Region_",1:8,sep = ""))), 
strand = rep("*", times = 8))

gr_S2

# GRanges List
list_ranges <- GRangesList(Sample1 = gr_S1, Sample2 = gr_S2)
length(list)
```

Many of the functions learnt for `GRanges` can also be applied to `GRangesList` objects even though the output will have to be interepreted accordingly:

```{r}
names(list_ranges)
length(list_ranges)
seqnames(list_ranges)
strand(list_ranges)
ranges(list_ranges)
start(list_ranges)
end(list_ranges)
width(list_ranges)
unlist(list_ranges)
```

To get the number of ranges in every object of the list use `elementNROWS`:

```{r}
elementNROWS(list_ranges) 
```

## Subsetting and looping over `Granges` list

Add an extra column to every `GRanges` in the `GRangesList`.

```{r}
addCols <- lapply(list_ranges, function(x){
elementMetadata(x) <- data.frame(NumberReads = rbinom(length(x),size = 100, prob = 0.5))
return(x)
})
class(addCols)
addCols <- GRangesList(addCols)
class(addCols)
```

In many cases `Granges` objects can be subsetted using the same rules that apply to normal lists or data frames in R:

```{r}
# Subsetting
# As a list
addCols[[1]]
addCols["Sample1"]
# As a data.frame object
addCols[1,"NumberReads"]
addCols["Sample1","NumberReads"]
```

```{r}
# Looping
lapply(addCols, length)
sapply(addCols, length)
```

Similar considerations apply for the other functions explored above like `findOverlaps()`, `countOverlaps()`, etc...

#### Challenge 2

> Use the code below to download the `cpg_islands` object containing the CpG islands for the Human chrY and the `txdb` gene annotation from `TxDb.Hsapiens.UCSC.hg19.knownGene`. 

```{r}
# Download CpG islands from UCSC genome browser
# in the same way you can extract many more genomic object
session <- browserSession("UCSC")
query <- ucscTableQuery(session, "CpG Islands",GRangesForUCSCGenome("hg19", "chrY"))
cpg_islands <- getTable(query)
# Download genes 
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
```

> 1. What is the class of `cpg_islands`? Convert `cpg_islands` to a `GRanges` object.
> 2. Using what you have previosly learnt, get the transcript for every gene from the `txdb` gene annotation and store them into the `transcriptHg19` object (suggestion: see `transcriptsBy`).
> 3. Subset `transcriptHg19` to extract only the transcript on chr21 (suggestion: first run unlist(transcriptHg19)).
> 4. Find the nearest CpG island in `cpg_islandsGR` to every transcript in `transcriptHg19` as well as their distance.
> 5. Optional: On average, every CpG island is close to how many genes?


\clearpage

# `Rtracklayer`

Now that the structure and manipulation of `GRanges` objects is clear, the `rtracklayer` package comes as a useful tool to import/export `GRanges` objects from and to different data formats commonly used in genomic analyses. `Rtracklayer` stands as an interface to mediate the crosstalk between R and genome browsers (UCSC built-in). Example of data files supported by this package are BED, bigWIG, wig, gff and gtf formats.

Another very important function of this package is to allow the visualisations of genomic annotation *tracks*. However, this part will not covered in the present tutorial. For more references see the [Rtracklayer Bioconductor page](http://bioconductor.org/packages/release/bioc/html/rtracklayer.html).

There is a very large number of different data formats used to store different sort of genomic data. In this tutororial only few of them will be covered as example. However, the UCSC genome growser offers a useful [FAQ page](https://genome.ucsc.edu/FAQ/FAQformat.html) where one can find the specification of all mandatory and optional fields for every data format.

## GTF and GFF 

The **GFF** and **GTF** are the preferred format used to store annotations and their exact specification is summarised in [GFF/GTF File Format](http://www.ensembl.org/info/website/upload/gff.html).

The function `import()` is used to import all the supported data types into a `GRanges` object in R. The function recognises the format from the extension of the file but the argument `format` can be used to expliciclty define it. The function `export()` works in the same way and it is used to export `GRanges` objects to files.

```{r }
# Rtracklayer
# Import/Export with Rtrcaklayer
gr_S1
export(gr_S1, file.path(dir,"transcriptHg19.gff"))
export(gr_S1, file.path(dir,"transcriptHg19.gtf"))

importHg19 <- import(file.path(dir,"transcriptHg19.gff"))
class(importHg19)
importHg19
```

#### Challenge 1

> 1. Create a simple `GRanges` object made out of 20 ranges and export it as a `gff` file. 
> 2. Read the file back into R.
> 3. Optional: Print the number of ranges of the `GRanges` object just imported and create a histogram of the widths of the ranges.


## Wiggle (WIG) and bigWIG file formats for graphing tracks

The [WIG](http://www.ensembl.org/info/website/upload/wig.html#tracklines) data format is used for the disply of dense continuous data, like scores. Every line in the file is an intervals and *WIG* files only allow equally spaced intervals. For intervals of variable width one should use the [bedGraph](http://www.ensembl.org/info/website/upload/bed.html#bedGraph) format. 
There are two main data formats: the *variableStep* and the *fixedStep*. The *WIG* file is constituted by one or more blocks separated by declaration lines. Let's go through them with two examples.

- **variableStep**

The following example is a *WIG* file with only one block of coordinates. The field *chrom* is required.

```{bash eval = FALSE, prompt = FALSE}
variableStep chrom=chr2
300701  12.5
300702  12.5
300703  12.5
300704  12.5
300705  12.5
```

The same file can be defined as follows using the *span* argument.

```{bash eval = FALSE, prompt = FALSE}
variableStep chrom=chr2 span=5
300701  12.5
```

- **fixedStep**

This format allows a more compact way of storing intervals. In this situation the *span* would not change the dimension of the file.

```{bash eval = FALSE, prompt = FALSE}
fixedStep chrom=chr3 start=400601 step=100
11
22
33
```

The [bigWig](https://genome.ucsc.edu/goldenpath/help/wiggle.html) file is created from a *WIG* file and in `R` this can be achieved through `wigToBigWig()`. *bigWig* is the recommended format for very large data tracks due to the way in which data are compressed and stored in an indexed binary format. Usually whole-genome coverage vectors are stored as *bigWig*. Howver, the loss is negligible when dealing with very large amount of data. 

```{r }
## Wiggle (WIG) and bigWIG
## Import
wig_path <- file.path(dir, "exampleWIG.wig")
import_wig <- import(con = wig_path,seqinfo = Seqinfo(genome="mm10"))
import_wig
?Seqinfo
seqinfo(import_wig)
```

Convert the *WIG* file into a *bigWig* file and import the *bigWig* file into R. 

```{r}
bigwig_path <- file.path(dir,"example_wigToBigWig.bw")
wigToBigWig(wig_path, dest = bigwig_path, seqinfo = seqinfo(import_wig))
import_bigwig <- import(con = bigwig_path)
```

## Import only a region of the bigWig file

The `import` functions allows to load into a R only a specific regions defined using a `GRanges` object.

```{r }
# Define a region
which_region <- GRanges(seqnames = "chr3", IRanges(start = 1, end = 500000))
import_bigwig_region <- import(con = bigwig_path, which = which_region)
```

# `liftOver` from different genome releases

This function allows to convert genomic coordinate between difference releases of the same genome or even between different species. 

```{r}
?AnnotationHub
ahub <- AnnotationHub()
ahub

# Chain files are data format used to convert coordinates between genome references 
ahub.chain <- subset(ahub, rdataclass == "ChainFile" & species == "Homo sapiens")
ahub.chain2 <- subset(ahub, rdataclass == "ChainFile" & species %in% c("Homo sapiens","Mus musculus"))
query(ahub.chain, c("hg18", "hg19"))

chain <- ahub.chain[ahub.chain$title == "hg19ToHg18.over.chain.gz"]
chain <- chain[[1]]
gr.hg18 <- liftOver(import_bigwig, chain)
gr.hg18
unlist(gr.hg18)
```

\clearpage

# `Rsamtools`

[Rsamtools](http://bioconductor.org/packages/release/bioc/html/Rsamtools.html) is an R interface for [Samtools](http://samtools.sourceforge.net/) which offers a large variety of utilities to manipulate *SAM* and *BAM* files. The main purose of [Rsamtools](http://bioconductor.org/packages/release/bioc/html/Rsamtools.html) is to import *BAM* files into R. This will allow the user to then create objects which can be used for several types of downstream analyses. *BAM* files are the standard way of storing 'short' aligned (and unaligned) reads produced by any type of experiment (RNA-Seq, ChIP-Seq, Methylation data, etc). *BAM* files are binary versions of *SAM* files and can be indexed allowing to access only localised regions of the genome, similar to what we have seen for the *bigWig* file. 
Each read is stored in a *BAM* file with several other measures like base quality, position of the 5' end of the read, read name, etc .... [Rsamtools](http://bioconductor.org/packages/release/bioc/html/Rsamtools.html) allows the user to specify which parameters of the reads to load into R. A detailed description of the fields in the *BAM* and their names is available at http://samtools.github.io/hts-specs/SAMv1.pdf. It has to be remembered that *BAM* files can be very large and there are limits in how much can be done into R and for more complex analyses is preferred to work with [Samtools](http://samtools.sourceforge.net/) is prefered. 

```{bash eval = FALSE}

D00626:239:CAFMCANXX:4:1201:16222:93271 403  chr1    17018   3  38M177N62M      =       16965   -330 TGGCCCAGGTCTGGCACATAGAAGTAGTTCTCTGGGACCTGCTGTTCCAGCTGCTCT  GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
RG:Z:9_CAFMCANXX_L004   NH:i:2  HI:i:2  NM:i:0  nM:i:1  AS:i:200

```

A very important field which is worth understanding is the *FLAG* of the read. This is a binary number which defines the type of read. In the example above it is the number *403*. For example, the *FLAG* defines whether the read is unmapped or mapped or with paired end (PE) reads if the first mate is mapped and the second one is not, and so on. It is not always trivial to decode a *FLAG* number and the [Broad Institute](https://www.broadinstitute.org/) offers a useful online tool to do that for you https://broadinstitute.github.io/picard/explain-flags.html.

## Input BAM into `R`

The function used to load BAM files into R is `scanBam()`. The function `ScanBamParam()` is used to specify the fields and the region/s of the *BAM* file to load into R. The genomic region to load is defined through `GRanges`. Below is an example on how to load a few regions of 1000 bp width from a *BAM* file containing reads on chromosome 21 from an RNA-Seq experiment. We also require only *uniquely mapped* and *properly paired* reads to be loaded as well as we do not load PCR duplicates.
The help page of `?scanBam` lists the name of the fields that can be accessed through the `what` argument in the `ScanBamParam()` function.

```{r}
# Path to the bamfile. In the same folder there has to be the chr1.bam.bai file
bamfile <- file.path(dir,"chr21.bam")
# Define parameters
which <- GRanges(seqnames = rep("chr21",3), 
                 IRanges(start = c(9827000, 16267000, 15890000), 
                         width = 10000)) 
what <- c("rname", "strand", "pos", "qwidth", "seq", "isize")
# ?scanBamFlag
flag = scanBamFlag(isDuplicate = FALSE, isProperPair = TRUE, isPaired = TRUE)

param <- ScanBamParam(which = which, what = what, flag = flag)
# load reads into R
reads = scanBam(bamfile, param = param)

# Explore the reads
names(reads)
region1 <- reads$`chr21:9827000-9836999`
names(region1)
head(region1$rname)
head(region1$pos)
head(region1$qwidth)
head(region1$seq)
head(region1$isize)
```

`readGAlignments` also loads *BAM* files into R directly into a `GRanges` object. The parameters are specified as in `RSamtools` with the argument `param`. 

```{r}
# GenomicAlignments
readsGA <- readGAlignments(bamfile, param = param)
```

#### Challenge 3

> 1. Define a new `ScanBamParam()` object that satisfies the following options:  use the same ranges as above, load both main alignments and PCR duplicates, properly paired and uniquely mapped reads.  
> 2. Plot an histogram of the fragment sizes of every read (Suggestion: look for the `isize` fields)

## Compute number of reads that falls within 100bp bins

It is often useful to summarise the number of reads falling within a bin to have an idea about the distribution of the reads across a region. There are several ways in which this can be accomplished and here two will be shown. Often, a read is counted in a bin if its 5' end falls within bin. This is the approach that we will consider.

```{r}
## Compute number of reads that falls within 100bp bins
bins <- GRanges(seqnames = rep("chr21",11), 
                 IRanges(start = seq(9827000,9827000+10000, by = 1000), 
                         width = 1000)) 
# Create vector of 0-1 positions
positions <- region1$pos
vector_positions <- rep(0,max(positions))
vector_positions[positions] <- 1
# Views Object
rangeViews <- Views(vector_positions, start = seq(min(positions),max(positions), length.out = 10), width = 60)
# viewSums
bin_sum <- viewSums(rangeViews)
# Add extra column to the ranges object
binned_counts <- ranges(rangeViews)
values(binned_counts) <- data.frame(ReadCounts = bin_sum)
```

