%\VignetteIndexEntry{NarrowPeaks Vignette. Shape-based splitting and trimming ChIP-Seq peaks using functional Principal Components} %\VignetteDepends{} %\VignetteKeywords{Visualization,ChIPseq,Transcription,Genetics,Sequencing,HighThroughputSequencing} %\VignettePackage{NarrowPeaks} \documentclass{article} <>= BiocStyle::latex(use.unsrturl=FALSE) @ \usepackage{graphicx} \usepackage{hyperref} \usepackage{natbib} \begin{document} \title{An Introduction to the \Biocpkg{NarrowPeaks} Package:\\ Narrowing Down Transcription Factor Binding Site Candidates from ChIP-Seq using Functional PCA} \date{Created: January, 2013. Last modified: June, 2014. Compiled: \today} \author{Pedro Madrigal} \maketitle \begin{small} \begin{flushleft} $^{1}$ Department of Biometry and Bioinformatics, Institute of Plant Genetics, Polish Academy of Sciences, Poznan, Poland \\ $^{2}$ \underline{Current address}: Wellcome Trust Sanger Institute, Wellcome Trust Genome Campus, Hinxton, Cambridge, UK \\ $^{3}$ \underline{Current address}: Anne McLaren Laboratory for Regenerative Medicine, Department of Surgery, University of Cambridge, Cambridge, UK \end{flushleft} \end{small} \tableofcontents <>= options(width=100) options(continue=" ") options(prompt="R> ") @ \section{Goal of this Vignette} We will show how to split and trim ChIP-Seq peaks in a WIG file by means of functional PCA. \section{Introduction} Comprehensive ChIP-Seq data analyses are carried out by many software tools \cite{Bailey2013}. State-of-the-art bioinformatic algorithms so-called peak finders (e.g., \cite{Laajala2009},\cite{Pepke2009},\cite{Wilbanks2010}), are used to detect transcription factor binding sites (TFBSs) in chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-Seq). Data analysis is usually based on peak-search criteria of the local maxima over the read-enriched candidate regions. For computation purposes, several assumptions are made regarding the distribution of sample and control reads \cite{Bailey2013}. 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 recognizing \emph{bona fide} TFBSs that were merged into a wide area of signal enrichment (broad peak). \\ We present here the R package \Biocpkg{NarrowPeaks} \cite{MadrigalSubmitted} able process data in WIG format (one of the most popular standard formats for visualisation of next-generation sequencing data is the \href{http://genome.ucsc.edu/FAQ/FAQformat}{wiggle track (WIG), and its indexed version bigWig}) data, and analyze it based on statistics of Functional Principal Component Analysis (FPCA) \cite{Ramsay2005}. Instructions to create WIG/bigWig coverage tracks can be found in \file{Text S1} in \cite{Bailey2013} \\ The aim of this novel approach is to extract the most significant regions of ChIP-Seq peaks 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 user-defined level of variance. Without the trimming mode, it also serves to describe peak shapes using a set of statistics (FPCA scores) directly liked to the the principal components of the dataset. It facilitates the posterior clustering of the peaks by their peak-shape, as we have done also for histone modification ChIP-Seq and other next-generation sequencing data in \cite{MadrigalSubmitted}. \section{Methodology} 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 \cite{MadrigalSubmitted}. \\ The algorithm converts a continuous signal of enrichment from a WIG file, and extracts signal profiles of candidate TFBSs. Afterwards, it characterizes the binding signals via B-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). A more detailed description of the method is available in \cite{MadrigalSubmitted}. \section{Example} We will use the example data set included in the \Biocpkg{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 \Biocpkg{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. \begin{small} <>= 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) @ \end{small} Next, the candidate binding site regions are extracted using the Bioconductor package \Biocpkg{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 \Biocpkg{GenomicRanges}). \textbf{Alternatively, ChIP-Seq peaks obtained using other peak-callers can be provided building an analogous \Robject{GRanges} object}. In this case, the metadata 'score' must represent a numeric value directly proportional to the confidence of the peak ($p$-value) or the strength of the binding (fold-change). \begin{small} <>= library(CSAR) candidates <- sigWin(experiment=wigScores$infoscores, t=1.0, g=30) head(candidates) @ \end{small} If \Biocpkg{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 with the function \Rfunction{narrowpeaks} to obtain shortened peaks , 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). \begin{small} <>= shortpeaksP0 <- narrowpeaks(inputReg=candidates, scoresInfo=wigScores$infoscores, lmin=0, nbf=25, 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=25, 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=25, rpenalty=0, nderiv=0, npcomp=2, pv=80, pmaxscor=100, ms=0) head(shortpeaksP100$broadPeaks) head(shortpeaksP100$narrowPeaks) @ \end{small} 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). \begin{small} <>= print(shortpeaksP0$reqcomp) print(shortpeaksP0$pvar) @ \end{small} 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: \begin{small} <>= shortpeaksP90 <- narrowpeaks(inputReg=candidates,scoresInfo=wigScores$infoscores, lmin=0, nbf=25, rpenalty=0, nderiv=0, npcomp=2, pv=80, pmaxscor=90, ms=0) shortpeaksP90ms20 <- narrowpeaks(inputReg=candidates,scoresInfo=wigScores$infoscores, lmin=0, nbf=25, rpenalty=0,nderiv=0, npcomp=2, pv=80, pmaxscor=90, ms=20) @ \end{small} We can make use of the class \Rclass{GRangesLists} in the package \Biocpkg{GenomicRanges} to create a compound structure: \begin{small} <>= library(GenomicRanges) exampleMerge <- GRangesList("narrowpeaksP90"=shortpeaksP90$narrowPeaks, "narrowpeaksP90ms20"=shortpeaksP90ms20$narrowPeaks); exampleMerge @ \end{small} Finally, we can export \Rclass{GRanges} objects or \Rclass{GRangesLists} into WIG, bedGraph, bigWig or other format files using the package \Biocpkg{rtracklayer}. For example: \begin{small} <>= 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") @ \end{small} \bibliography{NarrowPeaksVignette}{} \bibliographystyle{plain} \newpage %\section{Figures} % %<>= %library(fda) %plot.fd(shortpeaksP0$fdaprofiles,xlab="nucleotides",cex.lab=0.5,cex.axis=0.5,cex.main=0.5, ylab="Read-enrichment profiles of candidate peaks", main="Figure 1") %@ % %%\newpage %%<>= %%library(fda) %%plot.fd(shortpeaksP0$fdaprofiles[1],cex.lab=0.5,cex.axis=0.5, xlab="Chr1", ylab="Read-enrichment profile of the candidate peak #1", main="Figure 2") %%@ % %\newpage %<>= %library(fda) %dd=as.data.frame(candidates) %y=eval.fd(evalarg=1:(max(dd[,4])-1), fdobj=shortpeaksP0$fdaprofiles[1]) %cent=round((19142+18996)/2) %plot((cent-round(max(dd[,4])/2)):round(cent+(max(dd[,4])/2)-1),y, type='l',cex.lab=0.5,cex.axis=0.5,cex.main=0.5, xlab="Chr1", ylab="Read-enrichment profile of the candidate peak #1", main="Figure 2") %abline(v=18996, col="red") %abline(v=19142, col="red") %@ \section{Details} This document was written using: \begin{small} <>= sessionInfo() @ \end{small} \section{Acknowledgements} This work was supported by the EU Marie Curie Initial Training Network SYSFLO (agreement number 237909). Special thanks to the users who've provided feedback about early versions of NarrowPeaks. \end{document}