%\VignetteIndexEntry{NarrowPeaks Vignette I. Intra-sample variability: Splitting and narrowing down ChIP-seq peaks in a single experiment.}
%\VignetteDepends{NarrowPeaks}
%\VignetteKeywords{Visualization,ChIPseq,Transcription,Genetics}
%\VignettePackage{NarrowPeaks}
\documentclass{article}

\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Rmethod}[1]{{\texttt{#1}}}
\newcommand{\Rcode}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textsf{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}
\newcommand{\NarrowPeaks}{\Rpackage{NarrowPeaks}}
\newcommand{\CSAR}{\Rpackage{CSAR}}


\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{natbib}

\begin{document}
%
\title{An Introduction to the {\NarrowPeaks} Package: \\ Narrowing Down Transcription Factor Binding Site Candidates from Functional Data}
\date{January, 2013}
\author{Pedro Madrigal\footnote{pm@engineering.com}}
\maketitle
\begin{center}
Department of Biometry and Bioinformatics, Institute of Plant Genetics\\Polish Academy of Sciences\\Poznan, Poland
\end{center}

<<setup, echo=FALSE, results=hide>>=
options(width=60)
options(continue=" ")
options(prompt="R> ")
@ 

\section{Introduction}
State-of-the-art bioinformatic algorithms, so-called peak finders (see references  \cite{Laajala2009}, \cite{Pepke2009} and \cite{Wilbanks2010}), are used to detect transcription factor binding sites in high-throughput chromatin immunoprecipitation followed by sequencing (ChIP-seq). The data analysis is usually based on peak-search criteria of the local maxima over enriched candidate regions. For purposes of computation several assumptions are made regarding the distribution of sample and control reads. 

It has been shown that, although most sites reported by peak finders could be narrowed down to 100-400bp using merely visual inspection, this reduction is not typically reflected by the regions provided by current methods, therefore degrading the resolution \cite{Rye2011}. It is widely accepted that the subdivision of long regions into distinct subpeaks can further help to recognize individual true peaks that were merged into a wide area of signal enrichment.

We present here the R package {\NarrowPeaks} \cite{MadrigalSubmitted} able process WIG format\footnote{One of the most popular formats for ChIP-seq data visualization is the \href{http://genome.ucsc.edu/FAQ/FAQformat}{wiggle track (WIG)}.} data, and analyze it based on the theory of Functional Principal Component Analysis (FPCA) \cite{Ramsay2005}. 

The aim of this novel approach is to extract the most significant ChIP-seq enriched regions according to their primary modes of variation in the binding score profiles. It allows the user of this package to discriminate between binding regions in close proximity and shorten the length of the transcription factor binding sites preserving the information present in the the dataset at a desired level of variance.



\section{Methods}

The functional version of PCA establishes a method for estimating orthogonal basis functions (principal components or \emph{eigenfunctions}) from functional data \cite{Ramsay2005}, in order to capture as much of the variation as possible in as few components as possible. We can highlight the genomic locations contributing to maximum variation (measured by an aproximation of the variance-covariance function) from a list of peaks of a ChIP-seq experiment. 

The proposed algorithm converts a continuous signal of enrichment (from a WIG file into CSAR binary format), and extracts signal profiles of candidate
transcription factor binding sites. 
Afterwards, it characterizes the binding signals via spline basis functions expansion.
Finally, functional PCA is performed in order to measure the variation of the ChIP-seq signal profiles under study. The output consists of a score-ranked list of sites according to their contribution to the total variation, which is accounted for by the trimmed (narrowed) principal components (estimated from the data).


\section{Example}

We will use the example data set included in the {\NarrowPeaks} package for this demonstration. The data represents a small subset of a WIG file storing continuous value scores based on a Poisson test \cite{Muino2011} for the chromosome 1 of \emph{Arabidopsis thaliana} \cite{Kaufmann2010}. 

First, we load the {\NarrowPeaks} package and the data \Rclass{NarrowPeaks-dataset}, which contains a subsample of first 49515 lines of the original WIG file for the full experiment.
Using the function \Rfunction{wig2CSARScore} a set of binary files is constructed storing the enrichment-score profiles.


<<write score binaries from a WIG file>>=  
library(NarrowPeaks)
data("NarrowPeaks-dataset")
head(wigfile_test)
writeLines(wigfile_test, con="wigfile.wig")
wigScores <- wig2CSARScore(wigfilename="wigfile.wig", nbchr = 1, 
chrle=c(30427671))
print(wigScores$infoscores$filenames)
@

<<setup, echo=FALSE, results=hide>>=
gc(reset=TRUE)
@

Next, the candidate binding site regions are extracted using the R/Bioconductor package {\CSAR} \cite{Muino2011}. CSAR predictions are contiguous genomic regions separated by a maximum allowed of \Robject{g} base pairs, and score enrichment values greater than \Robject{t}. Candidate regions are stored in a \Robject{GRanges} object (see Bioconductor package \Rpackage{GenomicRanges}).

<<extract candidate regions of enrichment in GRanges object>>=
library(CSAR)
candidates <- sigWin(experiment=wigScores$infoscores, t=1.0, g=30)
head(candidates)
@


If {\CSAR} \cite{Muino2011} is used first to analyze ChIP-seq data, from its results we can obtain the false discovery rate (FDR) for a given threshold. For example, for the complete experiment described in \cite{Kaufmann2010}, \Robject{t = 10.81} corresponds to FDR = 0.01 and \Robject{t = 6.78} corresponds to FDR = 0.1. \\

Now we want to narrow down the candidate sites obtaining shortened peaks with the function \Rfunction{narrowpeaks}, representing each candidate signal as a linear combination of \Robject{nbf} B-spline basis functions with no derivative penalization \cite{Ramsay2005}. We can specify the amount of miminum variance \Robject{pv} we want to describe in form of \Robject{npcomp} principal components, and establish a cutoff \Robject{pmaxscor} for trimming of scoring functions of the candidate sites \cite{MadrigalSubmitted}.

We will run the function for three different values of the cutoff: \Robject{pmaxscor = 0} (no cutoff), \Robject{pmaxscor = 3} (cutoff is at 3\% level of the maximum value relative to the scoring PCA functions) and  \Robject{pmaxscor = 100} (cutoff is at the maximum value relative to the scoring PCA functions).


<<narrow down enriched regions by funtional PCA>>=
shortpeaksP0 <- narrowpeaks(inputReg=candidates, 
scoresInfo=wigScores$infoscores, lmin=0, nbf=150, rpenalty=0, 
nderiv=0, npcomp=2, pv=80, pmaxscor=0.0, ms=0)
head(shortpeaksP0$broadPeaks)
head(shortpeaksP0$narrowPeaks)
shortpeaksP3 <- narrowpeaks(inputReg=candidates, 
scoresInfo=wigScores$infoscores, lmin=0, nbf=150, rpenalty=0,
 nderiv=0, npcomp=2, pv=80, pmaxscor=3.0, ms=0)
head(shortpeaksP3$broadPeaks)
head(shortpeaksP3$narrowPeaks)
shortpeaksP100 <- narrowpeaks(inputReg=candidates, 
scoresInfo=wigScores$infoscores, lmin=0, nbf=150, rpenalty=0, 
nderiv=0, npcomp=2, pv=80, pmaxscor=100, ms=0)
head(shortpeaksP100$broadPeaks)
head(shortpeaksP100$narrowPeaks)
@

As we can see, there is no difference between \Robject{broadPeaks} and \Robject{narrowPeaks} for \Robject{pmaxscor = 0}, whereas for \Robject{pmaxscor = 100} just one punctual source of variation is reported. 
The number of components (\Robject{reqcomp}) required, as well as the variance (\Robject{pvar}) achieved, are the same for all three cases (\Robject{pmaxscor} of 0, 3 and 100). \\

<<final number of components and variance>>=
print(shortpeaksP0$reqcomp)
print(shortpeaksP0$pvar)
@

Now, we can do the same for \Robject{pmaxscor = 90} and the result consists of 3 peaks very close to each other. We can tune the parameter \Robject{ms} to merge the sites into a unique peak:

<<merge neighbouring narrow peaks>>=
shortpeaksP90 <- narrowpeaks(inputReg=candidates,
 scoresInfo=wigScores$infoscores, lmin=0, nbf=150, rpenalty=0, 
nderiv=0, npcomp=2, pv=80, pmaxscor=90, ms=0)
shortpeaksP90ms20 <- narrowpeaks(inputReg=candidates,
 scoresInfo=wigScores$infoscores, lmin=0, nbf=150, rpenalty=0,
 nderiv=0, npcomp=2, pv=80, pmaxscor=90, ms=20)
@

We can make use of the class \Rclass{GRangesLists} in the package \Rpackage{GenomicRanges} to create a compound structure:

<<create GRangesLists>>=
library(GenomicRanges)
exampleMerge <- GRangesList("narrowpeaksP90"=shortpeaksP90$narrowPeaks,
 "narrowpeaksP90ms20"=shortpeaksP90ms20$narrowPeaks);
exampleMerge
@

Finally, we can export \Rclass{GRanges} objects or \Rclass{GRangesLists} into WIG, bedGraph, bigWig or other format files using the package \Rpackage{rtracklayer}. For example:

<<export GRanges to annotation tracks in various formats>>=
library(GenomicRanges)
names(elementMetadata(shortpeaksP3$broadPeaks))[3] <- "score"
names(elementMetadata(shortpeaksP3$narrowPeaks))[2] <- "score"
library(rtracklayer)
export.bedGraph(object=candidates, con="CSAR.bed")
export.bedGraph(object=shortpeaksP3$broadPeaks, con="broadPeaks.bed")
export.bedGraph(object=shortpeaksP3$narrowPeaks, con="narrowpeaks.bed")
@


\bibliography{NarrowPeaksVignette}{}
\bibliographystyle{plain}

\section{Details}

This document was written using:

<<sessionInfo>>=
sessionInfo()
@ 

\end{document}