enrichedPeaks: Find peaks in sequencing experiments.

In htSeqTools: Quality Control, Visualization and Processing for High-Throughput Sequencing data

Description

Find peaks in significantly enriched regions found via enrichedRegions.

Usage

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

Arguments

regions	list 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 list 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.

See Also

enrichedRegions

Examples

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(as(sample1,'GRanges'))
coverinpeaks <- regionsCoverage(chr=seqnames(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

htSeqTools documentation built on May 6, 2019, 3:39 a.m.