### R code from vignette source 'vignettes/GenomicFiles/inst/doc/GenomicFiles.Rnw'

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


###################################################
### code chunk number 2: biocLite (eval = FALSE)
###################################################
## source("http://bioconductor.org/biocLite.R")
## biocLite("GenomicFiles")


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


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


###################################################
### code chunk number 5: class-slots
###################################################
getSlots("BamFileViews")


###################################################
### code chunk number 6: class-slots
###################################################
fileRange(bfv) <- GRanges("chr14", IRanges(c(1.9e7, 2e7), width = 900000))


###################################################
### code chunk number 7: pileups-new-ranges
###################################################
ranges1 <- GRanges("chr14", IRanges(c(19411677, 19659063, 105421963, 
                   105613740), width = 20))
bfv1 <- BamFileViews(fls, fileRange = ranges1)


###################################################
### code chunk number 8: pileups-MAP1
###################################################
MAP1 <- function(FILE, RANGE, ...) {
    require(deepSNV)
    ct <- bam2R(path(FILE), seqlevels(RANGE), start(RANGE), end(RANGE), q=0)
    ct[, c("A", "T", "C", "G", "a", "t", "c", "g")]
}


###################################################
### code chunk number 9: pileups-REDUCE1
###################################################
REDUCE1 <- function(MAPPED, ...) {
    Reduce("+", MAPPED)
}


###################################################
### code chunk number 10: pileups-reduceByRange
###################################################
res1 <- reduceByRange(bfv1, MAP1, REDUCE1) 
length(res1)


###################################################
### code chunk number 11: pileups-res1
###################################################
res1[[1]]


###################################################
### code chunk number 12: pileups-BamFileViews2
###################################################
ranges2 <- GRanges("chr14", 
               IRanges(c(19411677, 19659063, 105421963, 105613740),
                       c(19411747, 19659135, 105422990, 105614142)))
bfv2 <- BamFileViews(fls, fileRange = ranges2)


###################################################
### code chunk number 13: pileups-MAP2
###################################################
MAP2 <- function(FILE, RANGE, ...) {
    require(deepSNV)
    ct <- bam2R(path(FILE), seqlevels(RANGE), start(RANGE), end(RANGE), q=0)
    ct[, c("INS", "DEL", "ins", "del")]
}


###################################################
### code chunk number 14: pileups-REDUCE2
###################################################
REDUCE2 <- function(MAPPED, ...) {
    sapply(MAPPED, colSums)
}


###################################################
### code chunk number 15: pileups-reduceByFile
###################################################
res2 <- reduceByFile(bfv2, MAP2, REDUCE2) 
length(res2)


###################################################
### code chunk number 16: pileups-res2
###################################################
res2[[1]]


###################################################
### code chunk number 17: prebuilt-intro
###################################################
fl <- system.file("tests", "test.bw", package = "rtracklayer")
gr <- GRanges(Rle(c("chr2", "chr19"), c(4, 2)),
              IRanges(1 + c(200, 250, 500, 550, 1450, 1750), width=100))
bwfv <- BigWigFileViews(c(fl, fl), fileRange=gr)


###################################################
### code chunk number 18: prebuilt-coverage (eval = FALSE)
###################################################
## covSE <- coverage(bwfv)


###################################################
### code chunk number 19: prebuilt-summary (eval = FALSE)
###################################################
## sumSE <- summary(bwfv, type="mean")


###################################################
### code chunk number 20: groupfiles-group
###################################################
grp <- factor(rep(c("A","B"), each=ncol(bfv1)/2))


###################################################
### code chunk number 21: groupfiles-map
###################################################
MAP <- function(RANGE, FILE, ...) {
    require(GenomicAlignments)
    stopifnot(length(RANGE) == 1)
    param <- ScanBamParam(which=RANGE)
    algns <- readGAlignments(path(FILE), param=param)
    as.numeric(coverage(algns)[RANGE][[1]])
}


###################################################
### code chunk number 22: groupfiles-reduce
###################################################
REDUCE <- function(MAPPED, ..., grp) {
    m <- simplify2array(MAPPED)
    idx <- which(rowSums(m) != 0)
    df <- genefilter::rowttests(m[idx,], grp)
    cbind(offset=idx - 1, df)
}


###################################################
### code chunk number 23: groupfile-reduce
###################################################
fileRange(bfv1)$result <- reduceByRange(bfv1, MAP, REDUCE, grp=grp)
fileRange(bfv1)
head(fileRange(bfv1)$result[[1]])


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