reveal.js

## Genomic Intervals
---
## Terminology

- **interval**: a location on a sequence defined by a *start* and *end*
- **region**: often used as a synonym of *interval*
- **chromosome**: In the context of genomic intervals, *chromosome* is often used to indicate the *name of the sequence* on which an interval is defined, but this is really a fudging
- **interval set** or **region set**: a set of intervals/regions

---
<img src="images/genomic-intervals/intervals.svg">

---
## What can be represented as an interval?

- ChIP-seq or ATAC-seq peaks

---
## Peaks

![](images/genomic-intervals/peaks.svg)

Genomic intervals are often colloquially referred to as 'peaks'.

---
## What can be represented as an interval?

- ChIP-seq or ATAC-seq peaks
- Single-Nucleotide Polymorphisms (SNPs)

---
## SNPs

![](images/genomic-intervals/snps.svg)

SNPs are interval of width 1

---
## What can be represented as an interval?

- ChIP-seq or ATAC-seq peaks
- Single-Nucleotide Polymorphisms (SNPs)
- Genes and gene components (TSS, exons, introns, etc)

---
## Genes and gene components

![](images/genomic-intervals/brca2-gene-model.png)

---
## What can be represented as an interval?

- ChIP-seq or ATAC-seq peaks
- Single-Nucleotide Polymorphisms (SNPs)
- Genes and gene components (TSS, exons, introns, etc)
- Non-coding DNA annotation (promoters, enhancers)

---
## Non-coding DNA annotation

![](images/genomic-intervals/regulatory-build.png)

---
## What can be represented as an interval?

- ChIP-seq or ATAC-seq peaks
- Single-Nucleotide Polymorphisms (SNPs)
- Genes and gene components (TSS, exons, introns, etc)
- Non-coding DNA annotation (promoters, enhancers)
- Any sequence or subsequence (*e.g.* a chromosome)
- Chromosomes

---
![](images/genomic-intervals/chrom-sizes.png)

---
## Terminology

**coordinate system**: a set of sequence names and lengths

What constraints make two intervals comparable?

An interval comparison assumes a homology statement. A shared coordinate system is a prerequisite for interval comparison.

</div>

---
## What can be represented as an interval?

---
## Protein domains

![](images/genomic-intervals/ets-domains.jpeg)

---
## What can be represented as an interval?

---
# Key point

<span class="fragment">
Because of the linear nature of DNA and RNA, many biological entities can be conceptualized as genomic intervals.
</span>
<span class="fragment">Genomic intervals are often a simplified abstraction of genomic sequence.</span>

<span class="fragment">
$ \implies $
Interval operations are fundamental in genomics
</span>
---
## File formats for genomic intervals
- [BED](http://genome.ucsc.edu/FAQ/FAQformat.html#format1)
- [narrowPeak](http://genome.ucsc.edu/FAQ/FAQformat.html#format12)
- [broadPeak](http://genome.ucsc.edu/FAQ/FAQformat.html#format13)
- Genome Feature Format ([GTF](http://uswest.ensembl.org/info/website/upload/gff.html), [GFF2](http://gmod.org/wiki/GFF2), [GFF3](https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md))
  - Used for storing transcripts, exons, *etc.*
- Many custom variations!

---
### BED (Browser Extensible Data)

- originally designed by UCSC for the genome browser
- but, also often used for analysis
- a TSV file with one row per genomic region
- 3 required columns: `chrom`, `chromStart`, and `chromEnd`

Example of BED3 file:
```
chr1      15      250
chr2   38252    49502
```

---
## BED

- 9 additional optional fields:
	- name: name of region
	- score: range 0-1000
	- strand: `+`, `-`, or `.`
	- thickStart
	- thickEnd
	- itemRgb: display color (e.g. `55,0,0`)
	- blockCount
	- blockSizes
	- blockStarts

---
## BED

- BED3 is a BED with the first 3 columns
- BED6 is `chrom`, `chromStart`, and `chromEnd`, `name`, `score`, `strand`

---
## narrowPeak

- 10 Columns
- First 6 match BED: `chrom`, `chromStart`, and `chromEnd`, `name`, `score`, `strand`
- Final 4: `signalValue`, `pValue`, `qValue`, `peak`
- `peak`: Point-source (summit); offset from chromStart.

```
chr1    9356548 9356648 .       0       .       182     5.0945  -1  50
chr1    9358722 9358822 .       0       .       91      4.6052  -1  40
chr1    9361082 9361182 .       0       .       182     9.2103  -1  75
```

---
##  Related file formats
- [Wiggle](http://genome.ucsc.edu/goldenPath/help/wiggle.html) - continuous-valued data with constant size
	- Fixed-step - regular spacing between
	- Variable-step - variable spacing between
- BedGraph - continuous-valued data with variable size

---
## Wiggle (fixed step)

Continuous-valued data with constant size and regular spacing

Format:
```
fixedStep  chrom=chrN start=position  step=stepInterval [span=windowSize]
dataValue1
dataValue2
...
```

Example
```
fixedStep chrom=chr3 start=400601 step=100
11
22
33
```

---
## Wiggle (variable step)

Continuous-valued data with constant size and variable spacing

Format:
```
variableStep  chrom=chrN [span=windowSize]
  chromStartA  dataValueA
  chromStartB  dataValueB
  ... etc ...  ... etc ...
```

Example:
```
variableStep chrom=chr2
300701 12.5
300702 12.5
300703 12.5
300704 12.5
300705 12.5 
```

---
## BedGraph

Continuous-valued data with variable size

```
track type=bedGraph
chr19 49302000 49302300 -1.0
chr19 49302300 49302600 -0.75
chr19 49302600 49302900 -0.50
chr19 49302900 49303200 -0.25
chr19 49303200 49303500 0.0
chr19 49303500 49303800 0.25
```

---
## Why so many formats?

Sparsity. Wide vs tall structure on the data

- Widest: fixed-step wiggle
- Intermediate: variable-step wiggle
- Tallest: bedGraph
</div>

<div class="fragment">
Processing level. Raw signal vs. summary signal vs. regions
</div>

---
## Terminology

- **intersection**
- **union**
- **overlap**

---
### Biological questions, interval computations

- Is SNP $s$ located in a promoter element?
- What is the nearest gene downstream of SNP $s$?
- What is the average distance to genes for SNP set $S$?

---
### Biological questions, interval computations

- Which TF binding sites overlap region $r$?
- Are binding sites of TF $y$ enriched in enhancers?
- How many sites of TF $y$ are within 10 kb of gene $z$?

---
### Biological questions, interval computations

- What is the *shape* of signal across interval set $I$?
- What is the frequency of annotations for interval set $I$?
- How many aligned reads mapped to gene $g$?

---
# Overlap

Does $A$ overlap with $B$?

$  A \overset{?}\cap B > 0$

![](images/genomic-intervals/intervals-intersect.svg)

---
![](images/genomic-intervals/intervals-intersect.svg)

$A_{start} = 5$   $A_{end} = 25$   $B_{start} = 23$   $B_{end} = 30$

$cond 1: A_{end} \overset{?}> B_{start}$

<span class="fragment">
$cond 2: A_{start} \overset{?}< B_{end}$
</span>

---
Assume:
1. Intervals are well-formed ($X_{start} < X_{end}$)
2. Intervals are ordered by start ($A_{start} < B_{start}$)

$A_{start} < B_{start}  $ (order)  
$B_{start} < B_{end}  $  (well-formed)  
$\implies A_{start} < B_{end}$ (trans. property)

Then:

if $ A_{end} > B_{start}$

$ \implies A \cap B > 0$

---
# Counting overlap

![](images/genomic-intervals/intervals-intersect.svg)

$ overlap \overset{?}= A_{end} - B_{start} = 2$

$ = (A_{end} - B_{start}) - min(0, (A_{end} - B_{end})) = 2$

</span>

---
# Intersecting region

$ R = A \cap B $  
  
$R_{start} = max(A_{start}, B_{start})$  
$R_{end} = min(A_{end}, B_{end})$

---
# Union

$ R = A \cup B $  
  
if $ A \overset{?}\cap B$  
$R_{start} = min(A_{start}, B_{start})$  
$R_{end} = max(A_{end}, B_{end})$

---
###

If we make some assumptions, like pre-computing order, we can do things more efficiently.

Beware of containment

</div>

---
### Interval set overlap counting

Given:

A database intervals: $i_1,i_2,...\in\textbf{I}$,   
A query interval: $q$,

Count the number of intervals in $\textbf{I}$ overlapping $q$.

Application: What is the RNA-seq expression level of gene $g$?

<pre class="fragment"><code data-line-numbers="1-9">def overlap(query, i):
	return query.start < i.end && query.end > i.start

def count_overlaps(query, intervals_set):
	n_overlaps = 0
	for i in interval_set:
		if overlap(query, i):
			n_overlaps = n_overlaps +1
	return n_overlaps
</code></pre>

---

<pre><code data-line-numbers="6">def overlap(query, i):
	return query.start < i.end && query.end > i.start

def count_overlaps(query, intervals_set):
	n_overlaps = 0
	for i in interval_set:
		if overlap(query, i):
			n_overlaps = n_overlaps + 1
	return n_overlaps

</code></pre>

$O(N)$  
$N$ = Number of intervals in database

We may have billions of intervals.

Can we do better?

---
## Binary search

```Python
def binary_search(array, query, low, high):
	pivot = int((low + high)/2)
	if (low >= high):
		return False
	if (array[pivot] = query):
		return pivot
	if (array[pivot] > query):
		return binary_search(array, query, low, pivot-1)
	else:
		return binary_search(array, query, pivot+1, high)
```
> Binary search is one of the most powerful ideas in algorithm design... After only one comparison you can discard one half of all the names in the book. Twenty comparisons suffice to find any name in the million-name Manhattan phone book. -*Steven Skiena, The Algorithm Design Manual*
---
## Binary search

- Sequential search: $O(N)$; Binary search: $O(logN)$
- Requires sorted input (sequential search does not)

---

An $O(logN)$ option. Organize database as two lists:

![](images/genomic-intervals/subject-query-annotated.svg)

1. Binary search $q_{end}$ vs $\textbf{I}_{start}$  
2. Step backward until no more overlaps

![](images/genomic-intervals/binary-search.svg)
</div>

---

Say I have a series of regions in a BED file. Can I do a binary search on the file?

<div class="fragment">
Binary search also requires random access.

> Random access means the ability to access any element directly. This is contrast to *sequential access*, where I can only access elements in sequence.
</div>

---
### Binary search trees

<div class="mermaid">
graph TB;
    A((8))-->B((3))
    A-->C((10))
    B-->D((2))
    B-->E((6))
    C-->F((9))
    C-->G((12))
    E-->H((5))
    E-->I((7))
</div>

- Random access and array *insertion* are enemies.
- Insertion in a sorted array is $O(N)$.
- Binary search **trees** have $O(logN)$ search *and*  insertion.
- They do not require random access for binary search.

---

Hidden assumption: No containment

![](images/genomic-intervals/subject-query-containment.svg)

![](images/genomic-intervals/binary-search-fail.svg)

<div class="fragment">
How can we handle containment?
</div>

---
## Advanced data structures

- B trees
- R trees
- Nested containment lists
- Augmented interval lists

---
### B trees

- [B Tree at Wikipedia](https://en.wikipedia.org/wiki/B-tree)
- Generalization of a binary search tree:
	- Nodes have $m$ children
	- Node values are *boundaries* for subtrees
	- Boundaries provide overlap constraints

![](images/genomic-intervals/B-tree.svg)

---
### R trees - [Guttman 1984](https://doi.org/10.1145/971697.602266)
- Generalization of a B-tree beyond integer keys 
- Used by [BEDtools](https://bedtools.readthedocs.io/en/latest/)
- A flavor of [k-d trees](https://en.wikipedia.org/wiki/K-d_tree)
- Annotates tree nodes with a *minimum bounding rectangle* of elements. A query that does not intersect the bounding rectangle will not intersect any child element.

---

---
### Nested containment list
<a href="https://doi.org/10.1093/bioinformatics/btl647">Alekseyenko & Lee 2017</a>

- Used by Bioconductor GenomicRanges class
![](images/genomic-intervals/nclist.png)

---

Augmented Interval List - <a href="https://doi.org/10.1093/bioinformatics/btz407">Feng et al 2019</a>

![](images/genomic-intervals/subject-query-containment.svg)

![](images/genomic-intervals/binary-search-maxE.svg)

---
### Intro to [Bioconductor Genomic Ranges](https://bioconductor.org/packages/devel/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesIntroduction.html)

```R

if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("GenomicRanges")

library(GenomicRanges) # Load package

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

## GRanges object with 10 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        GC
##        <Rle> <IRanges>  <Rle> | <integer> <numeric>
##   a     chr1   101-111      - |         1  1.000000
##   b     chr2   102-112      + |         2  0.888889
##   c     chr2   103-113      + |         3  0.777778
##   d     chr2   104-114      * |         4  0.666667
##   e     chr1   105-115      * |         5  0.555556
##   f     chr1   106-116      + |         6  0.444444
##   g     chr3   107-117      + |         7  0.333333
##   h     chr3   108-118      + |         8  0.222222
##   i     chr3   109-119      - |         9  0.111111
##   j     chr3   110-120      - |        10  0.000000
##   -------
##   seqinfo: 3 sequences from an unspecified genome; no seqlengths
```

---

### Interval operations

```R
seqnames(gr)

## factor-Rle of length 10 with 4 runs
##   Lengths:    1    3    2    4
##   Values : chr1 chr2 chr1 chr3
## Levels(3): chr1 chr2 chr3

ranges(gr)

## IRanges object with 10 ranges and 0 metadata columns:
##         start       end     width
##     <integer> <integer> <integer>
##   a       101       111        11
##   b       102       112        11
##   c       103       113        11
##   d       104       114        11
##   e       105       115        11
##   f       106       116        11
##   g       107       117        11
##   h       108       118        11
##   i       109       119        11
##   j       110       120        11

sp <- split(gr, rep(1:2, each=5))
sp

## GRangesList object of length 2:
## $`1`
## GRanges object with 5 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        GC
##        <Rle> <IRanges>  <Rle> | <integer> <numeric>
##   a     chr1   101-111      - |         1  1.000000
##   b     chr2   102-112      + |         2  0.888889
##   c     chr2   103-113      + |         3  0.777778
##   d     chr2   104-114      * |         4  0.666667
##   e     chr1   105-115      * |         5  0.555556
##   -------
##   seqinfo: 3 sequences from an unspecified genome
## 
## $`2`
## GRanges object with 5 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        GC
##        <Rle> <IRanges>  <Rle> | <integer> <numeric>
##   f     chr1   106-116      + |         6  0.444444
##   g     chr3   107-117      + |         7  0.333333
##   h     chr3   108-118      + |         8  0.222222
##   i     chr3   109-119      - |         9  0.111111
##   j     chr3   110-120      - |        10  0.000000
##   -------
##   seqinfo: 3 sequences from an unspecified genome

flank(g, 10)

## GRanges object with 4 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        GC
##        <Rle> <IRanges>  <Rle> | <integer> <numeric>
##   a     chr1   112-121      - |         1  1.000000
##   b     chr2    92-101      + |         2  0.888889
##   c     chr2    93-102      + |         3  0.777778
##   j     chr3   121-130      - |        10  0.000000
##   -------
##   seqinfo: 3 sequences from an unspecified genome

shift(g, 5)

## GRanges object with 4 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        GC
##        <Rle> <IRanges>  <Rle> | <integer> <numeric>
##   a     chr1   106-116      - |         1  1.000000
##   b     chr2   107-117      + |         2  0.888889
##   c     chr2   108-118      + |         3  0.777778
##   j     chr3   115-125      - |        10  0.000000
##   -------
##   seqinfo: 3 sequences from an unspecified genome

resize(g, 30)

## GRanges object with 4 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        GC
##        <Rle> <IRanges>  <Rle> | <integer> <numeric>
##   a     chr1    82-111      - |         1  1.000000
##   b     chr2   102-131      + |         2  0.888889
##   c     chr2   103-132      + |         3  0.777778
##   j     chr3    91-120      - |        10  0.000000
##   -------
##   seqinfo: 3 sequences from an unspecified genome

g2 <- head(gr, n=2)
union(g, g2)

## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1   101-111      -
##   [2]     chr2   102-113      +
##   [3]     chr3   110-120      -
##   -------
##   seqinfo: 3 sequences from an unspecified genome

intersect(g, g2)

## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1   101-111      -
##   [2]     chr2   102-112      +
##   -------
##   seqinfo: 3 sequences from an unspecified genome

gr[2:3]  # subset

## GRanges object with 2 ranges and 2 metadata columns:
##     seqnames    ranges strand |     score        GC
##        <Rle> <IRanges>  <Rle> | <integer> <numeric>
##   b     chr2   102-112      + |         2  0.888889
##   c     chr2   103-113      + |         3  0.777778
##   -------
##   seqinfo: 3 sequences from an unspecified genome
```

---
### Region set enrichment analysis

Given a query region set (say, from a ChIP-seq experiment), what other, previously published region sets are most similar?

---
### Interval similarity metrics
- Overlap count
<li class="fragment">Jaccard index $\frac{A \cap B}{A \cup B}$</li>
<li class="fragment">Fisher's Exact Test (LOLA, <a href="https://doi.org/10.1093/bioinformatics/btv612">Sheffield and Bock 2016</a>)</li>

---
<img src="images/genomic-intervals/07-lola1.svg" />
---
<img src="images/genomic-intervals/08-lola2.svg" />
---
<img src="images/genomic-intervals/09-test.svg" />
---
<img src="images/genomic-intervals/10-results.svg" />

---

---
<img src="images/genomic-intervals/subject-query-one-vs-many-loop.svg">
<div class="fragment">
<img src="images/genomic-intervals/subject-query-one-vs-many-integrate.svg">
</div>

---
### B+ trees

- A B tree but nodes hold keys pointing to a linked list of values.
- Used by [GIGGLE](https://www.nature.com/articles/nmeth.4556)
- See [B+ Tree at Wikipedia](https://en.wikipedia.org/wiki/B%2B_tree)

![](images/genomic-intervals/Bplustree.png)
---
### Integrated Genome Database

- Introduced by [Feng and Sheffield 2021](https://doi.org/10.1093/bioinformatics/btaa1062)
- Uses a linear binning of an array

![](images/genomic-intervals/igd.jpg)

---
# Conclusion

- Genomic interval data is growing astronomically
- Genomic intervals are an extremely useful abstraction
- Diverse biological questions can be asked using intervals
- Data structures and algorithms enable large-scale analysis
- New methods for genomic intervals are being developed