% -*- mode: noweb; noweb-default-code-mode: R-mode; -*-
\documentclass[11pt]{article}
%% Set my margins
\setlength{\oddsidemargin}{0.0truein}
\setlength{\evensidemargin}{0.0truein}
\setlength{\textwidth}{6.5truein}
\setlength{\topmargin}{0.0truein}
\setlength{\textheight}{9.0truein}
\setlength{\headsep}{0.0truein}
\setlength{\headheight}{0.0truein}
\setlength{\topskip}{0pt}
%% End of margins

%%\pagestyle{myheadings}
%%\markboth{$Date$\hfil$Revision$}{\thepage}
\usepackage[pdftex,
bookmarks, 
bookmarksopen,
pdfauthor={Dongjun Chung, Pei Fen Kuan and Sunduz Keles},
pdftitle={mosaics Vignette}]
{hyperref}

\title{Analysis of ChIP-seq Data with `\texttt{mosaics}' Package}
\author{Dongjun Chung$^1$, Pei Fen Kuan$^2$ and S\"und\"uz Kele\c{s}$^{1,3}$\\
  $^1$Department of Statistics, University of Wisconsin\\ 
  Madison, WI 53706.\\
  $^2$Department of Biostatistics, University of North Carolina at Chapel Hill\\
  Chapel Hill, NC 27599.\\
  $^3$Department of Biostatistics and Medical Informatics, 
  University of Wisconsin\\Madison, WI 53706.}

\date{\today}

\SweaveOpts{engine=R, echo=TRUE, pdf=TRUE}

\begin{document}
%\VignetteIndexEntry{MOSAiCS}
%\VignetteKeywords{MOSAiCS}
%\VignettePackage{mosaics}
\maketitle

\section{Overview}

This vignette provides an introduction to the analysis of ChIP-seq data with
the `\texttt{mosaics}' package. R package \texttt{mosaics} implements MOSAiCS,
a statistical framework for the analysis of ChIP-seq data, proposed 
in \cite{mosaics}. MOSAiCS stands for ``\textbf{MO}del-based one and
two \textbf{S}ample \textbf{A}nalysis and \textbf{I}nference for
\textbf{C}hIP-\textbf{S}eq Data''. It implements a flexible parametric mixture
modeling approach for detecting peaks, i.e., enriched regions, in one-sample
(ChIP sample) or two-sample (ChIP and control samples) ChIP-seq data.
It accounts for mappability and GC content biases that arise in ChIP-seq data.

The package can be loaded with the command:

<<preliminaries,echo=FALSE,results=hide>>=
options(prompt = "R> ")
@

<<mosaics-prelim>>=
library("mosaics")
@

\section{Getting started}

We assume that you already have the aligned files for your samples
(ChIP and control)  such as files obtained from the ELAND aligner.
R package `\texttt{mosaics}' takes bin-level data as input and
these bin-level data can easily be generated from the aligned data 
using the preprocessing scripts we provide. You can download these 
preprocessing scripts from the `\texttt{mosaics}' package companion website,
\url{http://www.stat.wisc.edu/~keles/Software/mosaics/}. You also need to have
bin-level mappability, GC content, and sequence ambiguity 
score files for the reference genome you are working with. 

If you are working with organisms such as human (HG18 and HG19), mouse (MM9), 
rat (RN4), and Arabidopsis (TAIR9), you can download their corresponding 
preprocessed mappability, GC content, and sequence ambiguity score files 
at \url{http://www.stat.wisc.edu/~keles/Software/mosaics/}. 
If your reference genome of interest is not listed on our website, 
you can inquire about it at our Google group, 
\url{http://groups.google.com/group/mosaics_user_group}, and we would be happy 
to add your genome of interest to the list. The companion website also provides 
all the related scripts and easy-to-follow instructions to prepare these files. 
Please check \url{http://www.stat.wisc.edu/~keles/Software/mosaics/} 
for more details. We encourage  questions or requests regarding 
`\texttt{mosaics}' package to be posted on  our Google group
\url{http://groups.google.com/group/mosaics_user_group}.

\section{Workflow: Two-Sample Analysis}

\subsection{Reading  Bin-Level Data into the R Environment}

`\texttt{mosaics}' package assumes chromosome-wise analysis of ChIP-seq data. 
For the two-sample analysis, you need preprocessed bin-level ChIP data, 
control sample data, mappability score, GC content score, and sequence
ambiguity score. In this vignette, we use chromosome 21 data from a ChIP-seq
experiment of STAT1 binding in interferon-$\gamma$-stimulated HeLa S3 cells
\cite{stat1}. `\texttt{mosaicsExample}' package provides this example dataset.

<<mosaicsExample-prelim>>=
library(mosaicsExample)
@

Bin-level data can be imported to the R environment with the command:

<<ts-readbin>>=
exampleBinData <- readBins( type=c("chip","input","M","GC","N"),
    fileName=c( system.file( file.path("extdata","chip_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","input_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","M_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","GC_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","N_chr21.txt"), package="mosaicsExample") ) )
@

For the  `\texttt{type}' argument, \texttt{"chip"}, \texttt{"input"},
\texttt{"M"}, \texttt{"GC"}, and 	\texttt{"N"} indicate bin-level ChIP data, 
control sample data, mappability score, GC content score, and sequence
ambiguity score, respectively. You need to specify the corresponding file names 
in `\texttt{fileName}'. `\texttt{mosaics}' package assumes that each file name 
in `\texttt{fileName}' is provided in the same order as in `\texttt{type}'. 

R package \texttt{mosaics} provides functions for generating simple summaries
of the data. The following command prints out basic information about the
bin-level data, such as number of bins and total ``effective tag counts''.
``Total effective tag counts'' is defined as the sum of the tag counts of all
bins. This value is usually larger than the sequencing depth since tags are
counted after extension to average fragment length and an extended fragment can
contribute to multiple bins. 

<<ts-bindata-show>>=
exampleBinData
@

`\texttt{print}' method returns the bin-level data in data frame format.

<<ts-bindata-print>>=
print(exampleBinData)[51680:51690,]
@

`\texttt{plot}' method provides exploratory plots for the ChIP data. Different
type of plots can be obtained by varying the `\texttt{plotType}' argument.
`\texttt{plotType="M"}' and `\texttt{plotType="GC"}' generate plots of mean
ChIP tag counts versus mappability and GC content scores, respectively.
`\texttt{plotType="input"}' generates a plot of mean ChIP tag counts versus
control tag counts. Moreover, `\texttt{plotType="M|input"}' and
`\texttt{plotType="GC|input"}' generate plots of mean ChIP tag counts versus
mappability and GC content scores, respectively, conditional on control tag
counts. If `\texttt{plotType}' is not specified, this method plots the
histogram of ChIP tag counts.

<<ts-bindata-plot,eval=FALSE>>=
plot(exampleBinData)
plot( exampleBinData, plotType="M" )
plot( exampleBinData, plotType="GC" )
plot( exampleBinData, plotType="input" )
plot( exampleBinData, plotType="M|input" )
plot( exampleBinData, plotType="GC|input" )  
@

Figures \ref{fig:bindata-plot-hist}, \ref{fig:bindata-plot-M},
\ref{fig:bindata-plot-GC}, \ref{fig:bindata-plot-input},
\ref{fig:bindata-plot-M-input}, and \ref{fig:bindata-plot-GC-input} display
examples of different types of plots. As discussed in \cite{mosaics}, we
observe that mean ChIP tag count increases as mappability score increases
(Figure \ref{fig:bindata-plot-M}). Mean ChIP tag count depends on GC score
in a non-linear fashion (Figure \ref{fig:bindata-plot-GC}). The relationship
between mean ChIP tag counts and control tag counts seems to be linear,
especially for small  control tag counts (Figure \ref{fig:bindata-plot-input}).
When we condition on control tag counts (Figures \ref{fig:bindata-plot-M-input}
and \ref{fig:bindata-plot-GC-input} ), mean ChIP tag count versus mappability
and GC content relations exhibit similar patterns to that of marginal plots
given in Figures  \ref{fig:bindata-plot-M} and \ref{fig:bindata-plot-GC}.
MOSAiCS incorporates this observation by modeling ChIP tag counts from non-peak
regions  with a small number of control tag counts as a function of
mappability, GC content, and control tag counts.

\begin{figure}[tbh]
\begin{center}
<<fig-bindata-plot-hist,fig=TRUE,height=6,width=6,echo=FALSE>>=
plot(exampleBinData)
@
\caption{\label{fig:bindata-plot-hist} Histograms of the  count data 
from ChIP and control samples.}
\end{center}
\end{figure}

\begin{figure}[tbh]
\begin{center}
<<fig-bindata-plot-M,fig=TRUE,height=6,width=6,echo=FALSE>>=
plot( exampleBinData, plotType="M" )
@
\caption{\label{fig:bindata-plot-M} Mean ChIP tag count versus Mappability.}
\end{center}
\end{figure}

\begin{figure}[tbh]
\begin{center}
<<fig-bindata-plot-GC,fig=TRUE,height=6,width=6,echo=FALSE>>=
plot( exampleBinData, plotType="GC" )
@
\caption{\label{fig:bindata-plot-GC} Mean ChIP tag count versus GC content.}
\end{center}
\end{figure}

\begin{figure}[tbh]
\begin{center}
<<fig-bindata-plot-input,fig=TRUE,height=6,width=6,echo=FALSE>>=
plot( exampleBinData, plotType="input" )
@
\caption{\label{fig:bindata-plot-input} Mean ChIP tag count versus
Control tag count.}
\end{center}
\end{figure}

\begin{figure}[tbh]
\begin{center}
<<fig-bindata-plot-M-input,fig=TRUE,height=6,width=6,echo=FALSE>>=
plot( exampleBinData, plotType="M|input" )
@
\caption{\label{fig:bindata-plot-M-input} Mean ChIP tag count versus
Mappability, conditional on control tag counts.}
\end{center}
\end{figure}

\begin{figure}[tbh]
\begin{center}
<<fig-bindata-plot-GC-input,fig=TRUE,height=6,width=6,echo=FALSE>>=
plot( exampleBinData, plotType="GC|input" )
@
\caption{\label{fig:bindata-plot-GC-input} Mean ChIP tag count versus
GC content, conditional on control tag counts.}
\end{center}
\end{figure}

\subsection{Fitting MOSAiCS}

We are now ready to fit a MOSAiCS model using the bin-level data above
(\texttt{exampleBinData})   with the command:

<<ts-mosaicsfit,eval=FALSE>>=
exampleFit <- mosaicsFit( exampleBinData, analysisType="TS" )
@

<<ts-mosaicsfit-load-fit,echo=FALSE>>=
data(exampleFit)
@

`\texttt{analysisType="TS"}' indicates implementation of the two-sample
analysis.
`\texttt{mosaicsFit}' fits both one-signal-component and two-signal-component
models. When identifying peaks, you can choose the number of signal components
to be used for the final model. The optimal choice of the number of signal
components depends on the characteristics of data. In order to support users
in the choice of optimal signal model, \texttt{mosaics} package provides
Bayesian Information Criterion (BIC) values and Goodness of Fit (GOF) plots
of these signal models. 

The following command prints out BIC values of one-signal-component and
two-signal-component models, with additional information about the parameters
used in fitting the background (non-enriched) distribution. A lower BIC value
indicates a better model fit. For this dataset, we conclude that the
two-signal-component model has a lower BIC and hence it provides a better fit.

<<ts-mosaicsfit-show>>=
exampleFit
@

`\texttt{plot}' method provides the GOF plot.  This plots allows visual
comparisons of the fits of the background, one-signal-component, and
two-signal-component models with the actual data.  Figure
\ref{fig:mosaicsfit-plot} displays the GOF plot for our dataset and we 
conclude that the two-signal-component model provides a better fit as is also
supported by its lower BIC value compared to the one-signal component model. 


<<ts-mosaicsfit-plot,eval=FALSE>>=
plot(exampleFit)
@

\begin{figure}[tbh]
\begin{center}
<<fig-mosaicsfit-plot,fig=TRUE,height=6,width=6,echo=FALSE>>=
plot(exampleFit)
@
\caption{\label{fig:mosaicsfit-plot} Goodness of Fit (GOF) plot. Depicted are
actual data for ChIP and control samples with simulated data from the following
fitted models: (Sim:N): Background model; (Sim:N+S1): one-signal-component
model; (Sim:N+S1+S2): two-signal-component model. }
\end{center}
\end{figure}

In addition to `\texttt{analysisType}', `\texttt{mosaicsFit}' method provides
parameters to tune the background distribution of the MOSAiCS model. 
We specified appropriate default values for these parameters based on
computational experiments and analysis of diverse ChIP-seq datasets. Default
values work well in general but some tuning might be required for some cases.
You may need to consider parameter tuning if the fitted background model is 
too similar to the actual data in the GOF plot or you encounter some warning 
or error messages while running `\texttt{mosaicsFit}' method. Section
\ref{tuning} provides basic guidelines on parameter tuning. If you encounter
a fitting problem you need help with, feel free to contact us at our Google
group, \url{http://groups.google.com/group/mosaics_user_group}.

\subsection{Identifying Peaks Based on the Fitted Model}

Using BIC values and GOF plots in the previous section, we concluded that
two-signal-component model fits our data better. Next, we will identify peaks
with the two-signal-component model at a false discovery rate (FDR) of 0.05
using the command:

<<ts-mosaicspeak>>=
examplePeak <- mosaicsPeak( exampleFit, signalModel="2S", FDR=0.05, 
maxgap=200, minsize=50, thres=10 )
@

`\texttt{signalModel="2S"}' indicates two-signal-component model. Similarly,
one-signal-component model can be specified by `\texttt{signalModel="1S"}'. 
FDR can be controlled at the desired level by specifying `\texttt{FDR}'. 
In addition to these two essential parameters, you can also control three more
parameters, `\texttt{maxgap}', `\texttt{minsize}', and `\texttt{thres}'. These
parameters are for refining initial peaks called using specified signal model
and FDR.  Initial nearby peaks are merged if the distance (in bp) between them
is less than `\texttt{maxgap}'. Some initial peaks are removed if their lengths
are shorter than `\texttt{minsize}' or their ChIP tag counts are less 
than `\texttt{thres}'.

If you use a bin size shorter than the average fragment length in the
experiment, we recommend to set `\texttt{maxgap}'  to the  average fragment
length and `\texttt{minsize}'  to the bin size. This setting removes peaks 
that are too narrow (e.g., singletons). If you set the bin size to the average
fragment length (or maybe bin size is larger than the average fragment length),
we recommend  setting `\texttt{minsize}'  to a value  smaller than the average
fragment length while leaving `\texttt{maxgap}' the same as the average
fragment length. This is to prevent filtering using `\texttt{minsize}' because
initial peaks would already be at a reasonable width. `\texttt{thres}' is
employed to filter out initial peaks with very small ChIP tag counts because
such peaks might be false discoveries. Optimal choice of `\texttt{thres}'
depends on the sequencing depth of the ChIP-seq data to be analyzed. If you
don't wish to filter out initial peaks using ChIP tag counts, you can set
`\texttt{thres}'  to an arbitrary negative value.

The following command prints out a summary of identified peaks including 
the number of  peaks identified, median peak width, and the empirical false
discovery rate (FDR).

<<ts-mosaicspeak-show>>=
examplePeak
@

`\texttt{print}' method returns the peak calling results in data frame format.
This data frame can be used as an input for downstream analysis such as motif
finding. This output might have different number of columns, depending on
`\texttt{analysisType}' of `\texttt{mosaicsFit}'. For example, if
`\texttt{analysisType="TS"}', columns are peak start position, peak end
position, peak width, averaged posterior probability, minimum posterior
probability, averaged ChIP tag count, maximum ChIP tag count, averaged control
tag count, averaged control tag count scaled by sequencing depth, averaged log
base 2 ratio of ChIP over input tag counts, averaged mappability score, and
averaged GC content score for each peak. Here, the posterior probability of a
bin refers to  the probability that the bin is not a peak conditional on data.
Hence, smaller posterior probabilities provide more evidence that the bin is
actually a peak.

<<ts-mosaicspeak-print>>=
print(examplePeak)[1:15,]
@

You can export peak calling results to text files in diverse file formats.
Currently, `\texttt{mosaics}' package supports TXT, BED, and GFF file formats.
In the exported file, TXT file format (`\texttt{type="txt"}') includes all the
columns that `\texttt{print}' method returns. `\texttt{type="bed"}' and
`\texttt{type="gff"}' export peak calling results in standard BED and GFF file
formats, respectively, where score is the averaged ChIP tag counts in each
peak. Peak calling results can be exported in TXT, BED, and GFF file formats,
respectively, by the commands:

<<ts-mosaicspeak-export>>=
export( examplePeak, type="txt", fileLoc=".",
 fileName="TSpeakList.txt", chrID="chr21" )
export( examplePeak, type="bed", fileLoc=".",
 fileName="TSpeakList.bed", chrID="chr21" )
export( examplePeak, type="gff", fileLoc=".", 
 fileName="TSpeakList.gff", chrID="chr21" )
@

`\texttt{fileLoc}' and `\texttt{fileName}' indicate the directory and the name
of the exported file. `\texttt{chrID}' means chromosome ID and you need 
to specify this because `\texttt{mosaics}' package assumes chromosome-wise
analysis and does not require user to supply chromosome ID during the analysis.
The word specified in `\texttt{chrID}' will appear in the first column in the
exported file.

\section{One-Sample Analysis}

When control sample is not available, `\texttt{mosaics}' package accommodates
one-sample analysis of ChIP-seq data. Implementation of the MOSAiCS one-sample
model is very similar to that  of  the two-sample analysis. Bin-level data for
the one-sample analysis can be imported to the R environment with the command:

<<os-readbin>>=
OneSampleBinData <- readBins( type=c("chip","M","GC","N"),
    fileName=c( system.file( file.path("extdata","chip_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","M_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","GC_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","N_chr21.txt"), package="mosaicsExample") ) )
@

Note that you don't need to provide `\texttt{"input"}' in `\texttt{type}' and
the file name of  a control dataset in `\texttt{fileName}' here. In order to
fit a MOSAiCS model for the one-sample analysis, you need to specify
`\texttt{analysisType="OS"}' instead of `\texttt{analysisType="TS"}' 
when calling  the `\texttt{mosaicsFit}' method.

<<os-mosaicsfit,eval=FALSE>>=
OneSampleFit <- mosaicsFit( OneSampleBinData, analysisType="OS" )
@

Peak identification can be done exactly in the same way as in the case 
of the two-sample analysis.

<<os-mosaicspeak,eval=FALSE>>=
OneSamplePeak <- mosaicsPeak( OneSampleFit, signalModel="2S", FDR=0.05,
maxgap=200, minsize=50, thres=10 )
@

\section{Two-Sample Analysis Without Mappability and GC Content}

Application of MOSAiCS to multiple case studies showed that consideration of
mappability and GC content in the model improves sensitivity and specificity 
of peak identification even in the presence of a control sample \cite{mosaics}.
However, \texttt{mosaics} package accommodates a  two-sample analysis without
mappability and GC content by specification of `\texttt{analysisType="IO"}'
when calling the `\texttt{mosaicsFit}' method.

<<io-mosaicsfit,eval=FALSE>>=
inputOnlyFit <- mosaicsFit( exampleBinData, analysisType="IO" )
@

You can import bin-level data (for ChIP and control sample only) and fit
MOSAiCS model for the two-sample analysis without mappability and GC content
with the commands:

<<io-readbin>>=
inputOnlyBinData <- readBins( type=c("chip","input"),
    fileName=c( system.file( file.path("extdata","chip_chr21.txt"), package="mosaicsExample"),
    system.file( file.path("extdata","input_chr21.txt"), package="mosaicsExample") ) )
@
<<io-mosaicsfit-io-readbin,eval=FALSE>>=
inputOnlyFit <- mosaicsFit( inputOnlyBinData, analysisType="IO" )
@

\section{Tuning of MOSAiCS Parameters}\label{tuning}

In the two-sample analysis, users can control three tuning parameters:
`\texttt{s}', `\texttt{d}', and `\texttt{meanThres}'.  `\texttt{s}' and
`\texttt{d}'  are parameters of the background distribution and control the
functional form used for the control data.  Please see \cite{mosaics} 
for further details on these two model parameters. `\texttt{meanThres}' 
controls the number of strata used at the robust linear regression modelling
step of the background distribution fitting. `\texttt{mosaics}' package uses
the following parameter setting as default:

<<tuning-ts,eval=FALSE>>=
exampleFit <- mosaicsFit( exampleBinData, analysisType="TS", 
meanThres=1, s=2, d=0.25 )
@

Users might need to consider parameter tuning especially when the fitted
background model is too similar to the actual data, resulting in too few peaks.
If such cases are detected or predicted, `\texttt{mosaicsFit}' prints out
warning or error messages. You may also be able to detect this case using 
the GOF plot.  Using a higher `\texttt{s}' value and lower `\texttt{meanThres}'
often solves the problem, e.g., `\texttt{s = 6}' and `\texttt{meanThres = 0}'. 

\section{Conclusion and Ongoing Work}

R package \texttt{mosaics} provides effective tools to read and investigate
ChIP-seq data, fit MOSAiCS model, and identify peaks. We are continuously
working on improving \texttt{mosaics} package further, especially in supporting
more diverse genomes, automating fitting procedures, developing more friendly
and easy-to-use user interface, and providing more effective data investigation
tools. Please post any questions or requests regarding `\texttt{mosaics}'
package at \url{http://groups.google.com/group/mosaics_user_group}. Updates and
changes of `\texttt{mosaics}' package will be announced at our Google group and
the companion website (\url{http://www.stat.wisc.edu/~keles/Software/mosaics/}).

\begin{thebibliography}{99}
\bibitem{mosaics} Kuan, PF, D Chung, JA Thomson, R Stewart, and S Kele\c{s}
(2010), ``A Statistical Framework for the Analysis of ChIP-Seq Data'', submitted
(\url{http://works.bepress.com/sunduz_keles/19/}).
\bibitem{stat1} Rozowsky, J, G Euskirchen, R Auerbach, D Zhang, T Gibson, 
R Bjornson, N Carriero, M Snyder, and M Gerstein (2009), ``PeakSeq enables
systematic scoring of ChIP-Seq experiments relative to controls'',
\textit{Nature Biotechnology}, 27, 66-75.
\end{thebibliography}

\end{document}