findPeakHeight: FDR and optimal minHeight value estimation for ChIP-Seq peak...

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

Description

Uses False Discovery Rate to estimate the optimal value of the minHeight parameter in the call to enrichedPeaks. False Discovery Rate is computed by swapping IP and Input samples and calculating the ratio of 'false' peaks identified for a given set of minHeight values.

Usage

findPeakHeight(regions, sample1, sample2, hmin=5, hmax=200, myfdr=0.01, gridSize=25, space, mc.cores=1)
plotminHeight(x,...)

Arguments

regions	RangedData indicating the regions in which we wish to find peaks.
sample1	IRangesList or Rle object containing start and end of sequences in sample 1 (IP sample) or their coverage respectively.
sample2	Same for sample 2 (Control sample).
hmin	Minimum minHeight value to be considered for FDR estimation. Defaults to 5.
hmax	Maximum minHeight value to be considered for FDR estimation. Max coverage difference between sample 1 and sample 2 will be also calculated, and the minimum of the two will be used as hmax. This is done to avoid a skewed distribution of minHeight values to test for FDR estimation.
myfdr	Desired FDR cut-off.
gridSize	Number of intermediate steps of minHeight threshold to consider between hmin and hmax. Since FDR and optimal minHeight estimation is done by actually performing peak calls, selecting a high value for gridSize can come at a big computational cost. Default value of 25 or close is recommended.
space	Character text giving the name of the space for the RangedData object. Only used if sample1 and sample2 are of class Rle, for RangedData and IRangesList 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.
x	An object of class list as returned by the call to findPeakHeight.
...	Other graphical arguments passed to function plot.

Value

Object of class list with slots fdr, npeaks indicating peak calling FDR and number of peaks identified in the IP sample for each considered minHeight value and cut, opt with the desired FDR cut-off and the correponding minHeight value to be used in the call to enrichedPeaks.

Methods

signature(regions = "RangedData", sample1 = "RangedData", sample2 = "RangedData")

sample1 indicates the start/end of reads in sample 1 (IP Sample), and similarly for sample2 (Control sample). 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 = "Rle", sample2 = "Rle")

regions contains the regions of interest, sample1 and sample2 the coverage for the reads in sample 1 and sample 2, respectively. names(sample1) and names(sample2) must correspond to the space names used in regions.

See Also

enrichedPeaks

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)
minHeight <- findPeakHeight(regions,sample1=sample1,sample2=sample2)
plotminHeight(minHeight)
peaks <- enrichedPeaks(regions,sample1=sample1,sample2=sample2,minHeight=minHeight$opt)
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
plot(stdcoveringrid)

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