---
title: "ATACseqQC Guide"
author: "Jianhong Ou, Jun Yu, Michelle Kelliher, Lucio Castilla, Nathan Lawson, Lihua Julie Zhu"
date: "`r BiocStyle::doc_date()`"
package: "`r BiocStyle::pkg_ver('ATACseqQC')`"
bibliography: bibliography.bib
csl: nature.csl
vignette: >
  %\VignetteIndexEntry{ATACseqQC Vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document
---

```{r, echo=FALSE, results="hide", warning=FALSE, message=FALSE}
suppressPackageStartupMessages({
  library(ATACseqQC)
  library(ChIPpeakAnno)
  library(BSgenome.Hsapiens.UCSC.hg19)
  library(TxDb.Hsapiens.UCSC.hg19.knownGene)
  library(phastCons100way.UCSC.hg19)
  library(MotifDb)
})
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
```

# Introduction

Assay for Transposase-Accessible Chromatin using sequencing 
(ATAC-seq) is an alternative or complementary technique to MNase-seq DNase, 
and FAIRE-seq for chromatin accessibility analysis. The results 
obtained from ATAC-seq are similar to those from DNase and FAIRE-seq. 
ATAC-seq is gaining popularity because it does not require cross-linking, 
has higher signal to noise ratio, requires a much smaller amount of biological 
material and is faster and easier to perform, compared to other techniques [@buenrostro2013transposition]. 

To help researchers quickly assess the quality of  ATAC-seq data, we 
have developed ATACseqQC package for easily making diagnostic plots 
following the published guidelines [@buenrostro2013transposition]. 
In addition, it has functions to preprocess ATACseq data for 
subsequent peak calling.

# Quick start

Here is an example to using ATACseqQC with a subset of published ATAC-seq 
data[@buenrostro2013transposition]. Currently, only bam input file 
format is supported.

```{r}
## load the library
library(ATACseqQC)
## input is bamFile
bamfile <- system.file("extdata", "GL1.bam", 
                        package="ATACseqQC", mustWork=TRUE)
bamfile.labels <- gsub(".bam", "", basename(bamfile))
```

## Fragment size distribution

First, there should be a large proportion of reads with less 
than 100 bp, which represent the nucleosome-free region. 
Second, the fragment size distribution should have a clear 
periodicity, which is evident in the inset figure, indicative of 
nucleosome occupation (present in integer multiples).

```{r}
## generate fragement size distribution
fragSize <- fragSizeDist(bamfile, bamfile.labels)
```

## Nucleosome positioning

### Adjust the read start sites

Tn5 transposase has been shown to bind as a dimer and 
insert two adaptors separated by 9 bp [@adey2010rapid].

Therefore, for downstream analysis, such as peak-calling and footprinting,
all reads in input bamfile need to be shifted. 
The function `shiftGAlignmentsList` can be used to shift the reads. 
By default, all reads aligning to the positive strand are offset by +4bp, 
and all reads aligning to the negative strand are offset by -5bp [@buenrostro2013transposition]. 

The adjusted reads will be written into a new bamfile for peak calling or 
footprinting.

```{r}
## bamfile tags
tags <- c("AS", "XN", "XM", "XO", "XG", "NM", "MD", "YS", "YT")
## files will be output into outPath
outPath <- "splited"
dir.create(outPath)
## shift the bam file by the 5'ends
library(BSgenome.Hsapiens.UCSC.hg19)
seqlev <- "chr1" ## subsample data for quick run
which <- as(seqinfo(Hsapiens)[seqlev], "GRanges")
gal <- readBamFile(bamfile, tag=tags, which=which, asMates=TRUE)
gal1 <- shiftGAlignmentsList(gal)
shiftedBamfile <- file.path(outPath, "shifted.bam")
export(gal1, shiftedBamfile)
```

### Split reads

The shifted reads will be split into different bins, namely
nucleosome free, mononucleosome, dinucleosome, and trinucleosome.
Shifted reads that do not fit into any of the above bins will
be discarded. Splitting reads is a time-consuming step 
because we are using random forest to classify the fragments 
based on fragment length, GC content and conservation scores 
[@chen2013danpos].

By default, we assign the top 10% of short reads (reads below 100_bp) 
as nucleosome-free regions and the top 10% of intermediate length reads
as (reads between 180 and 247 bp) mononucleosome. 
This serves as the training set to classify the rest of the fragments
using random forest. The number of the tree will be set to 2 times 
of square root of the length of the training set.

```{r}
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txs <- transcripts(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(phastCons100way.UCSC.hg19)
## run program for chromosome 1 only
txs <- txs[seqnames(txs) %in% "chr1"]
genome <- Hsapiens
## split the reads into NucleosomeFree, mononucleosome, 
## dinucleosome and trinucleosome.
objs <- splitGAlignmentsByCut(gal1, txs=txs, genome=genome,
                              conservation=phastCons100way.UCSC.hg19)
```

Save the binned alignments into bam files.

```{r}
null <- writeListOfGAlignments(objs, outPath)
## list the files generated by splitBam.
dir(outPath)
```

You can also do shift, split and save bams in one step by calling `splitBam`.

```{r eval=FALSE}
objs <- splitBam(bamfile, tags=tags, outPath=outPath,
                 txs=txs, genome=genome,
                 conservation=phastCons100way.UCSC.hg19)
```

### Heatmap and coverage curve for nucleosome positions

By averaging the signal across all active TSSs, we should observe that 
nucleosome-free fragments are enriched at the TSS, 
whereas the nucleosome-bound fragments should be enriched both upstream 
and downstream 
of the active TSS and display characteristic phasing of upstream and 
downstream nucleosomes. Because ATAC-seq reads are concentrated at regions of 
open chromatin, users should see a strong nucleosome signal at the +1 
nucleosome that decreases at the +2, +3 and +4 nucleosomes. 

```{r fig.height=4, fig.width=4}
library(ChIPpeakAnno)
bamfiles <- file.path(outPath,
                     c("NucleosomeFree.bam",
                     "mononucleosome.bam",
                     "dinucleosome.bam",
                     "trinucleosome.bam"))
## Plot the cumulative percentage tag allocation in NucleosomeFree 
## and mononucleosome bams.
cumulativePercentage(bamfiles[1:2], as(seqinfo(Hsapiens)["chr1"], "GRanges"))
```
```{r fig.height=8, fig.width=4}
TSS <- promoters(txs, upstream=0, downstream=1)
TSS <- unique(TSS)
## estimate the library size for normalization
(librarySize <- estLibSize(bamfiles))
## calculate the signals around TSSs.
NTILE <- 101
dws <- ups <- 1010
sigs <- enrichedFragments(bamfiles, TSS=TSS,
                          librarySize=librarySize,
                          seqlev=seqlev,
                          TSS.filter=0.5,
                          n.tile = NTILE,
                          upstream = ups,
                          downstream = dws)
## log2 transformed signals
names(sigs) <- gsub(".bam", "", basename(names(sigs)))
sigs.log2 <- lapply(sigs, function(.ele) log2(.ele+1))
#plot heatmap
featureAlignedHeatmap(sigs.log2, reCenterPeaks(TSS, width=ups+dws),
                      zeroAt=.5, n.tile=NTILE)
```
```{r fig.show="hide"}
## get signals normalized for nucleosome-free and nucleosome-bound regions.
out <- featureAlignedDistribution(sigs, 
                                  reCenterPeaks(TSS, width=ups+dws),
                                  zeroAt=.5, n.tile=NTILE, type="l")
```
```{r}
## rescale the nucleosome-free and nucleosome signals to 0~1
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
out <- apply(out, 2, range01)
matplot(out, type="l", xaxt="n", 
        xlab="Position (bp)", 
        ylab="Fraction of signal")
axis(1, at=seq(0, 100, by=10)+1, 
     labels=c("-1K", seq(-800, 800, by=200), "1K"), las=3)
abline(v=seq(0, 100, by=10)+1, lty=2, col="gray")
```

## Footprints

ATAC-seq footprints infer factor occupancy genome-wide. The `factorFootprints` 
function uses `matchPWM` to predict the binding sites using the input position
weight matrix (PWM). 
Then it calculates and plots the accumulated coverage for those binding sites
to show the status of the occupancy genome-wide.
Unlike CENTIPEDE[@pique2011accurate], the footprints generated here 
do not take into consideration the conservation (PhyloP). 
`factorFootprints` function could also accept the possible 
binding sites as a GRanges object.

```{r}
## foot prints
library(MotifDb)
CTCF <- query(MotifDb, c("CTCF"))
CTCF <- as.list(CTCF)
print(CTCF[[1]], digits=2)
factorFootprints(shiftedBamfile, pfm=CTCF[[1]], 
                 genome=genome,
                 min.score="95%", seqlev=seqlev,
                 upstream=100, downstream=100)
```


Here is the CTCF footprints for the full dataset.
![CTCF footprints](CTCFfootprints.png)

# Session Info
```{r sessionInfo}
sessionInfo()
```

# References
