inst/doc/Streamer.R

In Streamer: Enabling stream processing of large files

### R code from vignette source 'Streamer.Rnw'

###################################################
### code chunk number 1: Load packages
###################################################
library(GenomicAlignments)
library(Streamer)


###################################################
### code chunk number 2: BamInput class
###################################################

.BamInput <-
    setRefClass("BamInput",
    contains="Producer",
    fields=list(
        file="character",
        ranges="GRanges",
        .seqNames="character"))

.BamInput$methods(
    yield=function()
    {
       "yield data from .bam file"
       if (verbose) msg("BamInput$.yield()")
       if(length(.self$.seqNames))
       {
            seq <- .self$.seqNames[1]
            .self$.seqNames <- .self$.seqNames[-1]
            idx <- as.character(seqnames(.self$ranges)) == seq
            param <- ScanBamParam(which=.self$ranges[idx],
                                  what=character())
            aln <- readGAlignments(.self$file, param=param)
            seqlevels(aln) <- seq
       } else {
           aln <- GAlignments()
       }
       list(aln)
    })



###################################################
### code chunk number 3: BamInput constructor
###################################################

BamInput <- function(file, ranges,...)
{
    .seqNames <- names(scanBamHeader(file)[[1]]$target)
    .BamInput$new(file=file, ranges=ranges, .seqNames=.seqNames, ...)
}



###################################################
### code chunk number 4: CountGOverlap class
###################################################

.CountGOverlap <-
    setRefClass("CountGOverlap",
    contains="Consumer",
    fields=list(ranges="GRanges",
                mode="character",
                ignore.strand="logical"))

.CountGOverlap$methods(
    yield=function()
    {
        "return number of hits"
        if (verbose) msg(".CountGOt$yield()")
        aln <- callSuper()[[1]]
        df <- DataFrame(hits=numeric(0))
        if(length(aln))
        {
            idx <- as.character(seqnames(.self$ranges)) == levels(rname(aln))
            which <- .self$ranges[idx]
            olap <- summarizeOverlaps(which, aln, mode=.self$mode,
                                      ignore.strand=.self$ignore.strand)
            df <- as(assays(olap)[[1]], "DataFrame")
            dimnames(df) <- list(rownames(olap), seqlevels(aln))
        }
        df
    })

CountGOverlap <-
    function(ranges,
             mode = c("Union", "IntersectionStrict",
               "IntersectionNotEmpty"),
             ignore.strand = FALSE, ...)
{
    values(ranges)$pos <- seq_len(length(ranges))
    .CountGOverlap$new(ranges=ranges, mode=mode,
                       ignore.strand=ignore.strand, ...)
}



###################################################
### code chunk number 5: BAM file and ranges
###################################################

galn_file <- system.file("extdata", "ex1.bam", package="Rsamtools")
gr <-
    GRanges(seqnames =
            Rle(c("seq2", "seq2", "seq2", "seq1"), c(1, 3, 2, 4)),
            ranges = IRanges(rep(10,1), width = 1:10,
              names = head(letters,10)),
            strand = Rle(strand(rep("+", 5)), c(1, 2, 2, 3, 2)),
            score = 1:10,
            GC = seq(1, 0, length=10))

bam <- BamInput(file = galn_file, ranges = gr)
olap  <- CountGOverlap(ranges=gr, mode="IntersectionNotEmpty")
s <- Stream(bam, olap)
yield(s)



###################################################
### code chunk number 6: overlap counter
###################################################

overlapCounter <- function(pr, cs) {
    s <- Stream(pr, cs)
    len <- length(levels(seqnames(pr$ranges)))
    lst <- vector("list", len)
    for(i in 1:len) {
        lst[[i]] <- yield(s)
        names(lst[[i]]) <- "Count"
    }
    do.call(rbind, lst)[names(cs$ranges), ,drop=FALSE ]
}

bam <- BamInput(file = galn_file, ranges = gr)
olap  <- CountGOverlap(ranges=gr, mode="IntersectionNotEmpty")
overlapCounter(bam, olap)