islandCounts {htSeqTools} | R Documentation

## Find genomic regions with high coverage and count number of reads overlapping each region in each sample

### 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)

[Package htSeqTools version 1.28.3 Index]