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

###################################################
### code chunk number 1: setup
###################################################
library(GenomicFeatures)
plotRanges <- function(x, xlim = x, main = deparse(substitute(x)),
                       col = "black", sep = 0.5, ...)
{
    height <- 1
    if (is(xlim, "Ranges"))
        xlim <- c(min(start(xlim)), max(end(xlim)))
    bins <- disjointBins(IRanges(start(x), end(x) + 1))
    plot.new()
    par(mai=c(0.6, 0.2, 0.6, 0.2))
    plot.window(xlim, c(0, max(bins)*(height + sep)))
    ybottom <- bins * (sep + height) - height
    rect(start(x)-0.5, ybottom, end(x)+0.5, ybottom + height, col = col, ...)
    title(main, cex.main=2.8, font.main=1)
    axis(1)
}


###################################################
### code chunk number 2: GRanges-genes
###################################################
genes <- GRanges(seqnames=c("3R", "X"),
                 ranges=IRanges(
                     start=c(19967117, 18962306),
                     end=c(19973212, 18962925)),
                 strand=c("+", "-"),
                 seqlengths=c(`3R`=27905053L, `X`=22422827L))


###################################################
### code chunk number 3: GRanges-display
###################################################
genes


###################################################
### code chunk number 4: GRanges-help (eval = FALSE)
###################################################
## ?GRanges


###################################################
### code chunk number 5: GRanges-vignettes (eval = FALSE)
###################################################
## vignette(package="GenomicRanges")


###################################################
### code chunk number 6: GRanges-basics
###################################################
genes[2]
strand(genes)
width(genes)
length(genes)
names(genes) <- c("FBgn0039155", "FBgn0085359")
genes  # now with names


###################################################
### code chunk number 7: ranges-ir
###################################################
ir <- IRanges(start=c(7, 9, 12, 14, 22:24),
              end=c(15, 11, 12, 18, 26, 27, 28))


###################################################
### code chunk number 8: ranges-ir-plot
###################################################
png("ranges-ir-plot.png", width=800, height=160)
plotRanges(ir, xlim=c(5, 35), main="Original")
dev.off()
png("ranges-shift-ir-plot.png", width=800, height=160)
plotRanges(shift(ir, 5), xlim=c(5, 35), main="Shift")
dev.off()
png("ranges-reduce-ir-plot.png", width=800, height=160)
plotRanges(reduce(ir), xlim=c(5, 35), main="Reduce")
dev.off()


###################################################
### code chunk number 9: ranges-shift-ir
###################################################
shift(ir, 5)


###################################################
### code chunk number 10: ranges-shift-ir-plot
###################################################
png("ranges-shift-ir-plot.png", width=800, height=160)
plotRanges(shift(ir, 5), xlim=c(5, 35), main="Shift")
dev.off()


###################################################
### code chunk number 11: ranges-reduce-ir
###################################################
reduce(ir)


###################################################
### code chunk number 12: ranges-reduce-ir-plot
###################################################
png("ranges-reduce-ir-plot.png", width=800, height=160)
plotRanges(reduce(ir), xlim=c(5, 35), main="Reduce")
dev.off()


###################################################
### code chunk number 13: coverage
###################################################
coverage(ir)


###################################################
### code chunk number 14: GRanges-mcols
###################################################
mcols(genes) <- DataFrame(EntrezId=c("42865", "2768869"),
                          Symbol=c("kal-1", "CG34330"))


###################################################
### code chunk number 15: GRanges-metadata
###################################################
metadata(genes) <-
    list(CreatedBy="A. User", Date=date())


###################################################
### code chunk number 16: GRangesList-eg-setup
###################################################
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene # alias
fbgn <- exonsBy(txdb, "gene")["FBgn0039155"]
seqlevels(fbgn) <- "chr3R"


###################################################
### code chunk number 17: GRangesList-eg
###################################################
fbgn


###################################################
### code chunk number 18: txdb
###################################################
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene # alias
ex0 <- exonsBy(txdb, "gene")
head(table(elementLengths(ex0)))
ids <- c("FBgn0002183", "FBgn0003360", "FBgn0025111", "FBgn0036449")
ex <- reduce(ex0[ids])


###################################################
### code chunk number 19: SAM
###################################################
fl <- system.file("extdata", "ex1.sam", package="Rsamtools")
strsplit(readLines(fl, 1), "\t")[[1]]


###################################################
### code chunk number 20: readGappedAlignments
###################################################
alnFile <- system.file("extdata", "ex1.bam", package="Rsamtools")
aln <- readGappedAlignments(alnFile)
head(aln, 3)


###################################################
### code chunk number 21: GappedAlignments-accessors
###################################################
table(strand(aln))
table(width(aln))
head(sort(table(cigar(aln)), decreasing=TRUE))


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


###################################################
### code chunk number 23: bam-ex-input
###################################################
## input
aln <- readGappedAlignments(fls[1])
xtabs( ~ seqnames + strand, as.data.frame(aln))


###################################################
### code chunk number 24: bam-ex-roi
###################################################
data(ex)             # from a later exercise


###################################################
### code chunk number 25: bam-ex-strand
###################################################
strand(aln) <- "*"   # protocol not strand-aware


###################################################
### code chunk number 26: bam-ex-hits
###################################################
hits <- countOverlaps(aln, ex)
table(hits)


###################################################
### code chunk number 27: bam-ex-cnt
###################################################
cnt <- countOverlaps(ex, aln[hits==1])


###################################################
### code chunk number 28: bam-count-fun
###################################################
counter <- 
    function(filePath, range)
{
    aln <- readGappedAlignments(filePath)
    strand(aln) <- "*"
    hits <- countOverlaps(aln, range)
    countOverlaps(range, aln[hits==1])
}


###################################################
### code chunk number 29: bam-count-all
###################################################
counts <- sapply(fls, counter, ex)


###################################################
### code chunk number 30: bam-count-mclapply (eval = FALSE)
###################################################
## if (require(parallel))
##     simplify2array(mclapply(fls, counter, ex))