getPeakReads: Get reads from indexed bam files for defined regions

In dhelekal/mmdiff3: Statistical Testing for ChIP-Seq data sets

Description

This function collects all short reads from bam files that map to pre-defined regions of interest. Note, that it fetches the exact start and end positions of mapped fragments, not the coverage. In the case of single-end reads, the left most postions of fragments mapping to the positive strands and the right most positions of fragments mapping to the negative strands are stored. To find centers of fragments use estimateFragmentCenters(). Positions are given relative to the start of the peak. Also computed are TotalCounts, i.e. number of fragments mapping to a peak region, as well as number of fragments mapping to forward and reverse strand.

Usage

getPeakReads(MD, PeakBoundary = 200, strand.specific = FALSE,
  filter.duplicates = FALSE, run.parallel = FALSE, BPPARAM = bpparam())

Arguments

MD	DBAmmd Object. This Object can be created using DBAmmd().
PeakBoundary	Defines flanking regions of peaks. The estimated fragment length is a good guess (DEFAULT: 200).
strand.specific	if strand specific reads should be kept apart. (DEFAULT:FALSE)
filter.duplicates	wether duplicate fragments should be removed. (DEFAULT: FALSE)
run.parallel	whether to run in parallel (currently no parallelization implemented) (DEFAULT: FALSE)
BPPARAM	an optional parameter object passed internally to bplapply when parallel=TRUE. If not specified, the parameters last registered with register will be used.

Value

DBAmmd object with updated slots

See Also

DBAmmd, estimateFragmentCenters

Examples

## Example using a small data set provided in the MMDiffBamSubset package

# setting the Experiment meta data:
ExpData <- list(dataDir=system.file("extdata", package="MMDiffBamSubset"),
           sampleSheet="Cfp1.csv")

MetaData <- list('ExpData' = ExpData)

# Creating a DBAmmd data set:
MMD <- DBAmmd(MetaData)

# defining a small Region for which to get reads:
Regions <- GRanges(seqnames=c('chr1'),
           IRanges(start = c(4560912,4677889), end = c(4562680,4679681)))
MMD <- setRegions(MMD,Regions)
MMD <- getPeakReads(MMD)

# To access Left ends of fragments mapping to positive strand:
Reads.L <- Reads(MMD,'Left.p')

# To access Right ends of fragments mapping to negative strand:
Reads.R <- Reads(MMD,'Right.n')

# To access Matrix of TotalCounts:
C.t <- Counts(MMD,whichCounts='T')

# Counts on positive strand:
C.p <- Counts(MMD,whichCounts='p')

# Counts on negative strand:
C.n <- Counts(MMD,whichCounts='n')

dhelekal/mmdiff3 documentation built on July 25, 2019, 8:18 p.m.