HMMtBroadPeak: Call broad peak by Hidden Markov Models with t emissions

View source: R/HMMtBroadPeak.R

HMMtBroadPeakR Documentation

Call broad peak by Hidden Markov Models with t emissions

Description

Call very broad peaks for data such as LAD domains, NAD domains. Reads will be count by each bins. Only bins with at least given reads (defined by background parameter) for all samples (pool all reads for each bin) will be subsequently normalized. These bins will be first normalized to CPM (count per million) reads and then do log2 transform for the ratio over control with a pseudocount. The peaks were defined by running a hidden markov model over the normalized values (using the R-package HMMt).

Usage

HMMtBroadPeak(
  treatment,
  control,
  binSize = 5000,
  background = 10,
  pseudocount = 1,
  gapwidth = binSize,
  ...
)

Arguments

treatment, control

Bam file of treatments and controls. Make sure the index file keep same prefix name with bam file.

binSize

The size of bins for count

background

Only bins with at least background reads (pool all reads) will be subsequently normalized.

pseudocount

default 1.

gapwidth

The Ranges of peaks separated by a gap less than gapwith positions will be merged.

...

parameters passed to BaumWelchT except m (fixed to 2).

Value

a list with elements counts and peaks. Bothe counts and peaks are GRanges objects.

Examples

treatment <- system.file("extdata", "LB1.KD.chr1_1_5000000.bam",
                          package = "HMMtBroadPeak",
                          mustWork = TRUE)
## call peak without control
res <- HMMtBroadPeak(treatment)
## call peak with control
control   <- system.file("extdata", "LB1.WT.chr1_1_5000000.bam",
                          package = "HMMtBroadPeak",
                          mustWork = TRUE)
called <- HMMtBroadPeak(treatment, control)
called$peaks
plotPeaks(called, seqname="chr1")


jianhong/HMMtBroadPeak documentation built on July 12, 2024, 6:11 p.m.