Re-estimate peak centres from coverage

Source: R/centrePeaks.R

Use coverage to estimate peak centres

centrePeaks(x, y, ...)

# S4 method for class 'GRanges,BamFileList'
centrePeaks(
  x,
  y,
  f = c("weighted.cov", "mean", "median"),
  BPPARAM = bpparam(),
  ...
)

# S4 method for class 'GRanges,BamFile'
centrePeaks(x, y, ...)

# S4 method for class 'GRanges,BigWigFileList'
centrePeaks(
  x,
  y,
  f = c("weighted.cov", "mean", "median"),
  BPPARAM = bpparam(),
  ...
)

# S4 method for class 'GRanges,BigWigFile'
centrePeaks(x, y, ...)

# S4 method for class 'GRanges,character'
centrePeaks(x, y, ...)

Arguments

x

A set of GRanges representing peaks

y

A suitable set of files with methods defined

...

Used to pass arguments between methods

f

The function to use when estimating a combined peak centre

BPPARAM

An object of class BPPARAM

Value

A GRanges object with all widths set to one

Details

Use coverage to estimate the centre of a set of peaks or GenomicRanges.

If using the mean or median, the point of maximum coverage for each sample will be found within each peak and these positions will be averaged to return a position representing an estimated peak centre.

If using weighted.cov, positions are weighted by the combined coverage across all samples to return the weighted mean position. In this case coverage will be scaled by total alignments within each bam file before summing across files

Examples

## Define some peaks
f <- system.file("extdata/peaks.bed.gz", package = "extraChIPs")
peaks <- importPeaks(f, type = "bed")[[1]]
peaks
#> GRanges object with 5 ranges and 1 metadata column:
#>       seqnames              ranges strand |     score
#>          <Rle>           <IRanges>  <Rle> | <integer>
#>   [1]    chr10 103865473-103865960      * |         0
#>   [2]    chr10 103877005-103877682      * |         0
#>   [3]    chr10 103879286-103881391      * |         0
#>   [4]    chr10 103892382-103893412      * |         0
#>   [5]    chr10 103899873-103900502      * |         0
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths

## Use a bam file to re-centre the regions using highest coverage
bf <- system.file("extdata/bam/ex1.bam", package = "extraChIPs")
centres <- centrePeaks(peaks, bf, BPPARAM = SerialParam())
#> Getting coverage across all 5 peaks...
#> done (0.11s)
#> Finding centres by chromosome...
#> done (0.11s)
centres
#> GRanges object with 5 ranges and 1 metadata column:
#>         seqnames    ranges strand |     score
#>            <Rle> <IRanges>  <Rle> | <integer>
#>   chr10    chr10 103865730      * |         0
#>   chr10    chr10 103877320      * |         0
#>   chr10    chr10 103880358      * |         0
#>   chr10    chr10 103892713      * |         0
#>   chr10    chr10 103900167      * |         0
#>   -------
#>   seqinfo: 1 sequence from an unspecified genome; no seqlengths