```{r init, include=FALSE}
library(knitr)
opts_chunk$set(message=FALSE)
```

```{r style, echo = FALSE, results = 'asis'}
library(BiocStyle)
BiocStyle::markdown()
```

# Getting started

Genomic ranges are a way of describing regions on the genome (or any other
linear object, such as a transcript, or even a protein). This functionality
is typically found in the `r Biocpkg("GenomicRanges")` package.

```{r loadlib}
library(GenomicRanges)
```

## Example GRanges

To get going, we can construct a `GRanges` object by hand as an example.

```{r mkexample}
gr <- GRanges(
    seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
    ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
    strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
    score = 1:10,
    GC = seq(1, 0, length=10))
gr
```

The `GRanges` class represents a collection of genomic ranges that each have a single start
and end location on the genome. It can be used to store the location of genomic features
such as contiguous binding sites, transcripts, and exons. These objects can be created by
using the GRanges constructor function. The following code just creates a `GRanges` object 
from scratch.




This creates a `GRanges` object with 10 genomic ranges. The output of the `GRanges` `show()` method
separates the information into a left and right hand region that are separated by `|` symbols.
The genomic coordinates (seqnames, ranges, and strand) are located on the left-hand side
and the metadata columns are located on the right. For this example, the
metadata is comprised of score and GC information, but almost anything can be stored in
the metadata portion of a GRanges object.

The components of the genomic coordinates within a GRanges object can be extracted using
the seqnames, ranges, and strand accessor functions.

```{r GRanges-location-accessors}
seqnames(gr)
ranges(gr)
strand(gr)
```

Note that the `GRanges` object has information to the "left" side of the `|` that has special 
accessors. The information to the right side of the `|`, when it is present, is the metadata
and is accessed using `mcols()`, for "metadata columns".

```{r mcols-example}
class(mcols(gr))
mcols(gr)
```

Since the `class` of `mcols(gr)` is `r class(mcols(gr))`, we can use our `DataFrame` approaches to 
work with the data.

```{r mcols-df-methods}
mcols(gr)$score
```

We can even assign a new column.

```{r new-mcols-column}
mcols(gr)$AT = 1-mcols(gr)$GC
gr
```


Another common way to create a `GRanges` object is to start with a `data.frame`, perhaps created
by hand like below or read in using `read.csv` or `read.table`. We can convert from a 
`data.frame`, when columns are named appropriately, to a `GRanges` object. 

```{r granges-from-data.frame}
df_regions = data.frame(chromosome = rep("chr1",10),
                        start=seq(1000,10000,1000),
                        end=seq(1100, 10100, 1000))
as(df_regions,'GRanges') # note that names have to match with GRanges slots
# fix column name
colnames(df_regions)[1] = 'seqnames'
gr2 = as(df_regions,'GRanges')
gr2
```

`GRanges` have one-dimensional-like behavior. For instance, we can check the `length` and 
even give `GRanges` names.

```{r names-and-length}
names(gr)
length(gr)
```

## Subsetting GRanges objects

While `GRanges` objects look a bit like a `data.frame`, they can be thought of
as vectors with associated ranges. Subsetting, then, works very similarly to
vectors. To subset a `GRanges` object to include only second and third regions:

```{r subsetgr1}
gr[2:3]
```

That said, if the `GRanges` object has metadata columns, we can also treat
it like a two-dimensional object kind of like a data.frame. Note that
the information to the left of the `|` is not like a `data.frame`, so
we *cannot* do something like `gr$seqnames`. Here is an example of subsetting
with the subset of one metadata column.

```{r subsetgr2}
gr[2:3, "GC"]
```

The usual `head()` and `tail()` also work just fine.

```{r othergrstuff}
head(gr,n=2)
tail(gr,n=2)
```

## Interval operations on one GRanges object

### Intra-range methods

The methods described in this section work *one-region-at-a-time* and are, therefore, called
"intra-region" methods. Methods that work across all regions are described below 
in the [Inter-range methods] section.

The `GRanges` class has accessors for the "ranges" part of the data. For example:

```{r basicgrfuncs}
# Make a smaller GRanges subset
g <- gr[1:3]
start(g) # to get start locations
end(g)   # to get end locations
width(g) # to get the "widths" of each range
range(g) # to get the "range" for each sequence (min(start) through max(end))
```

The GRanges class also has many methods for manipulating the ranges. The methods can
be classified as intra-range methods, inter-range methods, and between-range methods.
Intra-range methods operate on each element of a GRanges object independent of the other
ranges in the object. For example, the flank method can be used to recover regions flanking
the set of ranges represented by the GRanges object. So to get a GRanges object containing
the ranges that include the 10 bases *upstream* of the ranges:

```{r grflank}
flank(g, 10)
```

Note how `flank` pays attention to "strand". 
To get the flanking regions *downstream* of the ranges, we can do:

```{r grflank2}
flank(g, 10, start=FALSE)
```

Other examples of intra-range methods include `resize` and `shift`. The shift method will
move the ranges by a specific number of base pairs, and the resize method will extend the
ranges by a specified width.

```{r shiftandresize}
shift(g, 5)
resize(g, 30)
```

The `r Biocpkg("GenomicRanges")` help page `?"intra-range-methods"` summarizes these methods.

### Inter-range methods

Inter-range methods involve comparisons between ranges in a single GRanges object. For
instance, the reduce method will align the ranges and merge overlapping ranges to produce
a simplified set.

```{r grreduce}
reduce(g)
```

The reduce method could, for example, be used to collapse individual overlapping coding
exons into a single set of coding regions.

Sometimes one is interested in the gaps or the qualities of the gaps between the ranges
represented by your GRanges object. The gaps method provides this information:

```{r}
gaps(g)
```

In this case, we have not specified the lengths of the chromosomes, so Bioconductor is 
making the assumption (incorrectly) that the chromosomes end at the largest location
on each chromosome. We can correct this by setting the `seqlengths` correctly, but
we can ignore that detail for now.

The disjoin method represents a GRanges object as a collection of non-overlapping ranges:

```{r disjoin}
disjoin(g)
```


The `coverage` method quantifies the degree of overlap for all the ranges in a GRanges object.

```{r coveragegr1}
coverage(g)
```

The coverage is summarized as a `list` of coverages, one for each chromosome. 
The `Rle` class is used to store the values. Sometimes, one must convert these 
values to `numeric` using `as.numeric`. In many cases, this will happen automatically,
though. 

```{r coveragegr2}
covg = coverage(g)
covg_chr2 = covg[['chr2']]
plot(covg_chr2, type='l')
```

See the `r Biocpkg("GenomicRanges")` help page `?"intra-range-methods"` for more details.

## Set operations for GRanges objects

Between-range methods calculate relationships between different GRanges objects. Of central
importance are findOverlaps and related operations; these are discussed below. Additional
operations treat GRanges as mathematical sets of coordinates; union(g, g2) is the union of
the coordinates in g and g2. Here are examples for calculating the union, the intersect and
the asymmetric difference (using setdiff).




```{r intervalsgr1}
g2 <- head(gr, n=2)
union(g, g2)
intersect(g, g2)
setdiff(g, g2)
```

There is extensive additional help available or by looking at the vignettes
in at the `r Biocpkg("GenomicRanges")` pages.

```{r helppagesgr,eval=FALSE}
?GRanges
```

There are also many possible `methods` that work with `GRanges` objects. To
see a complete list (long), try:

```{r granges-methods-help, eval=FALSE}
methods(class="GRanges")
```


# GRangesList

Some important genomic features, such as spliced transcripts that are are comprised of exons,
are inherently compound structures. Such a feature makes much more sense when expressed
as a compound object such as a GRangesList. If we thing of each transcript as a set of exons,
each transcript would be summarized as a `GRanges` object. However, if we have multiple transcripts,
we want to somehow keep them separate, with each transcript having its own exons. The `GRangesList`
is then a list of `GRanges` objects that. Continuing with the transcripts thought, a `GRangesList`
can contain all the transcripts and their exons; each transcript is an element in the list. 

Whenever genomic features consist of multiple
ranges that are grouped by a parent feature, they can be represented as a GRangesList
object. Consider the simple example of the two transcript GRangesList below created using
the GRangesList constructor.


```{r example-granges-list-1}
gr1 <- GRanges(
    seqnames = "chr2", 
    ranges = IRanges(103, 106),
    strand = "+", 
    score = 5L, GC = 0.45)
gr2 <- GRanges(
    seqnames = c("chr1", "chr1"),
    ranges = IRanges(c(107, 113), width = 3),
    strand = c("+", "-"), 
    score = 3:4, GC = c(0.3, 0.5))
```

The `gr1` and `gr2` are each `GRanges` objects. We can combine them into a "named"
`GRangesList` like so:

```{r example-granges-list-2}
grl <- GRangesList("txA" = gr1, "txB" = gr2)
grl
```

The `show` method for a `GRangesList` object displays it as a named list of `GRanges` objects,
where the names of this list are considered to be the names of the grouping feature. In the
example above, the groups of individual exon ranges are represented as separate `GRanges`
objects which are further organized into a list structure where each element name is a transcript
name. Many other combinations of grouped and labeled `GRanges` objects are possible
of course, but this example is a common arrangement.

In some cases, `GRangesList`s behave quite similarly to `GRanges` objects.

## Basic _GRangesList_ accessors

Just as with GRanges object, the components of the genomic coordinates within a GRangesList
object can be extracted using simple accessor methods. Not surprisingly, the GRangesList
objects have many of the same accessors as GRanges objects. The difference is that many
of these methods return a list since the input is now essentially a list of GRanges objects.
Here are a few examples:

```{r basicGRLaccessors}
seqnames(grl)
ranges(grl)
strand(grl)
```

The length and names methods will return the length or names of the list and the seqlengths
method will return the set of sequence lengths.

```{r notthesameasgranges}
length(grl)
names(grl)
seqlengths(grl)
```


# Relationships between region sets

One of the more powerful approaches to genomic data integration is to ask
about the relationship between sets of genomic ranges. The key features of 
this process are to look at overlaps and distances to the nearest feature.
These functionalities, combined with the operations like `flank` and `resize`, 
for instance, allow pretty useful analyses with relatively little code.
In general, these operations are *very* fast, even on thousands to millions
of regions. 

## Overlaps



```{r findoverlaps1}
mtch <- findOverlaps(gr, g)
as.matrix(mtch)
```

If you are interested in only the `queryHits` or the `subjectHits`, there are 
accessors for those (ie., `queryHits(match)`).


As you might expect, the `countOverlaps` method counts the regions in the
second `GRanges` that overlap with those that overlap
with each element of the first.

```{r countoverlaps}
countOverlaps(gr, g)
```

The `subsetByOverlaps` method simply subsets the first `GRanges` object to include
*only* those that overlap the second. 

```{r subsetByOverlaps}
subsetByOverlaps(gr, g)
```

In some cases, you may be interested in only one hit when doing 
overlaps. Note the `select` parameter. See the help for `findOverlaps`

```{r select-first}
findOverlaps(gr, g, select="first")
findOverlaps(g, gr, select="first")
```

The `%over%` logical operator allows us to do similar things to 
`findOverlaps` and `subsetByOverlaps`. 

```{r}
gr %over% g
gr[gr %over% g]
```


## Nearest feature

There are a number of useful methods that find the nearest feature (region) in a second set
for each element in the first set. 

We can review our two `GRanges` toy objects:

```{r gr-review}
g
gr
```

- nearest: Performs conventional nearest neighbor finding. Returns an integer vector containing the index of the nearest neighbor range in subject for each range in x. If there is no nearest neighbor NA is returned. For details of the algorithm see the man page in the IRanges package (?nearest).

- precede: For each range in x, precede returns the index of the range in subject that is directly preceded by the range in x. Overlapping ranges are excluded. NA is returned when there are no qualifying ranges in subject.

- follow: The opposite of precede, follow returns the index of the range in subject that is directly followed by the range in x. Overlapping ranges are excluded. NA is returned when there are no qualifying ranges in subject.

Orientation and strand for precede and follow: Orientation is 5' to 3', consistent with the direction of translation. Because positional numbering along a chromosome is from left to right and transcription takes place from 5' to 3', precede and follow can appear to have ‘opposite’ behavior on the + and - strand. Using positions 5 and 6 as an example, 5 precedes 6 on the + strand but follows 6 on the - strand.

The table below outlines the orientation when ranges on different strands are compared. In general, a feature on * is considered to belong to both strands. The single exception is when both x and subject are * in which case both are treated as +.


```
       x  |  subject  |  orientation 
     -----+-----------+----------------
a)     +  |  +        |  ---> 
b)     +  |  -        |  NA
c)     +  |  *        |  --->
d)     -  |  +        |  NA
e)     -  |  -        |  <---
f)     -  |  *        |  <---
g)     *  |  +        |  --->
h)     *  |  -        |  <---
i)     *  |  *        |  --->  (the only situation where * arbitrarily means +)
```


```{r nearest}
res = nearest(g, gr)
res
```

While `nearest` and friends give the index of the nearest feature, the distance
to the nearest is sometimes also useful to have. The `distanceToNearest`
method calculates the nearest feature as well as reporting the distance.

```{r distanceToNearest}
res = distanceToNearest(g, gr)
res
```


# Session Info

```{r sessionInfo}
sessionInfo()
```



