%\VignetteIndexEntry{DiffBind: Differential binding analysis of ChIP-Seq peak data}
%\VignettePackage{DiffBind}
%\VignetteEngine{utils::Sweave}

\documentclass{article}

<<style, eval=TRUE, echo=FALSE, results=tex>>=
BiocStyle::latex()
savewd = getwd()
@ 

\newcommand{\DBA}{\Biocpkg{DiffBind}~}
\newcommand{\DiffBind}{\Biocpkg{DiffBind}~}
\newcommand{\edgeR}{\Biocpkg{edgeR}~}
\newcommand{\DESeq}{\Biocpkg{DESeq}~}

\newcommand{\reft}[1]{Table \ref{tab:#1}}
\newcommand{\reff}[1]{Figure \ref{fig:#1}}
\newcommand{\refsf}[1]{Figure \ref{fig:#1}}
\newcommand{\code}[1]{{\Rcode{#1}}}

\begin{document}
\SweaveOpts{concordance=TRUE}
\SweaveOpts{prefix.string=\Sexpr{savewd}/DiffBind}

\newcommand{\exitem}[3]{\item \Rcode{\textbackslash#1\{#2\}} #3 \csname#1\endcsname{#2}.}

\title{\Biocpkg{DiffBind}: Differential binding analysis of ChIP-Seq peak data}
\author{Rory Stark \\ \Rcode{rory.stark@cruk.cam.ac.uk}
\and
Gord Brown \\ \Rcode{gdbzork@gmail.com} 
\and University of Cambridge\\Cancer Research UK - Cambridge Institute
}

\date{Edited: February 12 2015; Compiled: \today}

\maketitle

\tableofcontents

\section{Introduction}

This document offers an introduction and overview of the
\R~Bioconductor package  \DBA, which
provides functions for processing ChIP-Seq data enriched for genomic loci where specific protein/DNA binding occurs, including peak sets identified by ChIP-Seq peak callers and aligned sequence read datasets. It is designed to work with multiple peak sets simultaneously, representing different ChIP experiments (antibodies, transcription factor and/or histone marks, experimental conditions, replicates) as well as managing the results of multiple peak callers. 

The primary emphasis of the package is on identifying sites that are differentially bound between two sample groups. 
It includes functions to support the processing of peak sets, including overlapping and merging peak sets, counting sequencing reads overlapping intervals in peak sets, and identifying statistically significantly differentially bound sites based on evidence of binding affinity (measured by differences in read densities). To this end it uses statistical routines developed in an RNA-Seq context (primarily the Bioconductor packages \edgeR and \DESeq). Additionally, the package builds on \R graphics routines to provide a set of standardized plots to aid in binding analysis.

This guide includes a brief overview of the processing flow, followed by three sections of examples: the first focusing on the core task of obtaining differentially bound sites based on affinity data, the second working through the main plotting routines, and the third revisiting occupancy data (peak calls) in more detail, as well as comparing the results of an occupancy-based analysis with an affinity-based one. Finally, some technical aspects of the how these analyses are accomplished are detailed.

\section{Processing overview}

\DBA works primarily with peaksets, which are sets of genomic intervals representing candidate protein binding sites. Each interval consists of a chromosome, a start and end position, and usually a score of some type indicating confidence in, or strength of, the peak. Associated with each peakset are metadata relating to the experiment from which the peakset was derived. Additionally, files containing mapped sequencing reads (BAM//BED) can be associated with each peakset (one for the ChIP data, and optionally another representing a control dataset).

Generally, processing data with \DBA involves five phases:
\begin{enumerate}
\item
\Rcode{Reading in peaksets}:
The first step is to read in a set of peaksets and associated metadata. Peaksets are derived either from ChIP-Seq peak callers, such as MACS (\cite{zhang2008model}), or using some other criterion (e.g. all the promoter regions in a genome). The easiest way to read in peaksets is using a comma-separated value (csv) sample sheet with one line for each peakset. (Spreadsheets in Excel\textsuperscript{\textregistered} format, with a \code{.xls} or \code{.xlsx} suffix, are also accepted.)  A single experiment can have more than one associated peakset, e.g. if multiple peak callers are used for comparison purposes, and hence have more than one line in the sample sheet. Once the peaksets are read in, a merging function finds all overlapping peaks and derives a single set of unique genomic intervals covering all the supplied peaks.
\item
\Rcode{Occupancy analysis}:
Peaksets, especially those generated by peak callers, provide an insight into the potential \emph{occupancy} of the protein being ChIPed for at specific genomic loci.  After the peaksets have been loaded, it can be useful to perform some exploratory plotting to determine how these occupancy maps agree with each other, e.g. between experimental replicates (re-doing the ChIP under the same conditions), between different peak callers on the same experiment, and within groups of samples representing a common experimental condition. \DBA provides functions to enable overlaps to be examined, as well as functions to determine how well similar samples cluster together. Beyond quality control, the product of an occupancy analysis may be a \emph{consensus peakset}, representing an overall set of candidate binding sites to be used in further analysis.
\item
\Rcode{Counting reads}:
Once a consensus peakset has been derived, \DBA can use the supplied sequence read files to count how many reads overlap each interval for each unique sample. 
The result of this is a \emph{binding affinity matrix} containing a (normalized) read count for each sample at every potential binding site. 
With this matrix, the samples can be re-clustered using affinity, rather than occupancy, data. 
The binding affinity matrix is used for QC plotting as well as for subsequent differential analysis.
\item
\Rcode{Differential binding affinity analysis}:
The core functionality of \DBA is the differential binding affinity analysis, which enables binding sites to be identified that are statistically significantly differentially bound between sample groups. 
To accomplish this, first a contrast (or contrasts) is established, dividing the samples into groups to be compared. 
Next the core analysis routines are executed, by default using \edgeR. 
This will assign a p-value and FDR to each candidate binding site indicating the significance of their being differentially bound.
\item
\Rcode{Plotting and reporting}:
Once one or more contrasts have been run, \DBA provides a number of functions for reporting and plotting the results. 
MA plots give an overview of the results of the analysis, while correlation heatmaps and PCA plots show how the groups cluster based on differentially bound sites. 
Boxplots show the distribution of reads within differentially bound sites corresponding to whether they gain or lose affinity between the two sample groups. 
A reporting mechanism enables differentially bound sites to be extracted for further processing, such as annotation and/or pathway analysis.
\end{enumerate}

\section{Example: obtaining differentially bound sites}

This section offers a quick example of how to use \DBA to identify significantly differentially bound sites using affinity (read count) data.

The dataset for this example consists of ChIPs against the transcription factor ERa using five breast cancer cell lines (\cite{ross2012differential}). 
Three of these cell lines are responsive to tamoxifen, while two others are resistant to tamoxifen treatment. 
There are at least two replicates for each of the cell lines, with one cell line having three replicates, for a total of eleven sequenced libraries. 
Note that this experiment includes two types of MCF7 cells: the regular tamoxifen responsive line as well as MCF7 cells specially treated with tamoxifen until a tamoxifen resistant cell line is obtained.
For each sample, we have one peakset originally derived using the MACS peak caller  (\cite{zhang2008model}), for a total of eleven peaksets. 
Note that to save space in the package, only data for chromosome 18 is used. The metadata and peak data are available in the \Rcode{extra} subdirectory of the \DBA package directory; you can make this your working directory by entering:

<<eval=TRUE,echo=FALSE,results=hide>>=
tmp =  tempfile(as.character(Sys.getpid()))
pdf(tmp)
savewarn = options("warn")
options(warn=-1)
@


<<eval=TRUE,echo=TRUE,results=verbatim>>=
library(DiffBind)
setwd(system.file("extra", package="DiffBind"))
@

Obtaining the sites significantly differentially bound (DB) between the samples that respond to tamoxifen and those that are resistant can be done in a five-step script:

<<eval=FALSE,echo=TRUE,results=hide>>=
print(savewd)
tamoxifen = dba(sampleSheet="tamoxifen.csv")
tamoxifen = dba.count(tamoxifen)
tamoxifen = dba.contrast(tamoxifen, categories=DBA_CONDITION)
tamoxifen = dba.analyze(tamoxifen)
tamoxifen.DB = dba.report(tamoxifen)
@

The following subsections describe these steps in more detail


\subsection{Reading in the peaksets}

The easiest way to set up an experiment to analyze is with a sample sheet. The sample sheet can be a \Robject{dataframe}, or it can be read directly from a \Rcode{csv} file. Here is the example sample sheet read into a \Rcode{dataframe} form a \Rcode{csv} file:

<<sampSheet, eval=TRUE,echo=TRUE,results=verbatim>>=
samples = read.csv(file.path(system.file("extra", package="DiffBind"),
                                         "tamoxifen.csv"))
samples
@

The peaksets are read in using the following \DBA function:

<<dbaContruct, eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen = dba(sampleSheet="tamoxifen.csv")
@

The result is a \Rclass{DBA} object; the metadata associated with this object can be displayed simply as follows:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen
@


This shows how many peaks are in each peakset, as well as (in the first line) total number of unique peaks after merging overlapping ones
(\Sexpr{nrow(tamoxifen$allvectors)})
and the default binding matrix of 11 samples by the
\Sexpr{nrow(tamoxifen$vectors)}
sites that overlap in at least two of the samples. This object is available for loading using \Rcode{data(tamoxifen\_peaks)}.

Using only this peak caller data, a correlation heatmap can be generated which gives an initial clustering of the samples using the cross-correlations of each row of the binding matrix:

<<tamox_occ_corhm, fig=TRUE, include=FALSE, width=9, height=9 >>=
plot(tamoxifen)
@
\incfig{DiffBind-tamox_occ_corhm}{.66\textwidth}{Correlation heatmap, using occupancy (peak caller score) data.}
  {Generated by: \Rcode{plot(tamoxifen)}; can also be generated by: \Rcode{dba.plotHeatmap(tamoxifen).}}
  
The resulting plot (Figure~\ref{DiffBind-tamox_occ_corhm}) shows that while the replicates for each cell line cluster together appropriately, the cell lines do not cluster into groups corresponding to those that are responsive (MCF7, T47D, and ZR75) vs. those resistant (BT474 and MCF7r) to tamoxifen treatment. It also shows that the two most highly correlated cell lines are the two MCF7-based ones, even though they respond differently to tamoxifen treatment. 

\subsection{Counting reads}
The next step is to calculate a binding matrix with scores based on read counts for every sample (affinity scores), rather than confidence scores for only those peaks called in a specific sample (occupancy scores). These reads are obtained using the \Rcode{dba.count} function:\footnote{Note that due to space limitations the reads are not shipped with the package. See Section 8 for options to obtains the full dataset. Alternatively, you can get the result of the \Rcode{dba.count} call by loading the supplied \R object by invoking \Rcode{data(tamoxifen\_counts)}}

<<eval=FALSE,echo=TRUE,results=hide>>=
tamoxifen = dba.count(tamoxifen, minOverlap=3) 
@

<<eval=TRUE,echo=FALSE,results=hide>>=
data(tamoxifen_counts)
@
<<tamox_aff_corhm, fig=TRUE, include=FALSE, width=9, height=9 >>=
plot(tamoxifen)
@
\incfig{DiffBind-tamox_aff_corhm}{.66\textwidth}{Correlation heatmap, using affinity (read count) data.}
  {Generated by: \Rcode{tamoxifen = dba.count(tamoxifen)}; can also be generated by: \Rcode{dba.plotHeatmap(tamoxifen)}}

If you do not have the raw reads available to you, this object is available loading using \Rcode{data(tamoxifen\_counts)}.
The \Rcode{dba.count} call plots a new correlation heatmap based on the affinity scores, seen in Figure~\ref{DiffBind-tamox_aff_corhm}. While this shows overall higher correlation levels (evidenced by the shift in the scale), and a slightly different clustering, responsiveness to tamoxifen treatment does not appear to form a basis for clustering when using all of the affinity scores. (Note that the clustering can change based on what scoring metric is used; see Section 4.4 for more details).

\subsection{Establishing a contrast}
Before running the differential analysis, we need to tell \DBA which cell lines fall in which groups. This is done using the \Rcode{dba.contrast} function, as follows:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen = dba.contrast(tamoxifen, categories=DBA_CONDITION)
@

The uses the \emph{Condition} metadata (Responsive vs. Resistant) to set up a a contrast with 4 samples in the Resistant group and 7 samples in the Responsive group.\footnote{This step is actually optional: if the main analysis function \Rcode{dba.analyze} is invoked with no contrasts established, \DBA will set up a default set of contrasts automatically, which will include the one we are interested in.}

\subsection{Performing the differential analysis}
The main differential analysis function is invoked as follows:

<<tamox_sdb_corhm, fig=TRUE, include=FALSE, width=9, height=9 >>=
tamoxifen = dba.analyze(tamoxifen) 
@
\incfig{DiffBind-tamox_sdb_corhm}{.66\textwidth}{Correlation heatmap, using only significantly differentially bound sites.}
  {Generated by: \Rcode{tamoxifen = dba.analyze(tamoxifen)}; can also be generated by: \Rcode{dba.plotHeatmap(tamoxifen, contrast=1)}}

This will run an \edgeR analysis (see subsequent section discussing the technical details of the \edgeR analysis) using the default binding matrix.
Displaying the resultant \Rcode{DBA object} shows that
\Sexpr{as.numeric(levels(dba.show(tamoxifen,bContrasts=T)[1,5]))}
of the
\Sexpr{nrow(tamoxifen$allvectors)}
sites are identified as being significantly differentially bound (DB) using the default threshold of FDR <= 0.1

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen    
@ 

By default, \Rcode{dba.analyze} plots a correlation heatmap if it finds any significantly differentially bound sites, shown in Figure~\ref{DiffBind-tamox_sdb_corhm}. 
Using only the differentially bound sites, we now see that the four tamoxifen resistant samples (representing two cell lines) cluster together, although the tamoxifen-responsive MCF7 replicates cluster closer to them than to the other tamoxifen responsive samples. Comparing Figure~\ref{DiffBind-tamox_aff_corhm}, which uses all 
\Sexpr{nrow(tamoxifen$allvectors)} 
consensus binding sites, with Figure~\ref{DiffBind-tamox_sdb_corhm}, which uses only the 
\Sexpr{as.numeric(levels(dba.show(tamoxifen,bContrasts=T)[1,5]))} 
differentially bound sites, demonstrates how the differential binding analysis isolates sites that help distinguish between the Resistant and Responsive sample groups. 

\subsection{Retrieving the differentially bound sites}
The final step is to retrieve the differentially bound sites as follows:

<<eval=TRUE,echo=TRUE,results=hide>>=
tamoxifen.DB = dba.report(tamoxifen)
@

These are returned as a \Rclass{GRanges} object, appropriate for downstream processing:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen.DB
@


The value columns show the mean read concentration over all the samples (the default calculation uses log2 normalized ChIP read counts with control read counts subtracted) and the mean concentration over the first (Resistant) group and second (Responsive) group. The Fold column shows the difference in mean concentrations between the two groups (Conc\_Resistant - Conc\_Responsive), with a positive value indicating increased binding affinity in the Resistant group and a negative value indicating increased binding affinity in the Responsive group. The final two columns give confidence measures for identifying these sites as differentially bound, with a raw p-value and a multiple testing corrected FDR in the final column.

\section{Example: plotting}

Besides the correlation heatmaps automatically generated by the core functions, a number of other plots are available using the affinity data. 
This sections covers Venn diagrams, MA plots, PCA plots, Boxplots, and Heatmaps.

\subsection{Venn diagrams}

Venn diagrams are useful for examining overlaps between peaksets, particularly when determining how best to derive consensus peaksets for further analysis. Section 6.2, which discusses consensus peaksets, shows a number of Venn plots in context, and the help page for \Rcode{dba.plotVenn} has a number of additional examples.

\subsection{MA plots}
MA plots are a useful way to visualize the effect of normalization on data, as well as seeing which of the datapoints are being identified as differentially bound. 
An MA plot can be obtained for the resistant-responsive contrast as follows:

<<tamox_sdb_ma, fig=TRUE, include=FALSE, width=12, height=9 >>=
data(tamoxifen_analysis)
dba.plotMA(tamoxifen)
@
\incfig{DiffBind-tamox_sdb_ma}{.66\textwidth}{MA plot of Resistant-Responsive contrast.}
  {Sites identified as significantly differentially bound shown in red. Generated by: \Rcode{dba.plotMA(tamoxifen)}}

The plot is shown in  Figure~\ref{DiffBind-tamox_sdb_ma}. Each point represents a binding site, with points in red representing sites identified as differentially bound. 
The plot shows how the differentially bound sites appear to have an absolute log fold difference of at least 2. 
This same data can also be shown with the concentrations of each sample groups plotted against each other plot using \Rcode{dba.plotMA(tamoxifen, bXY=TRUE)}.

\subsection{PCA plots}
While the correlation heatmaps already seen are good for showing clustering, plots based on principal components analysis can be used to give a deeper insight into how samples are associated. 
A PCA plot corresponding to Figure~\ref{DiffBind-tamox_aff_corhm}, which includes normalized read counts for all the binding sites, can be obtained as follows:

<<tamox_aff_pca, fig=TRUE, include=FALSE, width=9, height=9 >>=
dba.plotPCA(tamoxifen,DBA_TISSUE,label=DBA_CONDITION)
@
\incfig{DiffBind-tamox_aff_pca}{.66\textwidth}{PCA plot using affinity data for all sites.}
  { Generated by: \Rcode{dba.plotPCA(tamoxifen,DBA\_TISSUE,label=DBA\_CONDITION)}}
  
The resulting plot (Figure~\ref{DiffBind-tamox_aff_pca}) shows the four Resistant samples (black) not separable from the Responsive samples (red) in either the first (horizontal) or the second (vertical) components when looking at all the binding sites.
 
A PCA plot using only the differentially bound sites (corresponding to Figure~\ref{DiffBind-tamox_sdb_corhm}), 
using an FDR threshold of 0.05, can be drawn as follows:

<<tamox_sdb_pca, fig=TRUE, include=FALSE, width=9, height=9 >>=
dba.plotPCA(tamoxifen, contrast=1,th=.05,label=DBA_TISSUE)
@
\incfig{DiffBind-tamox_sdb_pca}{.66\textwidth}{PCA plot using affinity data for only differentially bound sites.}
  { Generated by: \Rcode{dba.plotPCA(tamoxifen,contrast=1,th=.05,label=DBA\_TISSUE)}}
  
This plot (Figure~\ref{DiffBind-tamox_sdb_pca}) shows that the differential analysis identifies sites than can be used to separate the sample groups along both the first and second components.

The \Rcode{dba.plotPCA} function is customizable. 
For example, if you want to see where the replicates for each of the unique cell lines lies, type
\Rcode{dba.plotPCA(tamoxifen, attributes=c(DBA\_TISSUE,DBA\_CONDITION),label=DBA\_REPLICATE}. 
If your installation of \R supports 3D graphics using the \Rpackage{rgl} package, try \Rcode{dba.plotPCA(tamoxifen, b3D=T)}. 
Seeing the first three principal components can be a useful exploratory exercise.

\subsection{Boxplots}

Boxplots provide a way to view how read distributions differ between classes of binding sites. Consider the example, where the 
\Sexpr{as.numeric(levels(dba.show(tamoxifen,bContrasts=T)[1,5]))} 
differentially bound sites are identified. 
The MA plot (Figure~\ref{DiffBind-tamox_sdb_ma}) shows that these are not distributed evenly between those that increase binding affinity in the Responsive group vs. those that increase binding affinity in the Resistant groups. 
This can be seen quantitatively using the sites returned in the report:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
sum(tamoxifen.DB$Fold<0)
sum(tamoxifen.DB$Fold>0)
@

But how are reads distributed amongst the different classes of differentially bound sites and sample groups? These data can be more clearly seen using a boxplot:

<<tamox_sdb_box, fig=TRUE, include=FALSE, width=9, height=9 >>=
pvals = dba.plotBox(tamoxifen)
@
\incfig{DiffBind-tamox_sdb_box}{.66\textwidth}{Box plots of read distributions for significantly differentially bound (DB) sites.}
  { Tamoxifen resistant samples are shown in red, and responsive samples are shown in blue. Left two boxes show distribution of reads over all DB sites in the Resistant and Responsive groups; middle two boxes show distributions of reads in DB sites that increase in affinity in the Responsive group; last two boxes show distributions of reads in DB sites that increase in affinity in the Resistant group. Generated by: \Rcode{dba.plotBox(tamoxifen)}}
  

The default plot (Figure~\ref{DiffBind-tamox_sdb_box}) shows in the first two boxes that amongst differentially bound sites overall, the Responsive samples have a somewhat higher mean read concentration. The next two boxes show the distribution of reads in differentially bound sites that exhibit increased affinity in the Responsive samples, while the final two boxes show the distribution of reads in differentially bound sites that exhibit increased affinity in the Resistant samples.

\Rcode{dba.plotBox} returns a matrix of p-values (computed using a two-sided Wilcoxon `Mann-Whitney' test, paired where appropriate) indicating which of these distributions are significantly different from another distribution.

<<eval=TRUE,echo=TRUE,results=verbatim>>=
pvals
@

The significance of the overall difference in distribution of concentrations amongst the differentially bound sites in the two groups is shown to be p-value=\Sexpr{sprintf("%1.6e",pvals[1,2])}, while those between the Resistant and Responsive groups in the individual cases (increased in Responsive or Resistant) have p-values computed as 
\Sexpr{sprintf("%1.6e",pvals[3,4])} and \Sexpr{sprintf("%1.6e",pvals[5,6])}.

\subsection{Heatmaps}

\DBA provides two types of heatmaps. This first, correlation heatmaps, we have already seen. 
For example, the heatmap shown in Figure~\ref{DiffBind-tamox_aff_corhm} can be generated as follows:
<<eval=TRUE,echo=TRUE,results=hide>>=
corvals = dba.plotHeatmap(tamoxifen)
@

The effect of different scoring methods (normalization) can be examined in these plots by setting the \Rcode{score} parameter to a different value.  The default value, \Rcode{DBA\_SCORE\_TMM\_MINUS\_FULL}, uses the TMM normalization procedure from \edgeR, with control reads subtracted first and using the full library size (total reads in library). Another scoring method is to use RPKM fold (RPKM of the ChIP reads divided by RPKM of the control reads; a correlation heatmap for all the data using this scoring method can be obtained by typing \Rcode{dba.plotHeatmap(tamoxifen, score=DBA\_SCORE\_RPKM\_FOLD)}.

Another way to view the patterns of binding affinity directly in the differentially bound sites is via a binding affinity heatmap. This can be plotted for the example case as follows:

<<tamox_sdb_hm, fig=TRUE, include=FALSE, width=9, height=9 >>=
corvals = dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE)
@
\incfig{DiffBind-tamox_sdb_hm}{.66\textwidth}{Binding affinity heatmap showing affinities for differentially bound sites.}
{Samples cluster first by whether they are responsive to tamoxifen treatment, then by cell line. Clusters of binding sites show distinct patterns of affinity levels. Generated by: \Rcode{dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE)}}

Figure~\ref{DiffBind-tamox_sdb_hm} shows the affinities and clustering of the differentially bound sites (rows), as well as the sample clustering (columns). This plot can be tweaked to get more contrast, for example by using row-scaling \Rcode{dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE, scale=row)}.

\section{Example: differential binding analysis using  a blocking factor}
     
The previous example showed how to perform a differential binding analysis using a single factor with two values; that is, finding the significantly differentially bound sites between two sets of samples. This section extends the example by including a second factor, potentially with multiple values, that represents a confounding condition. Examples of experiments where it is appropriate to use a blocking factor include ones where there are potential batch effects, with samples from the two conditions prepared together, or a matched design (e.g. matched normal and tumor pairs, where the primary factor of interest is to discover sites consistently differentially bound between normal and tumor samples. In the current example, the confounding effect we want to control for is the presence of two sets of samples, one tamoxifen responsive and one resistant, that are both derived from the same MCF7 cell line.

In the previous analysis, the two MCF7-derived cell lines tended to cluster together. While the differential binding analysis was able to identify sites that could be used to separate the resistance from the responsive samples, the confounding effect of the common ancestry could still be seen even when considering only the significantly differentially bound sites (Figure~\ref{DiffBind-tamox_aff_corhm}).

Using the generalized linear modelling (GLM) functionality included in \edgeR and \DESeq, the confounding factor can be explicitly modeled. This is done by specifying a blocking factor to \Rcode{dba.contrast}. There are a number of ways to specify this factor. If it is encapsulated in a piece of metadata (eg. DBA\_REPLICATE, or DBA\_TREATMENT etc.), simply specifying the metadata field is sufficient. In the current case, there is no specific metadata field that captures the factor we want to block (although an unused metadata field, such as DBA\_TREATMENT, could be used to specify this factor). An alternate way of specifying the confounded samples is to use a mask:

<<eval=TRUE,echo=TRUE,results=hide>>=
data(tamoxifen_counts)
tamoxifen = dba.contrast(tamoxifen,categories=DBA_CONDITION,
                         block=tamoxifen$masks$MCF7)
@

Now when the analysis is run, it will be run using both the single-factor comparison as well as fitting a linear model with the second, blocking factor, for comparison:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen = dba.analyze(tamoxifen)
tamoxifen
@

This indicates that where the standard, single-factor \edgeR analysis identifies 
\Sexpr{as.numeric(levels(dba.show(tamoxifen,bContrasts=T)[1,9]))} 
differentially bound sites, the analysis using the blocking factor finds 
\Sexpr{as.numeric(levels(dba.show(tamoxifen,bContrasts=T)[1,10]))} 
such sites. An MA plot shows how the analysis has changed:

<<tamox_block_ma, fig=TRUE, include=FALSE, width=12, height=9 >>=
dba.plotMA(tamoxifen,method=DBA_EDGER_BLOCK)
@
\incfig{DiffBind-tamox_block_ma}{.66\textwidth}{MA plot of Resistant-Responsive contrast, using MCF7 origin as a blocking factor.}
{Sites identified as significantly differentially bound shown in red. 
  Generated by: \Rcode{dba.plotMA(tamoxifen, method=DBA\_EDGER\_BLOCK)}}

The resulting plot is shown in Figure~\ref{DiffBind-tamox_block_ma}. Comparing this to Figure~\ref{DiffBind-tamox_sdb_ma}, at least two differences can be observed. The analysis has become more sensitive, with sites being identified as significantly differentially bound with lower magnitude fold changes (as low as twofold, as this plot is on a log2 scale). But it is not merely lowering a fold threshold: some sites with higher fold changes are no longer found to be significant. These were identified as significantly differentially bound in the earlier analysis because the confounding factor was not being modeled.

Consider the resulting separation and clustering using the newly discovered differentially bound sites:

<<tamox_block_corhm, fig=TRUE, include=FALSE, width=9, height=9 >>=
dba.plotHeatmap(tamoxifen,contrast=1,method=DBA_EDGER_BLOCK,
                attributes=c(DBA_TISSUE,DBA_CONDITION,DBA_REPLICATE))
@

\incfig{DiffBind-tamox_block_corhm}{.66\textwidth}{Correlation heatmap of using scores for significantly differentially bound sites for the Resistant-Responsive contrast, using MCF7 origin as a blocking factor.}
{Generated by: \Rcode{dba.plotHeatmap(tamoxifen, contrast=1, method=DBA\_EDGER\_BLOCK, attributes=c(DBA\_TISSUE,DBA\_CONDITION,DBA\_REPLICATE))}}
  
<<tamox_block_pca, fig=TRUE, include=FALSE, width=9, height=9 >>=
dba.plotPCA(tamoxifen,contrast=1,method=DBA_EDGER_BLOCK,
            attributes=DBA_CONDITION,label=DBA_TISSUE)
@
\incfig{DiffBind-tamox_block_pca}{.66\textwidth}{Plot of first two principal components, using scores for significantly differentially bound sites for the Resistant-Responsive contrast, using MCF7 origin as a blocking factor.}
{Generated by: \Rcode{dba.plotPCA(tamoxifen,contrast=1,method=DBA\_EDGER\_BLOCK, attributes=DBA\_CONDITION,label=DBA\_TISSUE)}}

Frequently, as more sites are included in these plots, the result is often worse clustering/separation along the grouping of primary interest. 
As the correlations used for the heatmap and the count scores used for the PCA plots do not take the blocking factor into account, the separation seen in these plots will not necessarily sharpen after completing a successful block analysis. 
Compare the resulting correlation heatmap (Figure~\ref{DiffBind-tamox_block_corhm} to Figure~\ref{DiffBind-tamox_sdb_corhm}).  
The PCA plot (Figure~\ref{DiffBind-tamox_block_pca}, compared to Figure~\ref{DiffBind-tamox_sdb_pca}), shows less distance between the two sample groups along the first component, as more binding sites common to the two MCF7 cell lines are able to be identified as being significantly differentially bound between the Resistant and Responsive sample groups.

It is also interesting to compare the performance of \edgeR with that of \DESeq and \Biocpkg{DESeq2} on this dataset:
<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen = dba.analyze(tamoxifen,method=DBA_ALL_METHODS)
tamoxifen
@

\DESeq has a much more conservative approach to the single factor analysis, identifying only
\Sexpr{as.numeric(levels(dba.show(tamoxifen,bContrasts=T)[1,11]))}
sites as differentially bound, while \Biocpkg{DESeq2} identifies a similar number of sites as \Biocpkg{edgeR}.
When modelling the confounding factor, the greater sensitivity results in many more sites being identified using all three packages. You can check this by looking at the identified sites using \Rcode{dba.report}, and performing MA, heatmap, and PCA plots.

We can also compare the sites identified using \Biocpkg{edgeR}, \Biocpkg{DESeq}, and \Biocpkg{DESeq2}. An easy way to do this is to use a special feature of the \Rcode{dba.report} function:

<<tamox_block_venn, fig=TRUE, include=FALSE, width=9, height=9 >>=
tam.block = dba.report(tamoxifen,method=DBA_ALL_BLOCK,bDB=TRUE,bAll=TRUE)
tam.block
dba.plotVenn(tam.block,1:3,label1="edgeR",label2="DESeq",label3="DESeq2")
@
\incfig{DiffBind-tamox_block_venn}{.66\textwidth}{Venn diagram showing overlap of differentially bound peaks identified using \Biocpkg{edgeR}, \Biocpkg{DESeq}, and \Biocpkg{DESeq2} to do a multi-factor (blocking) analysis.}
{Generated by plotting the result of: \Rcode{dba.plotVenn(tam.block,1:3,
label1="edgeR",label2="DESeq",label3="DESeq2")}}

The overlap is shown in Figure~\ref{DiffBind-tamox_block_venn}. Agreement is quite good, with \Biocpkg{DESeq2} identifying the fewest sites not identified by at least one other method.


\section{Example: occupancy analysis and overlaps}

In this section, we look at the tamoxifen resistance ER-binding dataset in some more detail, showing what a pure occupancy-based analysis would look like, and comparing it to the results obtained using the affinity data. For this we will start by re-loading the peaksets:

<<eval=TRUE,echo=TRUE,results=hide>>=
data(tamoxifen_peaks)
@

\subsection{Overlap rates}

One reason to do an occupancy-based analysis is to determine what candidate sites should be used in a subsequent affinity-based analysis. In the example so far, we took all sites that were identified in peaks in at least three of the eleven peaksets, reducing the number of sites from 
\Sexpr{nrow(tamoxifen$allvectors)}
overall to the 
\Sexpr{dba.overlap(tamoxifen,mode=DBA_OLAP_RATE)[3]}
sites used in the differential analysis.
We could have used a more stringent criterion, such as only taking sites identified in five or six of the peakset, or a less stringent one, such as including all 
\Sexpr{nrow(tamoxifen$allvectors)}
sites. In making the decision of what criteria to use many factors come into play, but it helps to get an idea of the rates at which the peaksets overlap 
(for more details on how overlaps are determined, see Section 7.2 on peak merging). 
A global overview can be obtained using the \Rcode{RATE} mode of the \Rcode{dba.overlap} function as follows:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
olap.rate = dba.overlap(tamoxifen,mode=DBA_OLAP_RATE)
olap.rate
@
The returned data in \Rcode{olap.rate} is a vector containing the number of peaks that appear in at least one, two, three, and so on up to all eleven peaksets.

These values can be plotted to show the overlap rate drop-off curve:

<<tamox_rate, fig=TRUE, include=FALSE, width=9, height=6 >>=
plot(olap.rate,type='b',ylab='# peaks', xlab='Overlap at least this many peaksets') 
@
\incfig{DiffBind-tamox_rate}{.66\textwidth}{Overlap rate plot.}
{Shows how the number of overlapping peaks decreases as the overlap criteria becomes more stringent. X axis shows the number of peaksets in which the site is identified, while the Y axis shows the number of overlapping sites. 
Generated by plotting the result of: \Rcode{dba.overlap(tamoxifen,mode=DBA\_OLAP\_RATE)}}

The rate plot is shown in Figure~\ref{DiffBind-tamox_rate}. These curves typically exhibit a roughly geometric drop-off, with the number of overlapping sites halving as the overlap criterion become stricter by one site. When the drop-off is extremely steep, this is an indication that the peaksets do not agree very well. For example, if there are replicates you expect to agree, there may be a problem with the experiment. In the current example, peak agreement is high and the curve exhibits a better than geometric drop-off.

\subsection{Deriving consensus peaksets}
When performing an overlap analysis, it is often the case that the overlap criteria are set stringently in order to lower noise and drive down false positives.
\footnote{It is less clear that limiting the potential binding sites in this way is appropriate when focusing on affinity data, as the differential binding analysis method will identify only sites that are significantly differentially bound, even if operating on peaksets that include incorrectly identified sites.} The presence of a peak in multiple peaksets is an indication that it is a "real" binding site, in the sense of being identifiable in a repeatable manner. The use of biological replicates (performing the ChIP multiple times), as in the tamoxifen dataset, can be used to guide derivation of a consensus peakset. Alternatively, an inexpensive but less powerful way to help accomplish this is to use multiple peak callers for each ChIP dataset and look for agreement between peak callers (\cite{li2011measuring}). 

Consider for example the standard (tamoxifen responsive) MCF7 cell line, represented by three replicates in this dataset. How well do the replicates agree on their peak calls? The overlap rate for just the positive MCF7 samples can be isolated using a \emph{sample mask}. A set of sample masks are automatically associated with a \Rcode{DBA object} in the \Rcode{\$masks} field:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
names(tamoxifen$masks)
@

Arbitrary masks can be generated using the \Rcode{dba.mask} function, or simply by specifying a vector of peakset numbers. In this case, a mask that isolates the MCF7 samples can be generated by combining to pre-defined masks (MCF7 and Responsive) and passed into the \Rcode{dba.overlap} function:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
dba.overlap(tamoxifen,tamoxifen$masks$MCF7 & tamoxifen$masks$Responsive,
            mode=DBA_OLAP_RATE)
@

There are 
\Sexpr{dba.overlap(tamoxifen,tamoxifen$masks$MCF7 & tamoxifen$masks$Responsive,mode=DBA_OLAP_RATE)[3]}
peaks (out of 
\Sexpr{dba.overlap(tamoxifen,tamoxifen$masks$MCF7 & tamoxifen$masks$Responsive,mode=DBA_OLAP_RATE)[1]})
identified in all three replicates. A finer grained view of the overlaps can be obtained with the \Rcode{dba.plotVenn} function:


<<tamox_mcf7_venn, fig=TRUE, include=FALSE, width=9, height=9 >>=
dba.plotVenn(tamoxifen, tamoxifen$masks$MCF7 & tamoxifen$masks$Responsive)
@
\incfig{DiffBind-tamox_mcf7_venn}{.66\textwidth}{Venn diagram showing how the ER peak calls for three replicates of responsive MCF7 cell line overlap.}
{Generated by plotting the result of: \Rcode{dba.venn(tamoxifen,tamoxifen\$masks\$MCF7 \& tamoxifen\$masks\$Responsive)}}

The resultant plot is shown as Figure~\ref{DiffBind-tamox_mcf7_venn}. This plot shows the 
\Sexpr{dba.overlap(tamoxifen,tamoxifen$masks$MCF7 & tamoxifen$masks$Responsive,mode=DBA_OLAP_RATE)[3]}
consensus peaks identified as common to all replicates, but further breaks down how the replicates relate to each other. 
The same can be done for each of the replicated cell line experiments, and rather than applying a global cutoff (3 of 11), each cell line could be dealt with individually in deriving a final peakset. 
A separate consensus peakset for each of the replicated sample types can be added to the DBA object using \Rcode{dba.peakset}:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen = dba.peakset(tamoxifen, consensus = c(DBA_TISSUE,DBA_CONDITION), 
            minOverlap=0.66)
tamoxifen
@

This adds a new consensus peakset for each set of samples that share the same Tissue and Condition values. The exact effect could be obtained by calling    
\Rcode{tamoxifen = dba.peakset(tamoxifen, consensus = -DBA\_REPLICATE)} on the original set of peaks; this tells \DiffBind to generate a consensus peakset for every set of samples that have identical metadata values \emph{except} the Replicate number.

From this, a new \Rclass{DBA} object can be generated consisting of only the five consensus peaksets (the \$Consensus mask filters peaksets previously formed using \Rcode{dba.peakset}) :

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen_consensus = dba(tamoxifen, mask = tamoxifen$masks$Consensus)
tamoxifen_consensus
@

Alternatively, a master consensus peakset could be generated, and reads counted, directly using \Rcode{dba.count}:
\Rcode{tamoxifen = dba.count(tamoxifen, peaks=tamoxifen\$masks\$Consensus)}

Finally, consider an analysis where we wished to treat all five MCF7 samples together to look for binding sites specific to that cell line irrespective of tamoxifen resistant/responsive status. We can create consensus peaksets for each cell type, and look at how the resultant peaks overlap (shown in Figure~\ref{DiffBind-tamox_lines_venn}):

<<tamox_lines_venn, fig=TRUE, include=FALSE, width=9, height=9 >>=
data(tamoxifen_peaks)
tamoxifen = dba.peakset(tamoxifen, consensus = DBA_TISSUE, minOverlap=0.66)
dba.plotVenn(tamoxifen, tamoxifen$masks$Consensus)
@
\incfig{DiffBind-tamox_lines_venn}{.66\textwidth}{Venn diagram showing how the consensus peaks for each cell type overlap.}
{Generated by plotting the result of: \Rcode{dba.venn(tamoxifen,tamoxifen\$masks\$Consensus)}}

\subsection{A complete occupancy analysis: identifying sites unique to a sample group}

Occupancy-based analysis, in addition to offering many ways of deriving consensus peaksets, can also be used to identify sites unique to a group of samples. This is analogous to, but not the same as, finding differentially bound sites. In these subsections, the two approaches are directly compared.

Returning to the original tamoxifen dataset:

<<eval=TRUE,echo=TRUE,results=hide>>=
data(tamoxifen_peaks)
@

We can derive consensus peaksets for the Resistant and Responsive groups. First we examine the overlap rates:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
dba.overlap(tamoxifen,tamoxifen$masks$Resistant,mode=DBA_OLAP_RATE)
dba.overlap(tamoxifen,tamoxifen$masks$Responsive,mode=DBA_OLAP_RATE)
@

Requiring that consensus peaks overlap in at least one third of the samples in each group results in
\Sexpr{dba.overlap(tamoxifen,tamoxifen$masks$Resistant,mode=DBA_OLAP_RATE)[2]}
sites for the Resistant group and 
\Sexpr{dba.overlap(tamoxifen,tamoxifen$masks$Responsive,mode=DBA_OLAP_RATE)[3]} 
sites for the Responsive group:

<<tamox_cons_venn, fig=TRUE, include=FALSE, width=9, height=9 >>=
tamoxifen = dba.peakset(tamoxifen, consensus =  DBA_CONDITION, minOverlap = 0.33)
dba.plotVenn(tamoxifen,tamoxifen$masks$Consensus)
@
\incfig{DiffBind-tamox_cons_venn}{.66\textwidth}{Venn diagram showing how the ER peak calls for two response groups overlap.}
{Generated by plotting the result of: \Rcode{dba.plotVenn(tamoxifen, tamoxifen\$masks\$Consensus)}}

Figure~\ref{DiffBind-tamox_cons_venn} shows that 
\Sexpr{nrow(dba.overlap(tamoxifen,tamoxifen$masks$Consensus,DataType=DBA_DATA_FRAME)$onlyA)} 
sites are unique to the Resistant group, and 
\Sexpr{nrow(dba.overlap(tamoxifen,tamoxifen$masks$Consensus,DataType=DBA_DATA_FRAME)$onlyB)} 
sites are unique to the Responsive group, with
\Sexpr{nrow(dba.overlap(tamoxifen,tamoxifen$masks$Consensus,DataType=DBA_DATA_FRAME)$inAll)} 
sites being identified in both groups (meaning in at least half the Resistant samples and at least three of the seven Responsive samples). 
If our primary interest is in finding binding sites that are different between the two groups, it may seem reasonable to consider the
\Sexpr{nrow(dba.overlap(tamoxifen,tamoxifen$masks$Consensus,DataType=DBA_DATA_FRAME)$inAll)} 
common sites to be uninteresting, and focus on the 
\Sexpr{nrow(dba.overlap(tamoxifen,tamoxifen$masks$Consensus,DataType=DBA_DATA_FRAME)$onlyA) + nrow(dba.overlap(tamoxifen,tamoxifen$masks$Consensus,DataType=DBA_DATA_FRAME)$onlyB)} 
sites that are unique to a specific group. These unique sites can be obtained using \Rcode{dba.overlap}:

<<eval=TRUE,echo=TRUE,results=hide>>=
tamoxifen.OL = dba.overlap(tamoxifen, tamoxifen$masks$Consensus)
@

The sites unique to the Resistant group are accessible in \Rcode{tamoxifen.OL\$onlyA}, with the Responsive-unique sites in \Rcode{tamoxifen.OL\$onlyB}:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen.OL$onlyA
tamoxifen.OL$onlyB
@


The scores associated with each site are derived from the peak caller confidence score, and are a measure of confidence in the peak call (occupancy), not a measure of how strong or distinct the peak is.

\subsection{Comparison of occupancy and affinity based analyses}

So how does this occupancy-based analysis compare to the previous affinity-based analysis? 

First, different criteria were used to select the overall consensus peakset. We can compare them to see how well they agree:

<<tamox_compare_venn, fig=TRUE, include=FALSE, width=9, height=9 >>=
tamoxifen = dba.peakset(tamoxifen,tamoxifen$masks$Consensus, 
            minOverlap=1,sampID="OL Consensus")
tamoxifen = dba.peakset(tamoxifen,!tamoxifen$masks$Consensus,
            minOverlap=3,sampID="Consensus_3") 
dba.plotVenn(tamoxifen,14:15)
@
\incfig{DiffBind-tamox_compare_venn}{.66\textwidth}{Venn diagram showing how the ER peak calls for two different ways of deriving consensus peaksets.}
{Generated by plotting the result of: \Rcode{dba.plotVenn(tamoxifen,14:15)}}

Figure~\ref{DiffBind-tamox_compare_venn} shows that the two sets agree on about 85\% of their sites, so the results should be directly comparable between the differing parameters used to establish the  consensus peaksets.
\footnote{Alternatively, we could re-run the analysis using the newly derived consensus peakset by passing it into the counting function: \Rcode{> tamoxifen = dba.count(tamoxifen, peaks = tamoxifen\$masks\$Consensus)}}

Next re-load the affinity analysis:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
data(tamoxifen_analysis)
@

To compare the sites unique to each sample group identified from the occupancy analysis with those sites identified as differentially bound based on affinity (read count) data, we use a feature of \Rcode{dba.report} that facilitates evaluating the occupancy status of sites. Here we obtain a report of all the sites (\Rcode{th=1}) with occupancy statistics (\Rcode{bCalled=T}):

<<eval=TRUE,echo=TRUE,results=hide>>=
tamoxifen.rep = dba.report(tamoxifen,bCalled=T,th=1)
@

The \Rcode{bCalled} option adds two columns to the report (\Rcode{Called1} and \Rcode{Called2}), one for each group, giving the number of samples within the group in which the site was identified as a peak in the original peaksets generated by the peak caller. We can use these to recreate the overlap criteria used in the occupancy analysis:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
onlyResistant = tamoxifen.rep$Called1>=2 & tamoxifen.rep$Called2<3
sum(onlyResistant )
onlyResponsive = tamoxifen.rep$Called2>=3 &  tamoxifen.rep$Called1<2
sum(onlyResponsive)
bothGroups = tamoxifen.rep$Called1>= 2 & tamoxifen.rep$Called2>=3
sum(bothGroups)
@

Comparing these numbers verifies the similarity with those seen in Figure~\ref{DiffBind-tamox_cons_venn}, showing again how the basic analysis is not oversensitive to differences in how the consensus peaksets are formed. 
This overlap analysis suggests that
\Sexpr{sum(onlyResistant)+sum(onlyResponsive)}
of the sites are uniquely bound in either the Responsive or Resistant groups, while
\Sexpr{sum(bothGroups)}
sites are common to both. 

Completing a full differential analysis and focusing on only those sites identified as significantly differentially bound (FDR <= 0.1), however, shows a different story than that obtainable using only occupancy data:

<<eval=TRUE,echo=TRUE,results=verbatim>>=
tamoxifen.DB = dba.report(tamoxifen,bCalled=T,th=.1)
onlyResistant.DB = tamoxifen.DB$Called1>=2 & tamoxifen.DB$Called2<3
sum(onlyResistant.DB)
onlyResponsive.DB = tamoxifen.DB$Called2>=3 & tamoxifen.DB$Called1<2
sum(onlyResponsive.DB)
bothGroups.DB = tamoxifen.DB$Called1>=2 & tamoxifen.DB$Called2>=3
sum(bothGroups.DB)
@

There are a number of notable differences in the results. First, overall there are many fewer sites identified as differentially bound 
(\Sexpr{sum(onlyResistant.DB)}+\Sexpr{sum(onlyResponsive.DB)}+\Sexpr{sum(bothGroups.DB)} = 
\Sexpr{length(tamoxifen.DB$Called1)})
than are unique to one condition
(\Sexpr{sum(onlyResistant)}+\Sexpr{sum(onlyResponsive)} =
\Sexpr{sum(onlyResistant)+sum(onlyResponsive)}).
Indeed, most of the sites identified in the occupancy analysis as unique to a sample group are not found to be significantly differentially bound using the affinity data. 
While partly this is a result of the stringency of the statistical tests, it shows how the affinity analysis can discriminate between sites where peak callers are making occupancy decisions that do not reflect significant differences in read densities at these sites. 
Note that only about 
\Sexpr{sum(sum(onlyResistant.DB), sum(onlyResponsive.DB))/sum(sum(onlyResistant),sum(onlyResponsive))*100}\%
of sites unique to one condition are identifiable as significantly differentially bound 
(\Sexpr{sum(onlyResistant.DB)}+\Sexpr{sum(onlyResponsive.DB)} =
\Sexpr{sum(onlyResistant.DB) + sum(onlyResponsive.DB)} 
out of 
\Sexpr{sum(onlyResistant)+sum(onlyResponsive)}).
Secondly, differentially bound sites are as likely to be called in the consensus of both response groups as they are to be unique to one group, as over a quarter of the total sites identified as significantly differentially bound (\Sexpr{sum(sum(onlyResistant.DB),sum(onlyResponsive.DB),sum(bothGroups.DB))}) are called as peaks in \emph{both} response groups (\Sexpr{sum(bothGroups.DB)}). 
Indeed, peaks identified binding sites in both sample groups are almost as likely to be identified as significantly differentially bound  
(\Sexpr{signif(sum(bothGroups.DB)/sum(bothGroups)*100,2)}\%, 
or 
\Sexpr{sum(bothGroups.DB)} out of 
\Sexpr{sum(bothGroups)}) 
as sites identified only in the Resistant sample group 
(\Sexpr{signif((sum(onlyResistant.DB)/sum(onlyResistant))*100,2)}\%, 
or \Sexpr{sum(onlyResistant.DB)} out of \Sexpr{sum(onlyResistant)}).

A final advantage of a quantitative analysis is that the differentially bound peaks identified using the affinity analysis are associated with significance statistics (p-value and FDR) that can be used to rank them for further examination, while the occupancy analysis yields a relatively unordered list of peaks, as the peak caller statistics refer only to the significance of occupancy, and not of differential binding.

\section{Technical notes}
This section includes some technical notes explaining some of the  technical details of \DBA processing.

\subsection{Loading peaksets}
There are a number of ways to get peaksets loaded into a DBA object. Peaksets can be read in from files or loaded from interval sets already stored in an R object. Samples can be specified either in a sample sheet  (using \Rfunction{dba}) or loaded one at a time (using \Rfunction{dba.peakset}).

When loading in peaksets from files, specifying what peak caller generated the file enables peaks from supported peak callers to be read in. See the help page for \Rfunction{dba.peakset} for a list of supported peak callers. Any string to indicate the peak caller; if it is not one of the supported callers, a default "raw" format is assumed, consisting of a text file with three or four columns (indicating the chromosome, start position, and end position, with a score for each interval found in the fourth column, if present). You can further control how peaks are read using the \Rcode{PeakFormat}, \Rcode{ScoreCol}, and \Rcode{bLowerBetter} fields if you want to override the defaults for the specified peak caller identifier. For example, with the tamoxifen dataset used in this tutorial, the peaks were called using the MACS peak caller, but the data are supplied as simple text files, no the expected MACS "xls" format. To maintain the peak caller in the metadata, we could specify the \Rcode{PeakCaller} as "macs" but the \Rcode{PeakFormat} as "raw". If  we wanted to use peak scores in a different column than the fourth, the \Rcode{scorecol} parameter could be set to indicate the appropriate column number.  When handling scoring, \DBA by default assumes that a higher score indicates a "better" peak. If this is not the case, for example if the score is a p-value or FDR, we could set  \Rcode{bLowerScoreBetter} to \Rcode{TRUE}.

When using a sample sheet,  values for fields missing in the sample sheet can be supplied when calling \Rfunction{dba}. In addition to the minimal sample sheet used for the tutorial, an equivalent sample sheet with all the metadata fields is included, called "tamoxifen\_allfields.csv". See the help page for \Rfunction{dba} for an example using this sample sheet.

\subsection{Merging peaks}
When forming the global binding matrix consensus peaksets, \DBA first identifies all unique peaks amongst the relevant peaksets. As part of this process, it merges overlapping peaks, replacing them with a single peak representing the narrowest region that covers all peaks that overlap by at least one base. There are at least two consequences of this that are worth noting.

First, as more peaksets are included in analysis, the average peak width tends to become longer as more overlapping peaks are detected and the start/end points are adjusted outward to account for them. Secondly, peak counts may not appear to add up as you may expect due to merging. For example, if one peakset contains two small peaks near to each other, while a second peakset includes a single peak that overlaps both of these by at least one base, these will all be replaced in the merged matrix with a single peak. As more peaksets are added, multiple peaks from multiple peaksets may be merged together to form a single, wider peak.

\subsection{\edgeR analysis}
When \Rfunction{dba.analyze} is invoked using the default \Rcode{method=DBA\_EDGER}, a standardized differential analysis is performed using the \edgeR package (\cite{Robinson:2010p249}). This section details the precise steps in that analysis.

For each contrast, a separate analysis is performed. 
First, a matrix of counts is constructed for the contrast, with columns for all the samples in the first group, followed by columns for all the samples in the second group. 
The raw read count is used for this matrix; if the \Rcode{bSubControl} parameter is set to \Rcode{TRUE} (as it is by default), the raw number of reads in the control sample (if available) will be subtracted (with a minimum final read count of 1). 
Next the library size is computed for each sample for use in subsequent normalization. 
By default, this is the total number of reads in the library (calculated from the source BAM//BED file).
Alternatively, if the \Rcode{bFullLibrarySize} parameter is set to FALSE,the total number of reads in peaks (the sum of each column) is used. 
Note that \"effective" library size (\Rcode{bFullLibrarySize = FALSE}) may be more appropriate for situations when the overall signal (binding rate) is expected to be directly comparable between the samples.
Next comes a call to \edgeR's \Rfunction{DGEList} function. The \Rclass{DGEList} object that results is next passed to \Rcode{calcNormFactors} with all other parameters retained as defaults  (\Rcode{method="TMM"}), returning an updated \Rcode{DGEList} object. 
This is passed to \Rcode{estimateCommonDisp} with default parameters. 

If the method is \Rcode{DBA\_EDGER\_CLASSIC}, then if \Rcode{ bTagwise} is TRUE (most useful when there are at least three members in each group of a contrast), the resulting \Rclass{DGEList} object is then passed to \Rcode{estimateTagwiseDisp}, with the prior set to 50 divided by two less than the total number of samples in the contrast, and \Rcode{trend="none"}. 
The final steps are to perform testing to determine the significance measure of the differences between the sample groups by calling \Rfunction{exactTest} (\cite{RobinsonandSmyth}) using the \Rclass{DGEList} with the \Rcode{dispersion} set based on the \Rcode{bTagwise} parameter.

If the method is \Rcode{DBA\_EDGER\_GLM} (the default), then a a design matrix is generated with two coefficients (the Intercept and one of the groups).
Next  \Rfunction{estimateGLMCommonDisp} is called; if \Rcode{bTagwise=TRUE},  \Rfunction{estimateGLMTagwiseDisp} is called as well. 
The model is fitted by calling \Rfunction{glmFit}, and the specific contrast fitted by calling \Rfunction{glmLRT}, specifying that the second coefficient be dropped.
Finally, an \Rfunction{exactTest} (\cite{McCarthyetal}) is performed, using either common or tagwise dispersion depending on the value specified for \Rcode{bTagwise}.

This final  \Rclass{DGEList} for contrast n is stored in the \Rclass{DBA} object as 

\Rcode{DBA\$contrasts[[n]]\$edgeR} 

and may be examined and manipulated directly for further customization. Note however that if you wish to use this object directly with \edgeR functions, then the \Rcode{bReduceObjects} parameter should be set to FALSE, otherwise the default value of TRUE will result in essential object fields being stripped.

If a blocking factor has been added to the contrast, an additional \edgeR analysis is carried out. 
This follows the \Rcode{DBA\_EDGER\_GLM} case detailed above, except a more complex design matrix is generated that includes all the unique values for the blocking factor. These coefficients are all included in the \Rfunction{glmLRT} call.  The resultant object is accessible as

\Rcode{DBA\$contrasts[[n]]\$edgeR\$block}. 


\subsection{\DESeq analysis}
When \Rfunction{dba.analyze} is invoked using \Rcode{method=DBA\_DESEQ}\footnote{Note that \DESeq can be made the default analysis method for a DBA object by setting \Rcode{DBA\$config\$AnalysisMethod=DBA\_DESEQ}.}, a standardized differential analysis is performed using the \DESeq package (\cite{Anders:2010p792}). This section details the  steps in that analysis.

For each contrast, a separate analysis is performed. 
First, a matrix of counts is constructed for the contrast, with columns for all the samples in the first group, followed by columns for all the samples in the second group. 
The raw read count is used for this matrix; if the \Rcode{bSubControl} parameter is set to \Rcode{TRUE} (as it is by default), the raw number of reads in the control sample (if available) will be subtracted. 
Next the library size is computed for each sample for use in subsequent normalization. 
By default, this is the total number of reads in the library (calculated from the source BAM//BED file).
Alternatively, if the \Rcode{bFullLibrarySize} parameter is set to FALSE,the total number of reads in peaks (the sum of each column) is used. 
Note that "effective" library size (\Rcode{bFullLibrarySize = FALSE}) may be more appropriate for situations when the overall signal (binding rate) is expected to be directly comparable between the samples.
The first step concludes with a call to \Biocpkg{DESeq}'s \Rfunction{newCountDataSet} function, which returns a \Rclass{CountDataSet} object. 
If \Rcode{bFullLibrarySize} is set to TRUE, then \Rfunction{sizeFactors} is called with the number of reads in the BAM/BED files for each ChIP sample, divided by the minimum of these; otherwise, \Rfunction{estimateSizeFactors} is invoked. 
Next, \Rfunction{estimateDispersions} is called with the \Rclass{CountDataSet} object and \Rcode{fitType} set to \Rcode{local}. If there are no replicates, (only one sample in each group), \Rcode{method} is set to \Rcode{blind}. 
Otherwise, if \Rcode{bTagwise} is TRUE, \Rcode{method} is set to \Rcode{per-condition}; if it is FALSE, \Rcode{method} is set to \Rcode{pooled} (or \Rcode{pooled-CR} for a blocking analysis). 

If the method is \Rcode{DBA\_DESEQ\_CLASSIC}, \Rfunction{nbinomTest} is called, and the result (reordered by adjusted p-value) saved for reporting.

If the method is \Rcode{DBA\_DESEQ\_GLM} (the default), two models are fitted using \Rfunction{fitNbinomGLMs}: a full model is fitted with all the coefficients, and a second model is fitted with the second coefficient dropped. These are tested against each other using \Rfunction{nbinomGLMTest}, with the resulting p values adjusted using \Rfunction{p.adjust} (with \Rcode{method="BH"}).

The final results are accessible within the \Rclass{DBA} object as 

\Rcode{DBA\$contrasts[[n]]\$DESeq1\$DEdata} 

and may be examined and manipulated directly for further customization. Note however that if you wish to use this object directly with \DESeq functions, then the \Rcode{bReduceObjects} parameter should be set to FALSE, otherwise the default value of TRUE will result in essential object fields being stripped.

If a blocking factor has been added to the contrast, an additional \DESeq analysis is carried out. 
This follows the \Rcode{DBA\_DESEQ\_GLM} case detailed above, except a more complex design is generated when \Rfunction{newCountDataSet} is called that includes all the unique values for the blocking factor. These coefficients are all included in the \Rfunction{fitNbinomGLMs} calls.  The resultant object is accessible as

\Rcode{DBA\$contrasts[[n]]\$DESeq1\$block\$DEdata}. 

\subsection{\Biocpkg{DESeq2} analysis}
When \Rfunction{dba.analyze} is invoked using \Rcode{method=DBA\_DESEQ2}\footnote{Note that \Biocpkg{DESeq2} can be made the default analysis method for a DBA object by setting \Rcode{DBA\$config\$AnalysisMethod=DBA\_DESEQ2}.}, a standardized differential analysis is performed using the \Biocpkg{DESeq2}package (\cite{Love2014}). This section details the precise steps in that analysis.

For each contrast, a separate analysis is performed. 
First, a matrix of counts is constructed for the contrast, with columns for all the samples in the first group, followed by columns for all the samples in the second group. 
The raw read count is used for this matrix; if the \Rcode{bSubControl} parameter is set to \Rcode{TRUE} (as it is by default), the raw number of reads in the control sample (if available) will be subtracted.
Next the library size is computed for each sample for use in subsequent normalization. 
By default, this is the total number of reads in peaks (the sum of each column). 
Alternatively, if the \Rcode{bFullLibrarySize} parameter is set to TRUE, the total number of reads in the library (calculated from the source BAM/BED file) is used. 
The first step concludes with a call to \Biocpkg{DESeq2}'s \Rfunction{DESeqDataSetFromMatrix} function, which returns a \Rclass{DESeqDataSet} object. 

If \Rcode{bFullLibrarySize} is set to TRUE, then \Rfunction{sizeFactors} is called with the number of reads in the BAM/BED files for each ChIP sample, divided by the minimum of these; otherwise, \Rfunction{estimateSizeFactors} is invoked. 

\Rfunction{estimateDispersions} is then called with the \Rclass{DESeqDataSet} object and \Rcode{fitType} set to \Rcode{local}. Next the model is fitted and tested using \Rfunction{nbinomWaldTest}.
The final results (as a \Rclass{DESeqDataSet}) are accessible within the \Rclass{DBA} object as 

\Rcode{DBA\$contrasts[[n]]\$DESeq2\$DEdata} 

and may be examined and manipulated directly for further customization. Note however that if you wish to use this object directly with \Biocpkg{DESeq2} functions, then the \Rcode{bReduceObjects} parameter should be set to FALSE, otherwise the default value of TRUE will result in essential object fields being stripped.

If a blocking factor has been added to the contrast, an additional \Biocpkg{DESeq2} analysis is carried out by setting the \Rcode{design} to include all the unique values for the blocking factor. This occurs before the dispersion values are calculated. The resultant \Rclass{DESeqDataSet} object is accessible as

\Rcode{DBA\$contrasts[[n]]\$DESeq2\$block\$DEdata}. 

\section{Vignette Data}
Due to space limitations, the aligned reads associated with the cell line data used in this vignette are not included as part of the \DiffBind package. 

If you would like to get access to the reads, you can email \Rcode{rory.stark@cruk.cam.ac.uk} and ask to be added to the Dropbox containing the vignette data.

The data are also available in the Short Read Archive (SRA), email for detailed instructions on how to retrieve them in the appropriate form.

Also included is a sample sheet that uses the downloaded data. This sample sheet can also be examined to see which GEO accession IDs correspond to which samples. The sample sheet can be found at:

<<eval=FALSE,echo=TRUE,results=hide>>=
file.path(system.file("extra", package="DiffBind"),"tamoxifen_GEO.csv")
@

\section{Using \Biocpkg{DiffBind} and \Biocpkg{ChIPQC} together}

\Biocpkg{DiffBind} and \Biocpkg{ChIPQC} are both packages that help manage and analyze ChIP-seq experiments, and are designed to be used together.

If you already have a project in \Biocpkg{DiffBind}, then \Biocpkg{ChIPQC} can accept a \Rclass{DBA} object in place of the sample sheet when creating a \Rclass{ChIPQCexperiment} object.

Once a \Rclass{ChIPQCexperiment} object has been constructed, it can be used in place of a \Rclass{DBA} object in most calls to \Biocpkg{DiffBind}.
All plotting, counting, and analysis functions are available from \Biocpkg{DiffBind}. 

It is also possible to extract a \Rclass{DBA} object from a \Rclass{ChIPQCexperiment} object using the \Rfunction{QCdba} method.
The resulting \Rclass{DBA} object can be used in \Biocpkg{DiffBind} without restriction, although neither it nor \Rclass{DBA} objects based on it can be re-attached to the original \Rclass{ChIPQCexperiment} object (although they can be used in lieu of a sample sheet when creating a new one.) 

In a typical workflow, the first step would be to run a \Biocpkg{ChIPQC} analysis before peak calling to assess library quality and establish what filtering should be done at the read level (mapping quality, duplicates, and blacklists).
Next peaks would be called externally, and read into a new \Rclass{ChIPQCexperiment} object to assess peak-based metrics, such as FRIP, peak profiles, and clustering.

At this point, \Biocpkg{DiffBind} could be used to perform occupancy analysis, derive consensus peak sets, re-count reads to form a binding matrix, and set up contrasts to carry out full differential binding analyses using the \Biocpkg{edgeR} and \Biocpkg{DESeq2} packages, along with  plotting and reporting functions.

\section{Acknowledgements}

This package was developed at Cancer Research UK's Cambridge Research Institute with the help and support of many people there. We wish to acknowledge everyone the Bioinformatics Core under the leadership of Matthew Eldridge, as well as the Nuclear Receptor Transcription Laboratory under the leadership of Jason Carroll. Researchers who contributed ideas and/or pushed us in the right direction include Caryn-Ross Innes, Vasiliki Theodorou, and Tamir Chandra among many others. We also thank  members of the Gordon Smyth laboratory at the WEHI, Melbourne, particularly Mark Robinson and Davis McCarthy, for helpful discussions.


\section{Session Info}
<<sessionInfo, results=tex, print=TRUE, eval=TRUE>>=
toLatex(sessionInfo())
@

<<eval=TRUE,echo=FALSE,results=hide>>=
setwd(savewd)
@
\bibliography{DiffBind}

\end{document}