enrichedPeaks {htSeqTools} | R Documentation

## 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` | `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

[Package htSeqTools version 1.28.3 Index]