## ----style, echo = FALSE, results = 'asis'--------------------------------------------------------
options(width=100)
knitr::opts_chunk$set(
    eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")),
    cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE")))

## ----setup, echo=FALSE----------------------------------------------------------------------------
suppressPackageStartupMessages({
    library(Biostrings)
    library(GenomicRanges)
    library(SummarizedExperiment)
    library(airway)
    library(rtracklayer)
    library(ShortRead)
    library(GenomicAlignments)
    library(RNAseqData.HNRNPC.bam.chr14)
    library(VariantAnnotation)
})

## ----install, eval=FALSE--------------------------------------------------------------------------
#  source("https://bioconductor.org/biocLite.R")
#  biocLite(c("DESeq2", "org.Hs.eg.db"))

## ----require--------------------------------------------------------------------------------------
library(GenomicRanges)

## ----help-bioc, eval=FALSE------------------------------------------------------------------------
#  help(package="GenomicRanges")
#  vignette(package="GenomicRanges")
#  vignette(package="GenomicRanges", "GenomicRangesHOWTOs")
#  ?GRanges

## ----five-lines-----------------------------------------------------------------------------------
x <- rnorm(1000)
y <- x + rnorm(1000)
df <- data.frame(X=x, Y=y)
plot(Y ~ X, df)
fit <- lm(Y ~ X, df)
anova(fit)
abline(fit)

## ----help-r, eval=FALSE---------------------------------------------------------------------------
#  class(fit)
#  methods(class=class(fit))
#  methods(plot)
#  ?"plot"
#  ?"plot.formula"

## ----classes-and-methods--------------------------------------------------------------------------
library(Biostrings)
dna <- DNAStringSet(c("AACAT", "GGCGCCT"))
reverseComplement(dna)

## ----classes-and-methods-discovery, eval=FALSE----------------------------------------------------
#  class(dna)
#  ?"DNAStringSet-class"
#  ?"reverseComplement,DNAStringSet-method"

## -------------------------------------------------------------------------------------------------
library(Biostrings)
data(phiX174Phage)
phiX174Phage
letterFrequency(phiX174Phage, c("A", "C", "G", "T"))
letterFrequency(phiX174Phage, "GC", as.prob=TRUE)

## ----ranges, message=FALSE------------------------------------------------------------------------
library(GenomicRanges)
gr <- GRanges("A", IRanges(c(10, 20, 22), width=5), "+")
shift(gr, 1)                            # intra-range
range(gr)                               # inter-range
reduce(gr)                              # inter-range
snps <- GRanges("A", IRanges(c(11, 17, 24), width=1))
findOverlaps(snps, gr)                  # between-range
setdiff(range(gr), gr)                  # 'introns'

## -------------------------------------------------------------------------------------------------
library(airway)
data(airway)
airway

## -------------------------------------------------------------------------------------------------
x <- assay(airway)
class(x)
dim(x)
head(x)
colData(airway)
rowRanges(airway)

## -------------------------------------------------------------------------------------------------
cidx <- colData(airway)$dex %in% "trt"
airway[, cidx]
## shortcut
airway[, airway$dex %in% "trt"]

## -------------------------------------------------------------------------------------------------
chr14 <- as(seqinfo(rowRanges(airway)), "GRanges")["14"]
ridx <- rowRanges(airway) %over% chr14
airway[ridx,]
## shortcut
chr14 <- as(seqinfo(airway), "GRanges")["14"]
airway[airway %over% chr14,]

## -------------------------------------------------------------------------------------------------
library(RNAseqData.HNRNPC.bam.chr14)
fls = RNAseqData.HNRNPC.bam.chr14_BAMFILES
basename(fls)

## -------------------------------------------------------------------------------------------------
library(GenomicAlignments)
bfls = BamFileList(fls)
bfl = bfls[[1]]

## -------------------------------------------------------------------------------------------------
ga = readGAlignments(bfl)
ga
table(strand(ga))

## -------------------------------------------------------------------------------------------------
tail(sort(table(cigar(ga))))
ga[cigar(ga) != "72M"]

## -------------------------------------------------------------------------------------------------
summarizeJunctions(ga)
junctions <- summarizeJunctions(ga, with.revmap=TRUE)
ga[ junctions$revmap[[1]] ]

## -------------------------------------------------------------------------------------------------
coverage(bfl)$chr14