alignPeaks: Align peaks in a ChIP-Seq experiment by removing the strand...

In htSeqTools: Quality Control, Visualization and Processing for High-Throughput Sequencing data

Description

Align peaks in a ChIP-Seq experiment by removing the shift between reads aligned to the plus and the minus strands.

Usage

alignPeaks(x, strand, npeaks = 1000, bandwidth = 150, mc.cores=1)

Arguments

x	A list, RangedData or an IRangesList object containing the aligned reads in each chromosome.
strand	Strand that each read was aligned to. If x is of class list, strand can be a character vector of length 1 indicating the name of the field in x indicating the strand, i.e. x[[1]][[strand]] contains the strand information.
npeaks	Number of peaks to be used to estimate the shift size.
bandwidth	Only reads with distance less than bandwidth between them and their closest gene are used to estimate the shift size.
mc.cores	Number of cores to be used for parallel computing (passed on to mclapply). Only used if x is of class list.

Details

The procedure detects the npeaks highest peaks (using reads from both strands simultaneously). Then it selects reads which are less than bandwidth base pairs away from any of the peaks. Then it computes (a) the average distance between reads on the plus strand and the closest peak, (b) the same distance for reads on the minus strand. The mean difference between (a) and (b) is the estimated shift size. Reads on the plus strand are shifted to the right, whereas reads on the minus strands are shifted to the left.

Value

A CompressedIRangesList object with all reads shifted so that the strand specific bias is no longer present.

Methods

signature(x = "IRangesList", strand = "list")

Each element in x corresponds to a chromosome, and each range gives the start/end of a sequence. strand indicates the strand for the ranges in x.

signature(x = "RangedData", strand = "character")

x gives read start and end positions, and strand gives the name of the variable in values(x) containing the strand information.

signature(x = "list", strand = "character")

The method for RangedData is applied to each element in x separately, as each element may have a different strand-specific bias.

Examples

#Generate 1000 reads containing strand-specific bias
st <- runif(1000,1,250)
strand <- rep(c('+','-'),each=500)
st[strand=='-'] <- st[strand=='-'] + runif(500,50,100)
x <- RangedData(IRanges(st,st+38),strand=strand)
#Estimate and remove the bias
xalign <- alignPeaks(x, strand='strand', npeaks=1)

htSeqTools documentation built on May 6, 2019, 3:39 a.m.