compHists: Compute Peak histograms

Description

This function computes histograms at pre-defined regions (peaks) from mapped fragments, i.e. fragment counts at genomic position. Note, in contrast to genomic coverage or density maps, this function uses a single position per fragment (usually its center) rather than the whole extend of the fragment. This results in a significant increase in resolution. The parameter whichPos determines whether fragment centers, start or end positions should be considered ('Center','Left','Right'). Results are stored as a list in the Hists slot of the DBAmmd Object, with one entry per peak. For each peak i, a (n x L_i) matrix is generated, where n is the number of samples and L_i is the number of bins used to cover the extend of the peak. Note, L_i varies between peaks of different lengths.

Usage

compHists(MD, bin.length = 20, whichPos = "Center")

Arguments

DBAmmd Object. This Object can be created using DBAmmd().

bin.length

size of binning window (in bp) (DEFAULT: 20)

whichPos

specifies which relative positions of mapped fragments should to be considered. Can be one of: 'Left.p', 'Right.p', 'Right.p' and 'Left.n': Start and end positions of fragments mapping to positive or negative strand, respectively ('Right.p' and 'Left.n' are not available for single-end reads). Additionally inferred positions: 'Center.n','Center.p','Center','Left','Right'. (DEFAULT: 'Center')

Value

DBAmmd object with updated slot Hists

Examples

Run this code

## Example using a small data set provided with this package:

data("MMD")
bin.length <- 20
MMD.1 <- compHists(MMD,bin.length)

# use \\code{plotPeak()} to plot indivdual peaks:
Peak.id <- '6'
plotPeak(MMD.1, Peak.id=Peak.id)

# or explicitly using the histograms:
H <- Hists(MMD.1, whichPos='Center')
Sample <- 'WT.AB2'
Peak.idx <- match(Peak.id, names(Regions(MMD.1)))
plot(H[[Peak.idx]][Sample,],t='l')

# add peak cooridnates:
Peak <- Regions(MMD.1)[Peak.idx]
meta <- metaData(MMD.1)
PeakBoundary <- meta$AnaData$PeakBoundary
x.coords <- as.integer(colnames(H[[Peak.idx]])) + start(Peak) - PeakBoundary
plot(x.coords,H[[Peak.idx]]['WT.AB2',],t='l',
    xlab=names(H)[Peak.idx], ylab='counts', main=Sample)

Run the code above in your browser using DataLab

Data Engineering and BI courses are free this week!