## ----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)


## ----results='hide'------------------------------------------------------
library(parathyroidSE)


## ----eval=FALSE----------------------------------------------------------
## data(package="parathyroidSE")


## ------------------------------------------------------------------------
data(parathyroidGenesSE)
parathyroidGenesSE
class(parathyroidGenesSE)
dim(parathyroidGenesSE)


## ------------------------------------------------------------------------
class(assay(parathyroidGenesSE))
dim(assay(parathyroidGenesSE))
assay(parathyroidGenesSE)[1:8, 1:8]


## ------------------------------------------------------------------------
colnames(parathyroidGenesSE)


## ------------------------------------------------------------------------
colSums(assay(parathyroidGenesSE))


## ------------------------------------------------------------------------
sum(rowSums(assay(parathyroidGenesSE)) != 0)


## ------------------------------------------------------------------------
rowData(parathyroidGenesSE)


## ------------------------------------------------------------------------
colData(parathyroidGenesSE)


## ------------------------------------------------------------------------
table(colData(parathyroidGenesSE)$treatment)


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


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


## ------------------------------------------------------------------------
galp0 <- readGAlignmentPairsFromBam(bampaths[1])
galp0


## ------------------------------------------------------------------------
length(galp0)


## ------------------------------------------------------------------------
first(galp0)
last(galp0)


## ------------------------------------------------------------------------
param <- ScanBamParam(flag=scanBamFlag(isDuplicate=FALSE,
                                       isNotPrimaryRead=FALSE))
readGAlignmentPairsFromBam(bampaths[1], param=param)


## ----makeTranscriptDbFromBiomart, eval=FALSE-----------------------------
## ## Requires INTERNET ACCESS and takes about 6 min. Please don't try to run this!
## library(GenomicFeatures)
## txdb <- makeTranscriptDbFromBiomart(biomart="ensembl",
##                                     dataset="hsapiens_gene_ensembl")


## ----exonsBy, eval=FALSE-------------------------------------------------
## ex_by_gene <- exonsBy(txdb, by="gene")  # GRangesList object


## ------------------------------------------------------------------------
data(exonsByGene)
exonsByGene


## ------------------------------------------------------------------------
length(exonsByGene)


## ----BamFileList---------------------------------------------------------
library(Rsamtools)
bamfile_list <- BamFileList(bampaths, index=character())


## ------------------------------------------------------------------------
read_count0 <- summarizeOverlaps(exonsByGene, bamfile_list,
                                 ignore.strand=TRUE,
                                 singleEnd=FALSE)
read_count0


## ------------------------------------------------------------------------
read_count <- summarizeOverlaps(exonsByGene, bamfile_list,
                                ignore.strand=TRUE,
                                singleEnd=FALSE,
                                param=param)
read_count


## ------------------------------------------------------------------------
colSums(assay(read_count0))
colSums(assay(read_count))