enrichedPeaks: Find peaks in sequencing experiments.

Description

Find peaks in significantly enriched regions found via enrichedRegions.

Usage

enrichedPeaks(regions, sample1, sample2, minHeight=100, space, mc.cores=1)

Arguments

regions

RangedDataList or RangedData indicating the regions in which we wish to find peaks.

sample1

IRangesList or IRanges object containing start and end of sequences in sample 1.

sample2

Same for sample 2. May be left missing, in which case only sample1 is used to find peaks.

minHeight

If sample2 is missing, peaks are defined as regions where the coverage in sample1 is greater or equal than minHeight. If sample2 is specified, the difference of coverage in sample1 minus sample2 must be greater or equal than minHeight.

space

Character text giving the name of the space for the RangedData object. Only used if sample1 and sample2 are of class RangedData, for RangedDataList this is set up automatically.

mc.cores

If mc.cores>1 computations for each element in the IRangesList objects are performed in parallel (using the parallel function from package parallel). Notice: this option launches as many parallel processes as there are elements in x, which can place strong demands on the processor and memory.

Value

Object of class RangedData indicating peaks higher than minHeight. Only peaks overlapping with regions are reported. The maximum of the coverage in each selected peak is reported in the column height (coverage in sample1 - sample2 when sample2 is specified). The column region.pvalue returns the p-value associated to the region that the peak belongs to (i.e. it is inherited from regions). Therefore, notice that all peaks corresponding to a single region will present the same region.pvalue.

Methods

signature(regions = "RangedData", sample1 = "IRanges", sample2 = "IRanges"): sample1 indicates the start/end of reads in sample 1, and similarly for sample2. Only the subset of regions indicated by the argument space will be used.
signature(regions = "RangedData", sample1 = "IRanges", sample2 = "missing"): sample1 indicates the start/end of reads in sample 1, and similarly for sample2. Only the subset of regions indicated by the argument space will be used.
signature(regions = "RangedData", sample1 = "IRangesList", sample2 = "IRangesList"): regions contains the regions of interest, sample1 and sample2 the reads in sample 1 and sample 2, respectively. names(sample1) and names(sample2) must correspond to the space names used in regions.
signature(regions = "RangedData", sample1 = "IRangesList", sample2 = "missing"): regions contains the regions of interest, sample1 the reads in sample 1. names(sample1) must correspond to the space names used in regions.
signature(regions = "RangedData", sample1 = "RangedData", sample2 = "missing"): space(sample1) indicates the chromosome, and start(sample1) and end(sample1) indicate the start/end of the reads in sample 1.
signature(regions = "RangedData", sample1 = "RangedData", sample2 = "RangedData"): space(sample1) indicates the chromosome, and start(sample1) and end(sample1) indicate the start/end of the reads in sample 1. Similarly for sample2.

Examples

Run this code

set.seed(1)
st <- round(rnorm(1000,500,100))
strand <- rep(c('+','-'),each=500)
space <- rep('chr1',length(st))
sample1 <- RangedData(IRanges(st,st+38),strand=strand,space=space)
st <- round(runif(1000,1,1000))
sample2 <- RangedData(IRanges(st,st+38),strand=strand,space=space)

#Find enriched regions and call peaks
mappedreads <- c(sample1=nrow(sample1),sample2=nrow(sample2))
regions <- enrichedRegions(sample1,sample2,mappedreads=mappedreads,minReads=50)
peaks <- enrichedPeaks(regions,sample1=sample1,sample2=sample2,minHeight=50)
peaks <- peaks[width(peaks)>10,]
peaks

#Compute coverage in peaks
cover <- coverage(sample1)
coverinpeaks <- regionsCoverage(chr=space(peaks),start=start(peaks),end=end(peaks),cover=cover)

#Evaluate coverage in regular grid and plot
#Can be helpful fo clustering of peak profiles
coveringrid <- gridCoverage(coverinpeaks)
coveringrid
plot(coveringrid)

#Standardize peak profiles dividing by max coverage
stdcoveringrid <- stdGrid(coveringrid, colname='maxCov')
stdcoveringrid