###################################################
### chunk number 1: options
###################################################
#line 33 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
options(width=60)


###################################################
### chunk number 2: preliminaries
###################################################
#line 37 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
library(Rsamtools)


###################################################
### chunk number 3: ScanBamParam
###################################################
#line 86 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
which <- RangesList(seq1=IRanges(1000, 2000),
                    seq2=IRanges(c(100, 1000), c(1000, 2000)))
what <- c("rname", "strand", "pos", "qwidth", "seq")
param <- ScanBamParam(which=which, what=what)


###################################################
### chunk number 4: scanBam
###################################################
#line 96 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
bamFile <- 
    system.file("extdata", "ex1.bam", package="Rsamtools")
bam <- scanBam(bamFile, param=param)


###################################################
### chunk number 5: scanBam-elts-1
###################################################
#line 105 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
class(bam)
names(bam)


###################################################
### chunk number 6: scanBam-elts-2
###################################################
#line 113 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
class(bam[[1]])
names(bam[[1]])


###################################################
### chunk number 7: scanBam-elts-class
###################################################
#line 119 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
sapply(bam[[1]], class)


###################################################
### chunk number 8: scanBam-to-IRanges
###################################################
#line 124 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
lst <- lapply(names(bam[[1]]), function(elt) {
    do.call(c, unname(lapply(bam, "[[", elt)))
})
names(lst) <- names(bam[[1]])


###################################################
### chunk number 9: lst-to-DataFrame
###################################################
#line 132 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
head(do.call("DataFrame", lst))


###################################################
### chunk number 10: indexed-file
###################################################
#line 138 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
list.files(dirname(bamFile), pattern="ex1.bam(.bai)?")


###################################################
### chunk number 11: bam-remote eval=FALSE
###################################################
## #line 158 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
## which <- RangesList("6"=IRanges(100000L, 110000L))
## param <- ScanBamParam(which=which)
## na19240bam <- scanBam(na19240url, param=param)


###################################################
### chunk number 12: summaryFunction
###################################################
#line 213 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
summaryFunction <- 
    function(seqname, bamFile, ...)
{
    param <- ScanBamParam(what=c('pos', 'qwidth'),
                          which=GRanges(seqname, IRanges(1, 1e7)),
                          flag=scanBamFlag(isUnmappedQuery=FALSE))
    x <- scanBam(bamFile, ..., param=param)[[1]]
    coverage(IRanges(x[["pos"]], width=x[["qwidth"]]))
}


###################################################
### chunk number 13: summaryByChromosome
###################################################
#line 228 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
seqnames <- paste("seq", 1:2, sep="")
cvg <- lapply(seqnames, summaryFunction, bamFile)
names(cvg) <- seqnames
cvg


###################################################
### chunk number 14: BamViews-parts
###################################################
#line 266 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
library(GenomicFeatures)
bamRanges <- local({
    fl <- system.file("extdata", "CaffeineTxdb.sqlite", 
                      package="Rsamtools")
    transcripts(loadFeatures(fl))
})
slxMaq09 <- local({
    fl <- system.file("extdata", "slxMaq09_urls.txt", 
                      package="Rsamtools")
    readLines(fl)
})


###################################################
### chunk number 15: BamViews-construct
###################################################
#line 286 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
bamExperiment <- 
    list(description="Caffeine metabolism views on 1000 genomes samples",
         created=date())
bv <- BamViews(slxMaq09, bamRanges=bamRanges, 
               bamExperiment=bamExperiment)
metadata(bamSamples(bv)) <- 
    list(description="Solexa/MAQ samples, August 2009",
         created="Thu Mar 25 14:08:42 2010")


###################################################
### chunk number 16: BamViews-query
###################################################
#line 300 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
bamExperiment(bv)


###################################################
### chunk number 17: bamIndicies eval=FALSE
###################################################
## #line 312 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
## bamIndexDir <- tempfile()
## indexFiles <- paste(bamPaths(bv), "bai", sep=".")
## dir.create(bamIndexDir)
## idxFiles <- mapply(download.file, indexFiles,
##                    file.path(bamIndexDir, basename(indexFiles)) ,
##                    MoreArgs=list(method="curl"))


###################################################
### chunk number 18: readBamGappedAlignments eval=FALSE
###################################################
## #line 326 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
## library(multicore)
## olaps <- readBamGappedAlignments(bv)


###################################################
### chunk number 19: olaps
###################################################
#line 335 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
load(system.file("extdata", "olaps.Rda", package="Rsamtools"))
olaps
head(olaps[[1]])


###################################################
### chunk number 20: read-lengths
###################################################
#line 345 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
head(t(sapply(olaps, function(elt) range(qwidth(elt)))))


###################################################
### chunk number 21: focal
###################################################
#line 359 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
rng <- bamRanges(bv)[1]
strand(rng) <- "*"
olap1 <- endoapply(olaps, subsetByOverlaps, rng)
olap1
head(olap1[[24]])


###################################################
### chunk number 22: olap-cvg
###################################################
#line 369 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
minw <- min(sapply(olap1, function(elt) min(start(elt))))
maxw <- max(sapply(olap1, function(elt) max(end(elt))))
cvg <- endoapply(olap1, coverage,
                 shift=-start(ranges(bamRanges[1])),
                 width=width(ranges(bamRanges[1])))
cvg
cvg[[1]]


###################################################
### chunk number 23: olap-cvg-as-m
###################################################
#line 383 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
m <- matrix(unlist(lapply(cvg, lapply, as.vector)),
            ncol=length(cvg))
summary(rowSums(m))
summary(colSums(m))


###################################################
### chunk number 24: sessionInfo
###################################################
#line 403 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
packageDescription("Rsamtools")
sessionInfo()


###################################################
### chunk number 25: caffeine-kegg
###################################################
#line 416 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
library(KEGG.db)
kid <- revmap(KEGGPATHID2NAME)[["Caffeine metabolism"]]
egid <- KEGGPATHID2EXTID[[sprintf("hsa%s", kid)]]


###################################################
### chunk number 26: caffeine-txdb eval=FALSE
###################################################
## #line 426 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
## library(biomaRt)
## mart <- useMart("ensembl", "hsapiens_gene_ensembl")
## ensid <- getBM(c("ensembl_transcript_id"), filters="entrezgene",
##                values=egid, mart=mart)[[1]]
## library(GenomicFeatures)
## txdb <- makeTranscriptDbFromBiomart(transcript_ids=ensid)


###################################################
### chunk number 27: caffeine-txdb-load
###################################################
#line 434 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
library(GenomicFeatures)
fl <- system.file("extdata", "CaffeineTxdb.sqlite", 
                  package="Rsamtools")
txdb <- loadFeatures(fl)


###################################################
### chunk number 28: txdb-transcripts
###################################################
#line 443 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
bamRanges <- transcripts(txdb)


###################################################
### chunk number 29: bam-avail eval=FALSE
###################################################
## #line 459 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
## library(RCurl)
## ftpBase <- 
##     "ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/pilot_data/data/"
## indivDirs <- 
##     strsplit(getURL(ftpBase, ftplistonly=TRUE), "\n")[[1]]
## alnDirs <- 
##     paste(ftpBase, indivDirs, "/alignment/", sep="")
## urls0 <- 
##     strsplit(getURL(alnDirs, dirlistonly=TRUE), "\n")


###################################################
### chunk number 30: bam-index eval=FALSE
###################################################
## #line 474 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
## urls <- urls[sapply(urls0, length) != 0]
## fls0 <- unlist(unname(urls0))
## fls1 <- fls0[grepl("bai$", fls0)]
## fls <- fls1[sapply(strsplit(fls1, "\\."), length)==7]


###################################################
### chunk number 31: slxMaq09 eval=FALSE
###################################################
## #line 485 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
## urls1 <- 
##     Filter(function(x) length(x) != 0,
##            lapply(urls, 
##                   function(x) x[grepl("SLX.maq.*2009_08.*bai$", x)]))
## slxMaq09.bai <- 
##    mapply(paste, names(urls1), urls1, sep="", USE.NAMES=FALSE)
## slxMaq09 <- sub(".bai$", "", slxMaq09.bai) #$


###################################################
### chunk number 32: bamIndicies eval=FALSE
###################################################
## #line 498 "vignettes/Rsamtools/inst/doc/Rsamtools-Overview.Rnw"
## bamIndexDir <- tempfile()
## dir.create(bamIndexDir)
## idxFiles <- mapply(download.file, slxMaq09.bai, 
##                    file.path(bamIndexDir, basename(slxMaq09.bai)) ,
##                    MoreArgs=list(method="curl"))