%\VignetteIndexEntry{An Introduction to the MMDiff2 method}
%\VignetteKeywords{MMDiff, MMDiff2, ChIP-Seq, kernel methods}
%\VignettePackage{MMDiff2}
%\VignetteEngine{knitr::knitr}

% To compile this document
% library('knitr'); rm(list
% !Rnw weave = knitr
\newcommand{\DBA}{\texttt{DiffBind}}
\newcommand{\DESeq}{\texttt{DESeq}}
\newcommand{\edgeR}{\texttt{edgeR}}
\newcommand{\MD}{\texttt{MMDiff2}}
\newcommand{\code}[1]{{\small\texttt{#1}}}
%\newcommand{\R}{\textsf{R}}

\newcommand{\tab}[1]{Table \ref{tab:#1}}
\newcommand{\fig}[1]{Figure \ref{fig:#1}}

\documentclass{article}

<<style-Sweave, eval=TRUE, echo=FALSE, results='asis'>>=
BiocStyle::latex()
@

\title{MMDiff2: Statistical Testing for ChIP-Seq Data Sets}
\author{Gabriele Schwikert $\mathsf{<}$\texttt{G.Schweikert@ed.ac.uk}$\mathsf{>}$  $^{1}$ and David Kuo $\mathsf{<}$\texttt{dkuo@cbio.mskcc.org}$\mathsf{>}$ $^{2}$\\[1em]
  \small{$^{1}$ The Informatics Forum, University of Edinburgh and} \\
  \small{The Wellcome Trust Center for Cell Biology, University of Edinburgh;} \\
  \small{$^{2}$ Department of Physiology, Biophysics and Systems Biology, Weill Cornell Medical College and} \\
  \small{Computational Biology Department, Memorial Sloan Kettering Cancer Center}}

\date{Modified: 25 Mar, 2016. Compiled: \today}

\begin{document}

\maketitle

\tableofcontents

\section{Introduction}

ChIP-seq is the standard technique for determining genome-wide
protein-DNA interactions and of particular importance for the fields
of regulatory biology and epigenetics. A common use-case for ChIP-seq is
determining where transcription factors bind on the genome.
The resulting ChIP-seq peaks for transcription factors are
typically \textit{narrow} and can be further queried for known DNA-binding motifs.
Another example is ascertaining chromatin states by profiling histone
modifications. Specific post-translational modification on histone proteins,
represented as \textit{broad} peaks, are highly informative of the gene
activity and regulatory status.

The utility and reproducibility of ChIP-seq led to its use in large-scale
projects such as \textit{ENCODE} and \textit{Epigenomics Roadmap}.
Thousands of ChIP-seq datasets were generated for various families
of transcription factors and histone modifications in cell lines and
primary tissues. These projects described the regulatory
landscape through the lens of ChIP-seq experiments and also established
a standard for comparison with other experimental datasets.

With an abundance of ChIP-seq data, there have also emerged many different
computational tools for the analysis and statistical testing of ChIP-seq
data between conditions, i\.e\., control cells against experimentally treated
cells. After mapping reads to the genome, ChIP-seq analysis begins by calling
peaks between a sample conjugated with an antibody against a
control with DNA input. Various computational tools exist for this
peak calling task but a popular tool is \textit{MACS2}. In another scenario,
regions of interest are pre-defined with regard to other annotated objects,
for example by applying windows around transcription start sites (TSS).
In querying all of the TSS for possible ChIP-seq reads, experimentalists can
observe the density of reads near regions of interest.

After regions/peaks are identified, a common analysis is determining regions
which are differentially bound by the studied transcription factor, or which
peaks demonstrate differential histone modifications between conditions.
For simplicity, we will call these regions differential regions. Often this is
treated as a read-count based test, and some statistical packages
such as \DESeq2 can be utilized for this task. ChIP-seq specific tools
such as \DBA also integrate count-based testing in their methods.

A caveat of these methods is that count-based methods disregard higher order
features of ChIP-seq reads such as peak shape. Between two conditions,
differential read density may not be significant but the distribution of reads
within the peak may imply differential protein binding or the presence of
protein complexes. Detection of these differential peak shapes is not possible
when simplifying peaks to single values, as is the case for the aforementioned
read-count based tests.

We therefore incorporate higher order ChIP-seq peak information into
a test for differential binding by adapting kernel-based
statistical tests to ChIP-Seq data. This package was initially released as
\textit{MMDiff}.

In this current release, \textit{MMDiff2}, we improve code stability, decrease analysis runtime and provide an interactive Shiny application for exploring peak profiles.

%ChIP-Seq has rapidly become the dominant experimental technique to
%determine the location of transcription factor binding sites and
%histone modifications. Typically, computational peak finders, such as
%MACS \citep{MACS}, are used to identify potential candidate regions,
%i.e. regions with significantly enriched read coverage relative to
%some background. In the following, we call these regions
%{\it peaks} and assume that their {\it genomic coordinates} are
%provided. Going beyond this basic analysis, it is often of interest to
%detect a subset of peaks where significant {\it changes of read
%coverage} occur in a treatment experiment relative to a control (see
%\fig{Example}). Despite the abundance of sequencing data,
%statistical analyses of ChIP-Seq data remains
%challenging, due to the highly structured nature of the data and the
%paucity of replicates.  Current approaches to detect differentially
%bound regions are mainly borrowed from RNA-Seq data analysis, thus
%focusing on total counts of fragments mapped to a region and
%ignoring any information encoded in the shape of the peak profile.
%
%Higher order features of ChIP-Seq peak enrichment profiles carry
%important and often complementary information to total counts, and
%hence are potentially important in assessing differential binding. For example,
%peaks of DNA-binding proteins that complex with other factors may have differentially
%shapped peaks compared to peaks without co-factor binding. Identification of
%differentially shaped peaks requires methods beyond count-based methods commonly
%associated with RNA-seq differential expression.

\section{Quick Start}
<<quickstart,message=FALSE,eval=FALSE>>=

library('MMDiff2')
library('MMDiffBamSubset')
ExperimentData <- list(genome='BSgenome.Mmusculus.UCSC.mm9',
                       dataDir=system.file("extdata", package="MMDiffBamSubset"),
                       sampleSheet="Cfp1.csv")
MetaData <- list('ExpData' = ExperimentData)
MMD <- DBAmmd(MetaData)
data("Cfp1-Peaks")
MMD <- setRegions(MMD,Peaks)
MMD <- getPeakReads(MMD)
MMD <- estimateFragmentCenters(MMD)
MMD <- compHists(MMD)
MMD <- compDists(MMD)
MMD <- setContrast(MMD,contrast='byCondition')
MMD <- compPvals(MMD)
res <- reportResults(MMD)
@

\section{Analysis Pipeline}

\MD requires position sorted and indexed \texttt{BAM} files of the sample and
input libraries and a set of regions in \texttt{GRanges} format, for example
\texttt{MACS2} called peaks. For this vignette, we utilize the
\texttt{MMDiffBamSubset} dataset that contains ChIP-seq peaks and reads from Clouaire \textit{et al.}, Genes \& Dev. 2012.

<<setup,message=FALSE>>=

# load software package
library('MMDiff2')
@

\subsection{Loading Data}
We load the \texttt{MMDiffBamSubset} data and specify an experimental genome
(mm9), a directory containing peaks and reads, and a \textit{csv} file that details
the paths to different samples. The \texttt{MetaData} object will be called
upon later for adding binding motif and gene annotation information.

<<dataload,message=FALSE>>=
# load data packages
library('MMDiffBamSubset')

# create metaData:

ExperimentData <- list(genome = 'BSgenome.Mmusculus.UCSC.mm9',
                    dataDir = system.file("extdata", package="MMDiffBamSubset"),
                    sampleSheet="Cfp1.csv")

MetaData <- list('ExpData' = ExperimentData)

@

\subsection{Processing Peaks from BAM Files}
A \texttt{GenomicRanges} object that summarizes the peaks has been provided
in this vignette. After loading this file, we instantiate the \texttt{DBAmmd}
class and read the \texttt{BAM} files for all peaks in all conditions.

<<regions,message=FALSE>>=

data('Cfp1-Peaks')
MMD <- DBAmmd(MetaData)
MMD <- setRegions(MMD,Peaks)
MMD <- getPeakReads(MMD)
@

Note that we are fetching the exact start and end positions of mapped fragments
and not the coverage. In the case of single-end reads, the left-most positions
of fragments mapping to the positive strand are stored in a list called
\texttt{Left.p}, and the right most positions of fragments mapping to the
negative strand are stored in \texttt{Right.n}. For paired-end reads,
\texttt{Right.p} and \texttt{Left.n} are additionally kept.

The next step in the pipeline is to estimate the fragment centers using
information from the forward and reverse strands. For paired-end reads, we
compute the exact Center position for each fragment. For single-strand
libraries, we estimate the Centers. We improve the resolution of the estimated
Center positions, relative to coverage-based approach, by utilizing the method
described in Gomes \textit{et al.} Genome Research, 2014.
%TODO: properly add refernce

<<fragC, message=FALSE, markup='hide'>>=
MMD <- estimateFragmentCenters(MMD)
@

Once the fragment centers have been estimated, histograms for each peak
are computed using the read information. The step is essential for plotting
regions and not essential in order to determine differentially shaped regions.
Two parameters as used to compute histograms:
\texttt{bin.length}, which determines the smoothing parameter over the
set of reads (Default: 20) and the \\{whichPos}, which determines the fragment
position from which histograms should be computed: \texttt{'Left.p'} for
left ends of fragments mapping to positive strand, \texttt{'Right.n'} for
right ends of fragments mapping to negative strands and  \texttt{'Center'} for
the Fragment centers for all positions (Default: \texttt{Center}).


<<hists,message=FALSE>>=
MMD <- compHists(MMD, bin.length=20, whichPos="Center")
@

\subsection{Examining the \texttt{DBAmmd} Object and Reporting}

% TODO : does this section make sense here: explain how to get slots:
% Regions(MMD), Samples(MMD), Hists(MMD), Counts(MMD) (see DBA-methods)

After peak histograms have been computed, we have a complete \texttt{DBAmmd}
object. This object contains several components including:
\begin{enumerate}
\item Peaks
\item Samples
\item Distances between Conditions (not yet computed)
\item Distances within Conditions (not yet computed)
\item P-values (not yet computed)
\end{enumerate}

Each of these components will be highlighted in the subsequent sections.

\subsection{Plotting Peaks}
Peaks can be visualized by providing a \texttt{Peak.id} parameter to the
\texttt{plotPeak} function. Line plots are drawn that signify the different
histograms from each of the samples.

<<peakplot1, markup='hide', message=FALSE, fig.width=7, fig.height=3, fig.align='center' >>=
plotPeak(MMD, Peak.id='241', plot.input = FALSE, whichPos="Center")
@

\subsubsection{Motif Enrichment within Peaks}
An additional feature to \MD is plotting motifs within peaks. We specify
motifs with the following command against the \texttt{MotifDb}:

<<motif specifcation>>=
library('MotifDb')
motifs <- query(query(MotifDb, 'Mmusculus'), 'E2F')
@

The query will return motifs that partially match the term "Pax", i\.e\., Pax4
and Pax5 could be plotted. Integration of the motifs is seamless into the
\texttt{plotPeak} function:

<<peakplot2,message=FALSE,fig.width=7, fig.height=4, fig.align='center'>>=
plotPeak(MMD, Peak.id='241', NormMethod=NULL,plot.input = FALSE,whichPos="Center",
  Motifs=motifs,Motifcutoff="80%")
@
\subsubsection{Gene Annotation}

In addition to plotting motifs of interest, we have enabled in \MD gene
annotation tracks that can also be plotted in a fully integrated plot
with the ChIP-seq peak by providing a \texttt{genomicRanges} object to
the \texttt{anno} parameter.

<<peakplot3,message=FALSE,fig.width=7, fig.height=5, fig.align='center'>>=
data("mm9-Genes")
names(GR) <- GR$tx_name
GR <- list(UCSCKnownGenes = GR)
plotPeak(MMD, Peak.id='241', NormMethod=NULL,plot.input = FALSE,
         whichPos="Center",Motifs=motifs, anno=GR)
@

\subsection{Computing per-region pair-wise distances between samples}
The core of \MD is the comparison of peak shapes between samples. This is
performed by using a kernel-based test that calculates 1) a distance between
treatment conditions and 2) within conditions. The \texttt{compDists} function
generates these per-region and pair-wise distances for all samples in the
experiment.

<<dists, message=FALSE, include=FALSE>>=
MMD <- compDists(MMD, dist.method = "MMD", run.parallel = FALSE)
@

\subsection{Differential Testing}
After distances are generated, \MD computes p-values based on the specified
conditions.

<<contrast,>>=
MMD <- setContrast(MMD,contrast='byCondition')
@

The contrast can also be manually set:
<<contrast2,message=FALSE, include=FALSE>>=
group1 <- Samples(MMD)$Condition=='1'
names(group1) <- Samples(MMD)$SampleID
group2 <- Samples(MMD)$Condition=='2'
names(group2) <-  Samples(MMD)$SampleID

contrast <- list(group1=group1,
                 group2=group2,
                 name1='WT,REsc',
                 name2='Null')

#setContrast(MMD,contrast=contrast)
@


<<diffTest,>>=
MMD <- compPvals(MMD,dist.method='MMD')
@

The \texttt{plotDists} function plots a summary of all of the peak distances
1) between groups and 2) within groups. Significant differentially shaped peaks
are highlighted in red.

\subsection{Analyzing Results}
<<>>=
plotDists(MMD, dist.method='MMD',whichContrast=1,
                          diff.method='MMD.locfit',
                          bUsePval=FALSE, th=0.1,
                          title=NULL, what=3,
                          xlim=NULL,ylim=NULL,Peak.IDs=NULL,
                          withLegend=TRUE)

@

% TODO chose a better peak, with higher coverage maybe

<<>>=

res <- reportResults(MMD)

Peak.ids <- names(res)

plotPeak(MMD, Peak.id=Peak.ids[1], NormMethod=NULL,plot.input = FALSE,
         whichPos="Center",Motifs=motifs, anno=GR,whichContrast = 1)
dev.off()
plotDISTS4Peak(MMD,Peak.id=Peak.ids[1],dist.method='MMD',
               whichContrast=1,Zoom=TRUE)

@

\subsection{Shiny Application Usage}

\MD~features an interactive Shiny application that displays
distances plots showing all peaks and in-depth peak plots.

Beginning in the top panels of the Shiny application are two
\texttt{plotDists} outputs. In the left panel, drawing a box dynamically zooms in
on the right panel. The drawn box can also be panned around the plot. To remove
the zoomed view, click once outside the box.

On the right panel, clicking anywhere on the plot will find the nearest peak.
The nearest peak to the mouse-click will then generate two plots below.

The bottom left panel is a histogram representation of all conditions
(\texttt{plotPeak}) while the
bottom right panel displays the distances between samples for that peak
(\texttt{plotDISTS4peak}).

Finally, a UCSC genome browser link is generated that matches the coordinates
of the selected peak for additional inspection.

<<shinyApp,eval=FALSE>>=
runShinyMMDiff2(MMD)
@

A demonstration of this is shown below.

\includegraphics[width=\linewidth]{ShinyMMDiff2Demo.png}

\section{Comparison with DiffBind}

To compare against the \DBA~package, we downloaded aligned \texttt{BAM} files
from the Ross-Innes \textit{et al.} 2012 paper and used the same called peaks
as analyzed in the \Biocpkg{DiffBind}~vignette.

\section{Setup}

This vignette was built using:
<<>>=
sessionInfo()
@
\end{document}