NarrowPeaks: Detect Narrow Peaks from Enrichment-Score Profiles forming Broad Peaks

Description

Detect narrow peaks from enrichment-score profiles forming broad peaks.

Usage

narrowpeaks(inputReg, scoresInfo, lmin = 0, nbf = 50, rpenalty= 0,  
nderiv= 0, npcomp = 5, pv = 80, pmaxscor = 0.0, ms = 0)

Arguments

inputReg

Output of the function sigWin in package CSAR.

scoresInfo

Output infoscores in the function wig2CSARScore, or the function ChIPseqScore after data analysis with package CSAR.

lmin

Minimum length of an enriched region from the WIG file to be processed. Integer value.

nbf

Number of order 4 B-spline basis functions that will represent the shape of each candidate site. Integer value.

rpenalty

Smoothing parameter for derivative penalization. Positive numeric value.

nderiv

Order of derivative penalization, if rpenalty>0. Integer value.

npcomp

Number of functional principal components. Integer value greater than or equal to nbf.

Minimum percentage of variation to take into account during the analysis. Numeric value in the range 0-100 (see the vignette and the paper).

pmaxscor

Cutoff for trimming of scoring function. Numeric value in the range 0-100.

Peaks closer to each other than ms nucleotides are to be merged in the final list. Integer value.

Value

fdaprofiles: A functional data object encapsulating the enrichment profiles (see fda package. To plot the data use plot.fd(fdaprofiles)).
broadPeaks: Description of the peaks prior to trimming. A GRanges object (see GenomicRanges package) with the information: seqnames (chromosome), ranges (start and end of the candidate site), strand (not used), max (maximum signal value for candidate site), average (mean signal value for candidate site), fpcaScore (sum of squares of the first reqcomp principal component scores for candidate site).
narrowPeaks: Description of the peaks after trimming. A GRanges object (see GenomicRanges package) with the information: seqnames (chromosome), ranges (start and end after trimming), strand (not used), broadPeak.subpeak, trimmedScore (see details), narrowedDownTo (length reduction relative to the candidate), merged (logical value).
reqcomp: Number of functional principal components used. Integer value.
pvar: Total proportion of variance accounted for by the reqcomp components used. Numeric value in the range 0-100 (always greater than or equal to argument pv).

Details

This function produces shortened sites from a list of candidate transcription factor binding sites of arbitrary extension and shape. First, the enrichment signal from each candidate site is represented by a smoothed function constructed using a linear combination of order 4 B-spline basis functions. The data values are fitted using either least squares (if $rpenalty = 0$), or penalized residuals sum of squares (spline smoothing if $rpenalty > 0$). Then, a functional principal component analysis for npcomp eigenfunctions is performed (Ramsay and Silverman, 2005), giving as a result a set of probe scores (principal component scores) which sum of squares is reported in elementMetadata(broadPeaks)[,"fpcaScore"]. The higher the value of fpcaScore, the higher the variance that candidate peak accounts for within the original data. Details on the usage of semi-metrics in functional PCA is described in Ferraty and Vieu, 2006. After that, we impose the condition that total scoring function for each reported narrow peak must be at least pmaxscor per cent of the maximum value. Max value is calculated from a set of scoring functions using only the eigenfunctions required to achieve pv percent of variance. A new set of scores is computed using trimmed versions of the eigenfunctions (Madrigal and Krajewski, submitted), and the root square is stored in elementMetadata(narrowPeaks)[,"trimmedScore"].

References

Madrigal P, Krajewski P (in preparation) Shape-based Dimensionality Reduction Analyses by Functional PCA Reveal Associations between First and Higher Order Components in Next-Generation Sequencing Coverage Profiles. Bailey T, Krajewski P, Ladunga I, Lefebvre C, Li Q, Liu T, Madrigal P, Taslim C, Zhang J (2013) Practical guidelines for the comprehensive analysis of ChIP-seq data. PLoS Comput Biol 9(11):e1003326. Muino JM, Kaufmann K, van Ham RC, Angenent GC, Krajewski P (2011) ChIP-seq analysis in R (CSAR): An R package for the statistical detection of protein-bound genomic regions. Plant Methods 7:11. Ramsay, J.O. and Silverman, B.W. (2005) Functional Data Analysis. New York: Springer. Ferraty, F. and Vieu, P. (2006) Nonparametric Functional Data Analysis. New York: Springer.

Examples

Run this code

owd <- setwd(tempdir())

##For this example we will use a subset of the AP1 ChIP-Seq data (Kaufmann et
##al., 2010)
##The data is obtained after analysis using the CSAR package available in 
##Bioconductor 
data("NarrowPeaks-dataset")
writeLines(wigfile_test, con="wigfile.wig")

##Write binary files with the WIG signal values for each chromosome 
##independently and obtain regions of read-enrichment with score values greater
##than 't', allowing a gap of 'g'. Data correspond to enriched regions found up
##to 105Kb in the Arabidopsis thaliana genome
wigScores <- wig2CSARScore(wigfilename="wigfile.wig", nbchr = 1, 
chrle=c(30427671))
gc(reset=TRUE) 
library(CSAR)
candidates <- sigWin(experiment=wigScores$infoscores, t=1.0, g=30)

##Narrow down ChIPSeq enriched regions by functional PCA
shortpeaks <- narrowpeaks(inputReg=candidates, 
scoresInfo=wigScores$infoscores, lmin=0, nbf=150, rpenalty=0, 
nderiv=0, npcomp=2, pv=80, pmaxscor=3.0, ms=0)

###Export GRanges object with the peaks to annotation tracks in various 
##formats. E.g.:
library(GenomicRanges)
names(elementMetadata(shortpeaks$broadPeaks))[3] <- "score"
names(elementMetadata(shortpeaks$narrowPeaks))[2] <- "score"
library(rtracklayer)
export.bedGraph(object=candidates, con="CSAR.bed")
export.bedGraph(object=shortpeaks$broadPeaks, con="broadPeaks.bed")
export.bedGraph(object=shortpeaks$narrowPeaks, con="narrowpeaks.bed")

setwd(owd)

Run the code above in your browser using DataLab