islandCounts: Find genomic regions with high coverage and count number of...

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

Description

Finds genomics regions where the coverage is above a user-specified threshold and counts the number of ranges in each sample overlapping each region.

Usage

islandCounts(x, minReads=10, mc.cores=1)

Arguments

x	RangedData or list containing the reads. If a list is provided, the overall coverage across all its elements is used to find the regions of interest, but individual counts are computed for each element in the list.
minReads	Only regions with coverage above minReads are considered.
mc.cores	If mc.cores>1 computations are performed in parallel, using function mclapply from package parallel.

Details

The output of islandCounts can be the input data for a number of downstream analysis methods. Although for a simple-minded analysis one could use enrichedRegions, one will usually want to use more sofisticated analyses (e.g. from packages DEseq, BayesPeak, limma etc.)

Value

Object of class RangedData indicating the regions with coverage above minReads and the number of reads overlapping each sample for those regions.

Methods

signature(x = "RangedData")

x is assumed to contain the reads from a single sample. Genomic regions with high coverage will be detected and the number of reads overlapping these regions will be computed.

signature(x = "list")

x is assumed to contain the reads for several samples, one sample in each element of the list. The overall coverage across all samples is computed by adding the coverage in the individual samples, and the regions with overall coverage above the user-specified threshold are selected. Then the number of reads overlapping each region is computed.

Examples

set.seed(1)
st <- round(rnorm(1000,500,100))
strand <- rep(c('+','-'),each=500)
space <- rep('chr1',length(st))
sample1 <- RangedData(IRanges(st,st+38),strand=strand,space=space)
st <- round(rnorm(1000,1000,100))
sample2 <- RangedData(IRanges(st,st+38),strand=strand,space=space)

regions <- islandCounts(list(sample1,sample2),minReads=50)
regions

#Plot coverage
plot(coverage(as(sample1,'GRanges'))[[1]],type='l',xlim=c(0,2000))
lines(coverage(as(sample2,'GRanges'))[[1]],col=2)

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