inst/doc/GenomicFiles.R

In GenomicFiles: Distributed computing by file or by range

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

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: install (eval = FALSE)
###################################################
## if (!require("BiocManager"))
##     install.packages("BiocManager")
## BiocManager::install("GenomicFiles")


###################################################
### code chunk number 3: quick_start-load
###################################################
library(GenomicFiles)


###################################################
### code chunk number 4: quick_start-ranges
###################################################
gr <- GRanges("chr14", IRanges(c(19411500 + (1:5)*20), width=10))


###################################################
### code chunk number 5: class-bam-data
###################################################
library(RNAseqData.HNRNPC.bam.chr14)
fls <- RNAseqData.HNRNPC.bam.chr14_BAMFILES


###################################################
### code chunk number 6: quick_start-MAP
###################################################
MAP <- function(range, file, ...) {
    requireNamespace("Rsamtools")
    Rsamtools::pileup(file, scanBamParam=Rsamtools::ScanBamParam(which=range))
}


###################################################
### code chunk number 7: quick_start-reduceByFile
###################################################
se <- reduceByFile(gr, fls, MAP, summarize=TRUE)
se


###################################################
### code chunk number 8: quick_start-assays
###################################################
dim(assays(se)$data)  ## ranges x files


###################################################
### code chunk number 9: quick_start-MAP-REDUCE-reduceByRange
###################################################
REDUCE <- function(mapped, ...) {
    cmb = do.call(rbind, mapped)
    xtabs(count ~ pos + nucleotide, cmb)
}

lst <- reduceByRange(gr, fls, MAP, REDUCE, iterate=FALSE)


###################################################
### code chunk number 10: quick_start-result
###################################################
head(lst[[1]], 3)


###################################################
### code chunk number 11: overview-GenomicFiles
###################################################
GenomicFiles(gr, fls)


###################################################
### code chunk number 12: pileups-ranges
###################################################
gr <- GRanges("chr14", IRanges(c(19411677, 19659063, 105421963, 
              105613740), width=20))


###################################################
### code chunk number 13: pileups-MAP
###################################################
MAP <- function(range, file, ...) {
    requireNamespace("deepSNV")
    ct = deepSNV::bam2R(file, 
                        GenomeInfoDb::seqlevels(range), 
                        GenomicRanges::start(range), 
                        GenomicRanges::end(range), q=0)
    ct[, c("A", "T", "C", "G", "a", "t", "c", "g")]
}


###################################################
### code chunk number 14: pileups-REDUCE
###################################################
REDUCE <- function(mapped, ...)
    Reduce("+", mapped)


###################################################
### code chunk number 15: pileups-reduceByRange
###################################################
pile2 <- reduceByRange(gr, fls, MAP, REDUCE) 
length(pile2)
elementNROWS(pile2)


###################################################
### code chunk number 16: pileups-res
###################################################
head(pile2[[1]])


###################################################
### code chunk number 17: ttest-ranges
###################################################
roi <- GRanges("chr14", IRanges(c(19411677, 19659063, 105421963,
               105613740), width=20))


###################################################
### code chunk number 18: ttest-group
###################################################
grp <- factor(rep(c("A","B"), each=length(fls)/2))


###################################################
### code chunk number 19: ttest-MAP
###################################################
MAP <- function(range, file, ...) {
    requireNamespace("GenomicAlignments")
    param <- Rsamtools::ScanBamParam(which=range)
    as.numeric(unlist(
        GenomicAlignments::coverage(file, param=param)[range], use.names=FALSE))
}


###################################################
### code chunk number 20: ttest-REDUCE
###################################################
REDUCE <- function(mapped, ..., grp) {
    mat = simplify2array(mapped)
    idx = which(rowSums(mat) != 0)
    df = genefilter::rowttests(mat[idx,], grp)
    cbind(offset = idx - 1, df)
}


###################################################
### code chunk number 21: ttest-results (eval = FALSE)
###################################################
## ttest <- reduceByRange(roi, fls, MAP, REDUCE, iterate=FALSE, grp=grp)


###################################################
### code chunk number 22: junctions-ranges
###################################################
gr <- GRanges("chr14", IRanges(c(19100000, 106000000), width=1e7))


###################################################
### code chunk number 23: junctions-MAP
###################################################
MAP <- function(range, file, ...) {
    requireNamespace("GenomicAlignments") ## for readGAlignments()
                                          ## ScanBamParam()
    param = Rsamtools::ScanBamParam(which=range)
    gal = GenomicAlignments::readGAlignments(file, param=param)
    table(GenomicAlignments::njunc(gal))
} 


###################################################
### code chunk number 24: junctions-GenomicFiles
###################################################
gf <- GenomicFiles(gr, fls)
gf


###################################################
### code chunk number 25: junctions-counts1
###################################################
counts1 <- reduceByFile(gf[,1:3], MAP=MAP)
length(counts1)        ## 3 files
elementNROWS(counts1)  ## 2 ranges


###################################################
### code chunk number 26: junctions-counts1-show
###################################################
counts1[[1]]


###################################################
### code chunk number 27: junctions-REDUCE
###################################################
REDUCE <- function(mapped, ...)
    sum(sapply(mapped, "[", "1"))

reduceByFile(gr, fls, MAP, REDUCE)


###################################################
### code chunk number 28: junctions-counts2
###################################################
counts2 <- reduceFiles(gf[,1:3], MAP=MAP)


###################################################
### code chunk number 29: junctions-counts2-show
###################################################
## reduceFiles returns counts for all ranges.
counts2[[1]]
## reduceByFile returns counts for each range separately.
counts1[[1]]


###################################################
### code chunk number 30: coverage1-tiles
###################################################
chr14_seqlen <- seqlengths(seqinfo(BamFileList(fls))["chr14"])
tiles <- tileGenome(chr14_seqlen, ntile=5)


###################################################
### code chunk number 31: coverage1-tiles-show
###################################################
tiles


###################################################
### code chunk number 32: coverage1-MAP
###################################################
MAP = function(range, file, ...) {
    requireNamespace("GenomicAlignments")  ## for ScanBamParam() and coverage()
    param = Rsamtools::ScanBamParam(which=range)
    rle = GenomicAlignments::coverage(file, param=param)[range]
    c(width = GenomicRanges::width(range), 
      sum = sum(S4Vectors::runLength(rle) * S4Vectors::runValue(rle))) 
}


###################################################
### code chunk number 33: coverage1-REDUCE
###################################################
REDUCE = function(mapped, ...) {
    Reduce(function(i, j) Map("+", i, j), mapped)
}


###################################################
### code chunk number 34: coverage1-results (eval = FALSE)
###################################################
## cvg1 <- reduceByFile(tiles, fls, MAP, REDUCE, iterate=TRUE)


###################################################
### code chunk number 35: coverage2-MAP
###################################################
MAP = function(range, file, ...) {
    requireNamespace("GenomicAlignments")  ## for ScanBamParam() and coverage()
    GenomicAlignments::coverage(
        file,
        param=Rsamtools::ScanBamParam(which=range))[range]
}


###################################################
### code chunk number 36: coverage2-REDUCE
###################################################
REDUCE = function(mapped, ...) {
    sapply(mapped, Reduce, f = "+")
}


###################################################
### code chunk number 37: coverage2-results
###################################################
cvg2 <- reduceFiles(unlist(tiles), fls, MAP, REDUCE)


###################################################
### code chunk number 38: coverage2-show
###################################################
cvg2[1]


###################################################
### code chunk number 39: coverage3-MAP
###################################################
MAP = function(range, file, ...) {
    requireNamespace("BiocParallel")  ## for bplapply()
    nranges = 2
    idx = split(seq_along(range), ceiling(seq_along(range)/nranges))
    BiocParallel::bplapply(idx, 
        function(i, range, file) {
            requireNamespace("GenomicAlignments")  ## ScanBamParam(), coverage()
            chunk = range[i]
            param = Rsamtools::ScanBamParam(which=chunk)
            cvg = GenomicAlignments::coverage(file, param=param)[chunk]
            Reduce("+", cvg)            ## collapse coverage within chunks
        }, range, file)
}


###################################################
### code chunk number 40: coverage3-REDUCE
###################################################
REDUCE = function(mapped, ...) {
    sapply(mapped, Reduce, f = "+")
}


###################################################
### code chunk number 41: coverage3-results (eval = FALSE)
###################################################
## cvg3 <- reduceFiles(unlist(tiles), fls, MAP, REDUCE)


###################################################
### code chunk number 42: reduceByYield-YIELD
###################################################
library(GenomicAlignments)
bf <- BamFile(fls[1], yieldSize=100000)
YIELD <- function(x, ...) readGAlignments(x)


###################################################
### code chunk number 43: reduceByYield-MAP-REDUCE
###################################################
gr <- unlist(tiles, use.names=FALSE)
MAP <- function(value, gr, ...) {
    requireNamespace("GenomicRanges")  ## for countOverlaps()
    GenomicRanges::countOverlaps(gr, value)
}
REDUCE <- `+` 


###################################################
### code chunk number 44: reduceByYield-DONE
###################################################
DONE <- function(value) length(value) == 0L


###################################################
### code chunk number 45: reduceByYield-bplapply
###################################################
FUN <- function(file, gr, YIELD, MAP, REDUCE, tiles, ...) {
    requireNamespace("GenomicAlignments")  ## for BamFile, readGAlignments()
    requireNamespace("GenomicFiles")       ## for reduceByYield()
    gr <- unlist(tiles, use.names=FALSE)
    bf <- Rsamtools::BamFile(file, yieldSize=100000)
    YIELD <- function(x, ...) GenomicAlignments::readGAlignments(x)
    MAP <- function(value, gr, ...) {
        requireNamespace("GenomicRanges")  ## for countOverlaps()
        GenomicRanges::countOverlaps(gr, value)
    }
    REDUCE <- `+` 
    GenomicFiles::reduceByYield(bf, YIELD, MAP, REDUCE, gr=gr, parallel=FALSE)
}


###################################################
### code chunk number 46: sessionInfo
###################################################
toLatex(sessionInfo())