---
title: "A Guide to using BiFET"
author: "Ahrim Youn"
date: '`r Sys.Date()`'
output:
rmarkdown::html_document:
toc: yes
toc_depth: 6
pdf_document:
toc: yes
vignette: >
%\VignetteIndexEntry{"A Guide to using BiFET"}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
## Introduction
`BiFET` (Bias-free Footprint Enrichment Test) aims to identify TFs whose
footprints are over-represented in target regions compared to background
regions after correcting for differences in read counts and GC contents between
target and background regions, where regions represent ATAC-seq or DNAse-seq
peaks. In this example, we used BiFET to detect TFs associated with
cell-specific open chromatin regions by analyzing ATAC-seq data from
human PBMCs *(Ucar, et al., 2017)* and pancreatic islets
*(Khetan, et al., 2017)*.
In order to use `BiFET`, the user will need an input file with
two `GRanges` objects:
* The first object `GRpeaks` represents a matrix of ATAC-seq or DNase-seq peaks in GRanges class where each row
represents the location of each peak. In addition there must be
3 metadata columns called "reads" (representing read counts in each peak),
"GC" (representing the GC content), and lastly "peaktype" which designates
each peak as either ("target","background","no").
* The second object `GRmotif` represents footprint calls from any footprint algorithms in GRanges class where each row represents the location
of each PWM occurrence. The footprint
calls in the forward strand and those in the backward strand from the same
PWM are not differentiated. The row names of GRmotif are the
motif IDs (e.g. MA01371 STAT1).
## 1. Obtaining a Peak File
```{r, eval=TRUE, echo=TRUE}
# Load necessary libraries
suppressPackageStartupMessages(library(BiFET))
suppressPackageStartupMessages(library(GenomicRanges))
library(BiFET)
library(GenomicRanges)
peak_file <- system.file("extdata", "input_peak_motif.Rdata",
package = "BiFET")
load(peak_file)
# Display the first few rows and columns of the peak file
head(GRpeaks)
```
In this example, the file *input_peak_motif.Rdata* was obtained
as follows: we used ATAC-seq data obtained from five human PBMC
*(Ucar, et al., 2017)* and five human islet samples *(Khetan, et al., 2017)*
and called peaks using **MACS version 2.1.0** *(Zhang, et al., 2008)* with
parameters **-nomodel -f BAMPE**. The peak sets from all samples were merged
to generate one consensus peak set (N = 57,108 peaks) by using
**R package DiffBind_2.2.5**. *(Ross-Innes, et al., 2012)*, where only the peaks
present at least in any two samples were included in the analysis.
We used the **summits** option to re-center each peak around the point of
greatest enrichment and obtained consensus peaks of same width **(200bp)**.
Out of these consensus peaks, we defined regions that are specifically
accessible in PBMC samples as regions where at least 4 PBMC samples have
a peak, whereas none of the islet samples have a peak
(n=4106 peaks; these regions are used as target regions in this example).
Similarly, we defined islet-specific peaks as those that were called as
a peak in at least 4 islet samples but none in any of the
PBMC samples (n=12886 peaks). The rest of the peaks excluding the
PBMC/islet-specific peaks were used as the background
(i.e., non-specific) peaks in our analyses (n=40116 peaks). For each peak,
GC content was obtained using peak annotation program
**annotatePeaks.pl** from the **HOMER** software *(Heinz. et al., 2010)*.
## 2. Obtaining a Matrix of Footprint Calls
```{r, eval=TRUE, echo=TRUE}
# Display the first few rows and columns of the motif file
head(GRmotif)
```
In this example, TF footprints were called using **PIQ** algorithm
*(Sherwood, et al., 2014)* using the pooled islet samples and pooled PBMC
samples to increase the detection power for TF footprints. We used only the
TF footprints that have a purity score greater than 0.9. The example file
contains footprint calls for only five PWMs from the
JASPAR database to reduce computing time.
## 3. Calculating enrichment p-value
The function `calculate_enrich_p` calculates p-value testing if footprints
of a TF are over-represented in the target set of peaks compared to the
background set of peaks correcting for the bias arising from the imbalance
of GC-content and read counts between target and background set. The function
requires two GRanges objects as input parameters: `GRpeaks` and `GRmotif`.
```{r, eval=TRUE, echo=TRUE}
# call the function “calculate_enrich_p†to return a list of
# parameter alpha_k, enrichment p values from BiFET algorithm
# and enrichment p values from the hypergeometric test :
result <- calculate_enrich_p(GRpeaks, GRmotif)
head(result)
```
## Citing BiFET
BiFET: A Bias-free Transcription Factor Footprint Enrichment Test, Ahrim Youn, Eladio J Marquez, Nathan Lawlor, Michael L Stitzel, Duygu Ucar, bioRxiv 2018:324277 (https://www.biorxiv.org/content/early/2018/05/16/324277)
---
## References
1. Khetan, S., et al. Chromatin accessibility profiling uncovers genetic-and
T2D disease state-associated changes in cis-regulatory element use in human
islets. bioRxiv 2017:192922.
2. Ross-Innes, C.S., et al. Differential oestrogen receptor binding is
associated with clinical outcome in breast cancer.
Nature 2012;481(7381):389-393.
3. Ucar, D., et al. The chromatin accessibility signature of human immune
aging stems from CD8+ T cells. Journal of
Experimental Medicine 2017;jem 20170416.
4. Sherwood, et al. Discovery of Directional and Nondirectional Pioneer
Transcription Factors by Modeling DNase Profile Magnitude and Shape.
Nature Biotechnology 32, no. 2 (February 2014): 171.
5. Heinz, et al. Simple Combinations of Lineage-Determining Transcription
Factors Prime Cis-Regulatory Elements Required for Macrophage and B Cell
Identities. Molecular Cell 38, no. 4 (May 28, 2010): 576–89.
6. Zhang, et al. Model-Based Analysis of ChIP-Seq (MACS).
Genome Biology 9 (2008): R137.
7. Ross-Innes, et al. Differential Oestrogen Receptor Binding Is Associated
with Clinical Outcome in Breast Cancer. Nature 481, no. 7381
(January 19, 2012): 389–93.