gr.peaks: Find peaks in a 'GRanges' over a given meta-data field

In mskilab/skitools: Various mskilab R utilties

Find peaks in a GRanges over a given meta-data field

Description

Finds "peaks" in an input GRanges with value field y. first piles up ranges according to field score (default = 1 for each range) then finds peaks. If peel > 0, then recursively peels segments contributing to top peak, and recomputes nextpeak "peel" times if peel>0, bootstrap controls whether to bootstrap peak interval nbootstrap times if id.field is not NULL will peel off with respect to unique (sample) id of segment and not purely according to width if FUN preovided then will complex aggregating function of piled up values in dijoint intervals prior to computing "coverage" (FUN must take in a single argument and return a scalar) if id.field is not NULL, AGG.FUN is a second fun to aggregate values from id.field to output interval

Usage

gr.peaks(
  gr,
  field = "score",
  minima = FALSE,
  peel = 0,
  id.field = NULL,
  bootstrap = TRUE,
  na.rm = TRUE,
  pbootstrap = 0.95,
  nbootstrap = 10000,
  FUN = NULL,
  AGG.FUN = sum,
  peel.gr = NULL,
  score.only = FALSE,
  verbose = peel > 0
)

Arguments

gr	GRanges with some meta-data field to find peaks on
field	character field specifying metadata to find peaks on, default "score, can be NULL in which case the count is computed
minima	logical flag whether to find minima or maxima
id.field	character denoting field whose values specifyx individual tracks (e.g. samples)
bootstrap	logical flag specifying whether to bootstrap "peel off" to statistically determine peak boundaries
na.rm	remove NA from data
pbootstrap	quantile of bootstrap boundaries to include in the robust peak boundary estimate (i.e. essentially specifies confidence interval)
FUN	function to apply to compute score for a single individual
AGG.FUN	function to aggregate scores across individuals
nboostrap	number of bootstraps to run

Examples

## outputs example gene rich hotspots from example_genes GRanges
pk = gr.peaks(example_genes)

## now add a numeric quantity to example_genes and compute
## peaks with respect to a numeric scores, e.g. "exon_density"
example_genes$exon_density = example_genes$exonCount / width(example_genes)
pk = gr.peaks(example_genes, field = 'exon_density')

## can quickly find out what genes lie in the top peaks by agggregating back with
## original example_genes
pk[1:10] %$% example_genes[, 'name']

mskilab/skitools documentation built on Aug. 31, 2023, 1:13 p.m.