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

###################################################
### code chunk number 1: setup
###################################################
library(RNAseqData.HNRNPC.bam.chr14)
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(Rsamtools)
library(ShortRead)
library(parallel); options(mc.cores=detectCores())
library(lattice)


###################################################
### code chunk number 2: restriction-what
###################################################
param <- ScanBamParam(what=c("rname", "pos", "cigar"))


###################################################
### code chunk number 3: restriction-which
###################################################
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
exByGn <- exonsBy(TxDb.Dmelanogaster.UCSC.dm3.ensGene, "gene")
seqlevels(exByGn, force=TRUE) <- "chr3L"
gns <- unlist(range(exByGn))
param <- ScanBamParam(which=gns)


###################################################
### code chunk number 4: restrictions-several
###################################################
param <- ScanBamParam(what=c("rname", "pos", "cigar"), which=gns)


###################################################
### code chunk number 5: bigdata
###################################################
bigdata <- "~/BigData"
bamfls <- dir(file.path(bigdata, "bam"), ".bam$", full=TRUE) #$
names(bamfls) <- sub("_subset.bam", "", basename(bamfls))


###################################################
### code chunk number 6: scanBam-restriction
###################################################
param <- ScanBamParam(what="rname")
bam <- scanBam(bamfls[1], param=param)[[1]]
table(bam$rname)


###################################################
### code chunk number 7: readGappedAlignments-restriciton
###################################################
gal <- readGappedAlignments(bamfls[1])
head(gal, 3)
param <- ScanBamParam(what="seq")  ## also input sequence
gal <- readGappedAlignments(bamfls[1], param=param)
head(mcols(gal)$seq)


###################################################
### code chunk number 8: GC-sampling
###################################################
library(Morgan2013)
gcFunction


###################################################
### code chunk number 9: fastq-files
###################################################
bigdata <- "~/BigData"
fqfl <- dir(file.path(bigdata, "fastq"), ".fastq$", full=TRUE) #$


###################################################
### code chunk number 10: fastq-sampling
###################################################
sampler <- FastqSampler(fqfl, 100000)
fq <- yield(sampler)                    # 100,000 reads
lattice::densityplot(gcFunction(sread(fq)), plot.points=FALSE)
fq <- yield(sampler)                    # a different 100,000 reads


###################################################
### code chunk number 11: fastq-paired-end (eval = FALSE)
###################################################
## ## NOT RUN
## set.seed(123)
## end1 <- yield(FastqSampler("end_1.fastq"))
## set.seed(123)
## end2 <- yield(FastqSampler("end_2.fastq"))


###################################################
### code chunk number 12: readLines-chunks (eval = FALSE)
###################################################
## ## NOT RUN
## con <- file("<hypothetical-file>.txt")
## open(f)
## while (length(x <- readLines(f, n=10000))) {
##     ## work on character vector 'x'
## }
## close(f)


###################################################
### code chunk number 13: iteration-bam
###################################################
library(RNAseqData.HNRNPC.bam.chr14)
bamfl <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]
countBam(bamfl)
bf <- BamFile(bamfl, yieldSize=200000) # could be larger, e.g., 2-10 million


###################################################
### code chunk number 14: BamFile-iter (eval = FALSE)
###################################################
## ## initialize, e.g., for step 3 ...
## open(bf)
## while (length(gal <- readGappedAlignments(bf))) {
##     ## step 2: do work...
##     ## step 3: aggregate results...
## }
## close(bf)


###################################################
### code chunk number 15: counter-iter
###################################################
counter <- 
    function(query, subject, ..., ignore.strand=TRUE)
    ## query: GRanges or GRangesList
    ## subject: GappedAlignments
{
    if (ignore.strand)
        strand(subject) <- "*"
    hits <- countOverlaps(subject, query)
    countOverlaps(query, subject[hits==1])
}


###################################################
### code chunk number 16: counter-roi
###################################################
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
query <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")


###################################################
### code chunk number 17: counter-work (eval = FALSE)
###################################################
## ## initialize, e.g., for step 3 ...
## open(bf)
## while (length(gal <- readGappedAlignments(bf))) {
##     ## step 2: do work...
##     count0 <- counter(query, gal, ignore.strand=TRUE)
##     ## step 3: aggregate results...
## }
## close(bf)


###################################################
### code chunk number 18: counter-work
###################################################
counter1 <- 
    function(bf, query, ...)
{
    ## initialize, e.g., for step 3 ...
    counts <- integer(length(query))
    open(bf)
    while (length(gal <- readGappedAlignments(bf, ...))) {
        ## step 2: do work...
        count0 <- counter(query, gal, ignore.strand=TRUE)
        ## step 3: aggregate results...
        counts <- counts + count0
    }
    close(bf)
    counts
}


###################################################
### code chunk number 19: counter1
###################################################
bf <- BamFile(bamfl, yieldSize=500000)
counts <- counter1(bf, query)


###################################################
### code chunk number 20: count-bfl-make
###################################################
fls <- RNAseqData.HNRNPC.bam.chr14_BAMFILES  # 8 BAM files
bamfls <- BamFileList(fls, yieldSize=500000) # use a larger yieldSize


###################################################
### code chunk number 21: count-bams (eval = FALSE)
###################################################
## counts <- simplify2array(lapply(bamfls, counter1, query))


###################################################
### code chunk number 22: count-bams-mclapply
###################################################
options(mc.cores=detectCores())         # use all cores
counts <- simplify2array(mclapply(bamfls, counter1, query))
head(counts[rowSums(counts) != 0,], 3)