chip-seq.Rmd

---
title: 'Bioinformatics for Big Omics Data: ChIP-seq data analysis'
author: "Raphael Gottardo"
date: "February 12, 2015"
output:
  ioslides_presentation:
    fig_caption: yes
    fig_retina: 1
    keep_md: yes
    smaller: yes
---

## Setting up some options

Let's first turn on the cache for increased performance and improved styling
```{r, cache=FALSE}
# Set some global knitr options
library("knitr")
opts_chunk$set(tidy=TRUE, tidy.opts=list(blank=FALSE, width.cutoff=60), cache=TRUE, messages=FALSE)
```


## Outline

During this lecture, we will discuss another application of next generation sequencing to 
study DNA-binding proteins using a technique called chromatin immunoprecipitation. You should read the following papers:

1. Park, P. J. ChIP-seq: advantages and challenges of a maturing technology. Nature Reviews Genetics 10, 669–680 (2009). 
2. Zhang, X. et al. PICS: probabilistic inference for ChIP-seq. Biometrics 67, 151–163 (2011).
3. Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, R137 (2008).

We will focus on the following topics:
- High throughput sequencing and ChIP-Seq
- Computational and statistical challenges

## Experimental protocol

<img src="http://www.nature.com/nrg/journal/v10/n10/images/nrg2641-f1.jpg" width=400>


Source: [Park Nat Rev Gen (2009)](http://www.nature.com/nrg/journal/v10/n10/full/nrg2641.html)

## ChIP-seq vs chip-chip

For a good comparison of chip-chip and chip-seq, have a look at [Park's paper](http://www.nature.com/nrg/journal/v10/n10/fig_tab/nrg2641_T1.html).

## Motivation for data analysis

<img src="http://www.nature.com/nrg/journal/v10/n10/images/nrg2641-f5.jpg" width=300>

Source: [Park Nat Rev Gen (2009)](http://www.nature.com/nrg/journal/v10/n10/fig_tab/nrg2641_F5.html)

Given the figure on the left, how would you analyze ChIP-seq data?

## ChIP-Seq summary

Couple ChIP with NGS: 

- A typical ChIP-Seq experiment generates millions of short reads
- Read lengths are in the order of 50-150bps, mostly single-end, but some paired-end
- Because of chromatin, antibodies and alignment biases, a control sample is still recommended
- Numerous types of controls: Input DNA, Mock IP, Naked DNA, etc. 

## Doing some analysis in R

There exist multiple packages for data analysis in R. Here we will explore one of the 
simpler packages available on Bioconductor. As usual, you first need to install the package

```{r eval=FALSE}
library(BiocInstaller)
biocLite("chipseq")
```

and then load the package.
```{r}
library(chipseq)
```

## Getting stated with the chipseq package

The `chipseq` package contains a small dataset, `cstest`, that we will use here. 
The `cstest` dataset contains data for three chromosomes from Solexa lanes, one from a CTCF mouse ChIP-Seq, and one from a GFP mouse ChIP-Seq. The raw reads were already aligned to the mouse genome and read into R using `ShortReads`. The resulting data are stored in a `GRanges` object.


```{r}
data(cstest)
```

let's look what's inside

```{r}
cstest
```

## Extending reads

The first step is to reconstruct the orginal DNA fragments by extending the shortreads

```{r}
# Estimate the fragment length
fraglen <- estimate.mean.fraglen(cstest$ctcf)
fraglen
```

```{r}
table(width(cstest$ctcf))
```

```{r}
# Now we can extend the fragments
ctcf.ext <- resize(cstest$ctcf, width = 180)
```

## Coverage

Next we calculate the coverage

```{r}
ctcf_cov <- coverage(ctcf.ext)
```


## Islands

The regions of interest are contiguous segments of non-zero coverage, also known as islands.


```{r}
ctcf_peaks <- slice(ctcf_cov, lower = 50)
```

## Islands

For each island, we can compute the number of reads in the island, and the maximum coverage depth within that island.

```{r}
viewSums(ctcf_peaks)
viewMaxs(ctcf_peaks)
```

## Finding peaks

```{r}
ctcf_summary <- peakSummary(ctcf_peaks)
# Subset chr10
ctcf_summary_chr10 <- ctcf_summary["chr10"]
# Order by score (max)
ctcf_summary_chr10 <- ctcf_summary_chr10[order(-ctcf_summary_chr10$max),]
ctcf_summary_chr10
```

## Validating peaks

After identifying peaks, it might be useful to see if the peaks contain the expected binding site for the transcription factor in question. We first need to download the Mmusculus genome sequence, to extract the relevant sequences.

```{r eval=FALSE}
biocLite("BSgenome.Mmusculus.UCSC.mm10")
```

```{r}
library(BSgenome.Mmusculus.UCSC.mm10)
Mmusculus[["chr10"]]
# Create views for our peaks
seq_views <- Views(Mmusculus[["chr10"]],start=ctcf_summary_chr10$maxpos-100, end=ctcf_summary_chr10$maxpos+100)
```

## Validating peaks

Let's then scan for our expected motif:

```{r}
# CTCF logo
CTCF_logo <- DNAString("GGCG")
# matchPWM
matchPattern(CTCF_logo, seq_views, max.mismatch=0)
matchPattern(rev(CTCF_logo), seq_views, max.mismatch=0)
```

## Accounting for biases

- Biases in ChIP-seq experiment can arise due to technical or biological reasons. Such biases can significantly affect the ranking of the ChIP-seq peaks. Biases can happen because of non-specificity of the antibody, DNA accessibility (open chromatin), or even GC content

- As we have discussed earlier, one potential way to control for these sources of bias is to use a control

- Enrichment scores can be adjusted for non-specific binding by computing a relative score (fold change) between the IP and control samples

- By swapping the IP/Control samples, we can get an estimate of the number of false positives for a given threshold, and thus compute an estimate of the FDR

- This requires an analysis of the control sample

## Peak calling

- Here we've seen a fairly simple approach to ChIP-seq analysis. There exist many other (more elaborated) approaches:
     - MACS: Zhang, Y. et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9, R137 (2008).
     - PICS: Zhang, X. et al. PICS: probabilistic inference for ChIP-seq. Biometrics 67, 151–163 (2011).
     - cisGenomes: Ji, H. et al. An integrated software system for analyzing ChIP-chip and ChIP-seq data. Nature Biotechnology 26, 1293–1300 (2008).
     - etc, etc
     
Although these approaches use different statistical methods or models, they all follow the same basic idea. They mostly differ in the way they correct for background binding and their FDR calculation. 


## PICS - Our approach

- Probabilistic model
    - Model binding events
    - Use prior information (Fragment length distribution)
    - Measures of uncertainty
- Bidirectional reads
- Estimation of missing reads
- PICS: Probabilistic Inference of ChIP-Seq data

## PICS

- Divide the genomic reads into regions with “enough” F/R reads
- Scan the genome every 10 bps  with a sliding window of size 150 bps
- Minimum number of F reads on the left and R reads on the right
- Merge overlapping regions 
- N disjoint regions 
- Model each region separately

## PICS model

<img src="./Images/PICS1.png" width=600>

## PICS model

<img src="./Images/PICS2.png" width=600>

## PICS model (prior)


<img src="./Images/PICS3.png" width=600>

## PICS prior

<img src="./Images/PICS4.png" width=600>

## PICS prior

<img src="./Images/PICS5.png" width=600>

## PICS prior

<img src="./Images/PICS6.png" width=600>

## PICS prior

<img src="./Images/PICS7.png" width=600>

## Handling missing reads

- Use an idea similar to that of McLachlan and Jones (1998) for grouped and truncated data

- We introduce additional missing data, corresponding to the number of F/R missing reads, and the missing reads themselves

- Assume the number of missing reads follows a multinomial distribution over the possible intervals

- Closed formed E-steps, and the EM remains essentially the same

## PICS and missing reads

<img src="./Images/PICS8.png" width=600>

## PICS and missing reads

<img src="./Images/PICS9.png" width=600>

## ChIP, histones and nucleosome positioning

ChIP-seq can also be used to study histone marks and nucleosome positioning. See for example, the following papers:

- Zhang, X., Robertson, G., Woo, S., Hoffman, B. G. & Gottardo, R. Probabilistic inference for nucleosome positioning with MNase-based or sonicated short-read data. PLoS ONE 7, e32095 (2012).

- Woo, S., Zhang, X., Sauteraud, R., Robert, F. & Gottardo, R. PING 2.0: an R/Bioconductor package for nucleosome positioning using next-generation sequencing data. Bioinformatics 29, 2049–2050 (2013).