## ----style, eval=TRUE, echo=FALSE, results='asis'------------------------
options(max.print=1000)
stopifnot(BiocInstaller::biocVersion() == "2.13")
BiocStyle::latex()
library(knitr)
opts_chunk$set(cache=TRUE, tidy=FALSE)

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.5, 0.2, 1.2, 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)
}


## ----echo=FALSE----------------------------------------------------------
suppressPackageStartupMessages(library(GenomicRanges))
gr1 <- GRanges(Rle(c("chr1", "chr19", "chrX", "chrM", "chr19"),
                   c(2, 1, 2, 1, 1)),
               IRanges(c(10001, 5508, 25644, 685, 142501, 33077, 8001),
                       width=c(400, 285, 1620, 333, 1, 2000, 421)),
               strand=c("-", "*", "+", "+", "-", "*", "+"))


## ------------------------------------------------------------------------
gr1


## ----echo=FALSE----------------------------------------------------------
suppressPackageStartupMessages(library(BSgenome.Scerevisiae.UCSC.sacCer2))
tiles <- tileGenome(seqinfo(Scerevisiae), tilewidth=50000,
                    cut.last.tile.in.chrom=TRUE)


## ------------------------------------------------------------------------
tiles


## ----echo=FALSE----------------------------------------------------------
gr2 <- gr1
mcols(gr2) <- DataFrame(desc=c("exon",
                               "repeat region",
                               "gene",
                               "exon",
                               "exon",
                               "snp",
                               "repeat region"),
                        GC_content=c(0.55, 0.495, 0.562, 0.53, 0.552, NA, 0.51))


## ------------------------------------------------------------------------
gr2


## ----echo=FALSE----------------------------------------------------------
suppressPackageStartupMessages(library(TxDb.Hsapiens.UCSC.hg19.knownGene))
hg19_cds <- cds(TxDb.Hsapiens.UCSC.hg19.knownGene,
                columns="gene_id")


## ------------------------------------------------------------------------
hg19_cds


## ----echo=FALSE----------------------------------------------------------
while (search()[2L] != "package:parallel")
    detach(search()[2L], character.only=TRUE)
detach("package:parallel")


## ----results='asis'------------------------------------------------------
library(GenomicRanges)


## ----eval=FALSE----------------------------------------------------------
## class?GRanges


## ----results='hide'------------------------------------------------------
example(GRanges)


## ------------------------------------------------------------------------
gr


## ------------------------------------------------------------------------
seqnames(gr)


## ------------------------------------------------------------------------
class(seqnames(gr))


## ------------------------------------------------------------------------
width(gr)


## ------------------------------------------------------------------------
all(width(gr) == end(gr) - start(gr) + 1)


## ------------------------------------------------------------------------
names(gr)


## ------------------------------------------------------------------------
names(gr) <- toupper(names(gr))


## ------------------------------------------------------------------------
mcols(gr)


## ------------------------------------------------------------------------
mcols(gr)$GC


## ------------------------------------------------------------------------
gr[mcols(gr)$GC >= 0.3 & mcols(gr)$GC <= 0.6]


## ------------------------------------------------------------------------
subset(gr, GC >= 0.3 & GC <= 0.6)


## ------------------------------------------------------------------------
mcols(gr)$id <- sprintf("ID%03d", seq_along(gr))
gr


## ----ranges-gr0-plot,results='hide',echo=FALSE,fig.keep='none'-----------
gr0 <- IRanges(start=c(7, 9, 12, 14, 22:24),
               end=c(15, 11, 12, 18, 26, 27, 28))

png("ranges-gr0-plot.png", width=800, height=170)
plotRanges(gr0, xlim=c(5, 35), main="gr0", col="blue")
dev.off()


## ----ranges-shift-gr0-plot,results='hide',echo=FALSE,fig.keep='none'-----
png("ranges-shift-gr0-plot.png", width=800, height=170)
plotRanges(shift(gr0, 5), xlim=c(5, 35), main="shift(gr0, 5)", col="blue")
dev.off()


## ----ranges-disjoin-gr0-plot,results='hide',echo=FALSE,fig.keep='none'----
png("ranges-disjoin-gr0-plot.png", width=800, height=170)
plotRanges(disjoin(gr0), xlim=c(5, 35), main="disjoin(gr0)", col="blue")
dev.off()


## ----ranges-reduce-gr0-plot,results='hide',echo=FALSE,fig.keep='none'----
png("ranges-reduce-gr0-plot.png", width=800, height=170)
plotRanges(reduce(gr0), xlim=c(5, 35), main="reduce(gr0)", col="blue")
dev.off()


## ------------------------------------------------------------------------
gr2 <- shift(gr, 500)


## ------------------------------------------------------------------------
flank(gr2, width=100)


## ------------------------------------------------------------------------
shift


## ------------------------------------------------------------------------
showMethods(shift)


## ------------------------------------------------------------------------
selectMethod(shift, "GRanges")


## ------------------------------------------------------------------------
is(gr, "GenomicRanges")
class(gr)


## ----eval=FALSE----------------------------------------------------------
## ?`shift,GenomicRanges-method`


## ------------------------------------------------------------------------
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
TxDb.Hsapiens.UCSC.hg19.knownGene


## ------------------------------------------------------------------------
ex_by_gn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, by="gene")
class(ex_by_gn)


## ------------------------------------------------------------------------
length(ex_by_gn)


## ------------------------------------------------------------------------
"1008" %in% names(ex_by_gn)


## ------------------------------------------------------------------------
ex_by_gn[["1008"]]


## ------------------------------------------------------------------------
seqnames(ex_by_gn[["1008"]])
strand(ex_by_gn[["1008"]])


## ------------------------------------------------------------------------
nex_per_gn <- elementLengths(ex_by_gn)


## ------------------------------------------------------------------------
max(nex_per_gn)


## ------------------------------------------------------------------------
which(nex_per_gn == max(nex_per_gn))


## ------------------------------------------------------------------------
which.max(nex_per_gn)


## ------------------------------------------------------------------------
library(parathyroidSE)
bamdir <- system.file("extdata", package="parathyroidSE")
bampaths <- list.files(bamdir, pattern="bam$", full.names=TRUE)
bampaths


## ------------------------------------------------------------------------
library(Rsamtools)
reads <- readGAlignmentsFromBam(bampaths[1])
reads


## ------------------------------------------------------------------------
library(parathyroidSE)
bamdir <- system.file("extdata", package="parathyroidSE")
bampaths <- list.files(bamdir, pattern="bam$", full.names=TRUE)
reads <- readGAlignmentsFromBam(bampaths[3])
class(reads)
reads


## ------------------------------------------------------------------------
length(grep("N", cigar(reads)))


## ------------------------------------------------------------------------
max(ngap(reads))
table(ngap(reads))


## ------------------------------------------------------------------------
read_ranges <- as(reads, "GRangesList")


## ------------------------------------------------------------------------
table(elementLengths(read_ranges))


## ------------------------------------------------------------------------
as(reads, "GRanges")


## ------------------------------------------------------------------------
cvg <- coverage(read_ranges)
class(cvg)
cvg


## ------------------------------------------------------------------------
cvg[["12"]]
max(cvg[["12"]])


## ------------------------------------------------------------------------
max(cvg)