---
title: "ATACseqTFEA Guide"
author: "Jianhong Ou"
bibliography: bibliography.bib
csl: nature.csl
vignette: >
  %\VignetteIndexEntry{ATACseqTFEA Vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    toc_depth: 4
---

```{r, echo=FALSE, results="hide", warning=FALSE, message=FALSE}
suppressPackageStartupMessages({
  library(ATACseqTFEA)
  library(BSgenome.Drerio.UCSC.danRer10)
  library(Rsamtools)
  library(ATACseqQC)
})
knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.width=5, fig.height=3.5)
```

# Introduction

ATAC-seq, an assay for Transposase-Accessible Chromatin using sequencing, is a
widely used technique for chromatin accessibility analysis. 
Detecting differential activation of transcription factors between two 
different experiment conditions provides the possibility of decoding the 
key factors in a phenotype. Lots of tools have been developed to detect 
the differential activity of TFs (DATFs) for different groups of samples. 
Those tools can be divided into two groups. One group detects DATFs from 
differential accessibility analysis, such as MEME[@bailey2006meme],
HOMER[@heinz2010simple], enrichr[@chen2013enrichr], and ChEA[@lachmann2010chea].
Another group finds the DATFs by enrichment tests, such as
BiFET[@youn2019bifet], diffTF[@berest2019quantification], and TFEA[@rubin2020transcription].
For single-cell ATAC-seq analysis, Signac and chromVar are widely used tools.

# Motivation

All of these tools detect the DATF by only considering the open status of 
chromatin. None of them take the TF footprint into count. The open status 
provides the possibility of TF can bind to that position. 
The TF footprint by ATAC-seq shows the status of TF bindings. 

To help researchers quickly assess the differential activity of hundreds of TFs
by detecting the difference in TF footprint via enrichment
score[@subramanian2005gene],
we have developed the _ATACseqTFEA_ package. 
The _ATACseqTFEA_ package is a robust and reliable computational tool to
identify the key regulators responding to a phenotype. 

![schematic diagram of ATACseqTFEA](ATACseqTFEA.png)

# Quick start

Here is an example using _ATACseqTFEA_ with a subset of ATAC-seq data.

## Installation
First, install _ATACseqTFEA_ and other packages required to run 
the examples.
Please note that the example dataset used here is from zebrafish. 
To run an analysis with a dataset from a different species or different assembly, 
please install the corresponding Bsgenome and "TxDb".
For example, to analyze mouse data aligned to "mm10", 
please install "BSgenome.Mmusculus.UCSC.mm10", 
and "TxDb.Mmusculus.UCSC.mm10.knownGene". 
You can also generate a TxDb object by 
functions `makeTxDbFromGFF` from a local "gff" file,
or `makeTxDbFromUCSC`, `makeTxDbFromBiomart`, and `makeTxDbFromEnsembl`, 
from online resources in the _GenomicFeatures_ package.

```{r, eval=FALSE}
library(BiocManager)
BiocManager::install(c("ATACseqTFEA",
                       "ATACseqQC",
                       "Rsamtools",
                       "BSgenome.Drerio.UCSC.danRer10",
                       "TxDb.Drerio.UCSC.danRer10.refGene"))
```

## Load library

```{r}
library(ATACseqTFEA)
library(BSgenome.Drerio.UCSC.danRer10) ## for binding sites search
library(ATACseqQC) ## for footprint
```

## Prepare binding sites

To do TFEA, there are two inputs, the binding sites, and the change ranks.
To get the binding sites, the _ATACseqTFEA_ package provides the 
`prepareBindingSites` function. Users can also try to get the binding sites
list by other tools such as "fimo"[@grant2011fimo].

The `prepareBindingSites` function request a cluster of position weight matrix 
(PWM) of TF motifs. _ATACseqTFEA_ prepared a merged `PWMatrixList` for 
405 motifs. The `PWMatrixList` is a collection of jasper2018, jolma2013 and
cisbp_1.02 from package motifDB (v 1.28.0) and merged by distance smaller than
1e-9 calculated by MotIV::motifDistances function (v 1.42.0). 
The merged motifs were exported by motifStack (v 1.30.0).

```{r}
motifs <- readRDS(system.file("extdata", "PWMatrixList.rds",
                               package="ATACseqTFEA"))
```

The `best_curated_Human` is a list of TF motifs
downloaded from [TFEA github](https://github.com/Dowell-Lab/TFEA)[@rubin2020transcription].
There are 1279 human motifs in the data set.

```{r}
motifs_human <- readRDS(system.file("extdata", "best_curated_Human.rds",
                                    package="ATACseqTFEA"))
```

Another list of non-redundant TF motifs are also available by downloading
the data from [DeepSTARR](https://github.com/bernardo-de-almeida/motif-clustering)[@de2021deepstarr]. There are 6502 motifs in the data set.
```{r}
MotifsSTARR <- readRDS(system.file("extdata", "cluster_PWMs.rds",
                                      package="ATACseqTFEA"))
```

To scan the binding sites along a genome, a `BSgenome` object is required by
the `prepareBindingSites` function.

```{r}
# for test run, we use a subset of data within chr1:5000-100000
# for real data, use the merged peaklist as grange input.
# Drerio is the short-link of BSgenome.Drerio.UCSC.danRer10
seqlev <- "chr1" 
bindingSites <- 
  prepareBindingSites(motifs, Drerio, seqlev,
                      grange=GRanges("chr1", IRanges(5000, 100000)),
                      p.cutoff = 5e-05)#set higher p.cutoff to get more sites.
```


## TFEA

The correct insertion site is the key to the enrichment analysis of TF binding
sites. The parameter `positive` and `negative` in the function of `TFEA`
are used to shift the 5' ends of the reads to the correct insertion positions.
However, this shift does not consider the soft clip of the reads.
The best way to generate correct shifted bam files is using
ATACseqQC::shiftGAlignmentsList[@ou2018atacseqqc] for paired-end or 
shiftGAlignments for single-end of the bam file.
The samples must be at least biologically duplicated for the one-step `TFEA`
function.

```{r}
bamExp <- system.file("extdata",
                      c("KD.shift.rep1.bam",
                        "KD.shift.rep2.bam"),
                      package="ATACseqTFEA")
bamCtl <- system.file("extdata",
                      c("WT.shift.rep1.bam",
                        "WT.shift.rep2.bam"),
                      package="ATACseqTFEA")
res <- TFEA(bamExp, bamCtl, bindingSites=bindingSites,
            positive=0, negative=0) # the bam files were shifted reads
```

## View results

The results will be saved in a `TFEAresults` object. We will use multiple
functions to present the results.
The `plotES` function will return a `ggplot` object for single TF input and
no `outfolder` is defined.
The `ESvolcanoplot` function will provide an overview of all the TFs enrichment.
And we can borrow the `factorFootprints` function from `ATACseqQC` package
to view the footprints of one TF.

```{r}
TF <- "Tal1::Gata1"
## volcanoplot
ESvolcanoplot(TFEAresults=res, TFnameToShow=TF)
### plot enrichment score for one TF
plotES(res, TF=TF, outfolder=NA)
## footprint
sigs <- factorFootprints(c(bamCtl, bamExp), 
                         pfm = as.matrix(motifs[[TF]]),
                         bindingSites = getBindingSites(res, TF=TF),
                         seqlev = seqlev, genome = Drerio,
                         upstream = 100, downstream = 100,
                         group = rep(c("WT", "KD"), each=2))
## export the results into a csv file
write.csv(res$resultsTable, tempfile(fileext = ".csv"), 
          row.names=FALSE)
```

The command-line scripts are available at `extdata` named as `sample_scripts.R`. 

# Do TFEA step by step.

The one-step `TFEA` is a function containing multiple steps, which include:

1. count the reads in binding sites, proximal region, and distal region;
2. filter the binding site not open;
3. normalize the count number by the width of the count region;
4. calculate the binding scores and weight the binding scores by open scores;
5. differential analysis by limma for the binding score
6. filter the differential results by P-value and fold change
7. TF enrichment analysis

If you want to tune the parameters, it will be much better to do it step by step
to avoid repeating the computation for the same step.
Here are the details for each step.

## Counting reads

We will count the insertion site in binding sites, proximal and distal regions
by counting the 5' ends of the reads in a shifted bam file.
Here we suggest keeping the `proximal` and `distal` the same value.

```{r}
# prepare the counting region
exbs <- expandBindingSites(bindingSites=bindingSites,
                           proximal=40,
                           distal=40,
                           gap=10)
## count reads by 5'ends
counts <- count5ends(bam=c(bamExp, bamCtl),
                     positive=0L, negative=0L,
                     bindingSites = bindingSites,
                     bindingSitesWithGap=exbs$bindingSitesWithGap,
                     bindingSitesWithProximal=exbs$bindingSitesWithProximal,
                     bindingSitesWithProximalAndGap=
                         exbs$bindingSitesWithProximalAndGap,
                     bindingSitesWithDistal=exbs$bindingSitesWithDistal)
```

## Filter the counts

We filter the binding sites by at least there is 1 reads in proximal region.
Users may want to try filter the sites by more stringent criteria such as
"proximalRegion>1".

```{r}
colnames(counts)
counts <- eventsFilter(counts, "proximalRegion>0")
```

## Normalize the counts by width of count region

We will normalize the counts to count per base (CPB).

```{r}
counts <- countsNormalization(counts, proximal=40, distal=40)
```

## Get weighted binding scores

Here we use the open score to weight the binding score. Users can also define
the weight for binding score via parameter `weight` in
the function `getWeightedBindingScore`.

```{r}
counts <- getWeightedBindingScore(counts)
```


## Differential analysis

Here we use `DBscore`, which borrows the power of the `limma` package,
to do differential binding analysis.

```{r}
design <- cbind(CTL=1, EXPvsCTL=c(1, 1, 0, 0))
counts <- DBscore(counts, design=design, coef="EXPvsCTL")
```

## Filter the DB results

We can filter the binding results to decrease the data size by
the function `eventsFilter`.
For the sample data, we skip this step.

## TF enrichment analysis

Last, we use the function `doTFEA` to get the enrichment scores.

```{r}
res <- doTFEA(counts)
res
plotES(res, TF=TF, outfolder=NA) ## will show same figure as above one
```

# SessionInfo

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


# References
