### R code from vignette source 'RBiocForEveryone-lab.Rnw'

###################################################
### code chunk number 1: rbioc-pdata (eval = FALSE)
###################################################
## pdataFile <- "/home/martin/RBiocForEveryone/pData.csv"


###################################################
### code chunk number 2: rbioc-pdata-build
###################################################
## used when building the vignette
pdataFile <- "pData.csv"


###################################################
### code chunk number 3: rbioc-pdata-csv
###################################################
pdata <- read.table(pdataFile)  
dim(pdata)
names(pdata)


###################################################
### code chunk number 4: rbioc-pdata-subset
###################################################
head(pdata$sex) # same as pdata[,"sex"], pdata[["sex"]]
sapply(pdata, class)


###################################################
### code chunk number 5: rbioc-pdata-sextab
###################################################
table(pdata$sex, useNA="ifany")


###################################################
### code chunk number 6: rbioc-pdata-molbiol
###################################################
with(pdata, table(mol.biol, useNA="ifany"))


###################################################
### code chunk number 7: rbbioc-pdata-bcrabl
###################################################
ridx <- pdata$mol.biol %in% c("BCR/ABL", "NEG")


###################################################
### code chunk number 8: rbioc-pdata-molbiol-selected
###################################################
table(ridx)
sum(ridx)


###################################################
### code chunk number 9: rbioc-pdata-subset
###################################################
pdata1 <- pdata[ridx,]


###################################################
### code chunk number 10: rbioc-pdata-subset-levels
###################################################
levels(pdata$mol.biol)


###################################################
### code chunk number 11: rbioc-pdata-subset-recode
###################################################
pdata1$mol.biol <- factor(pdata1$mol.biol)
table(pdata1$mol.biol)


###################################################
### code chunk number 12: rbioc-pdata-age-molbiol
###################################################
with(pdata1, t.test(age ~ mol.biol))


###################################################
### code chunk number 13: rbioc-pdata-boxplot (eval = FALSE)
###################################################
## ## not evaluated
## boxplot(age ~ mol.biol, pdata1)


###################################################
### code chunk number 14: setup
###################################################
library(pasillaBamSubset)
library(Rsamtools)  # scanBam
library(ShortRead)  # alphabetByCycle


###################################################
### code chunk number 15: param
###################################################
flag <- scanBamFlag(isMinusStrand=FALSE)
param <- ScanBamParam(what=c("seq", "qual"), flag=flag)


###################################################
### code chunk number 16: query
###################################################
fl <- untreated1_chr4()
res <- scanBam(fl, param=param)[[1]]


###################################################
### code chunk number 17: abcplot
###################################################
abc <- alphabetByCycle(res$seq)
matplot(t(abc[1:4,]), type="l", lty=1, lwd=2, 
        xlab="Cycle", ylab="Count")
legend("topright", legend=rownames(abc)[1:4], 
       lty=1, lwd=2, col=1:4)


###################################################
### code chunk number 18: qualplot
###################################################
qual <- as(res$qual, "matrix")
boxplot(qual ~ col(qual), outline=FALSE,
        xlab="Cycle", ylab="Quality")


###################################################
### code chunk number 19: annotations
###################################################
library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)  # genome coordinates
exByGn <- exonsBy(TxDb.Dmelanogaster.UCSC.dm3.ensGene, "gene")
chr4 <- exByGn[ all(seqnames(exByGn) == "chr4") ]


###################################################
### code chunk number 20: files
###################################################
fls <- c(untreated1_chr4(), untreated3_chr4())
names(fls) <- sub("_chr4.bam", "", basename(fls))
bfl <- BamFileList(fls)


###################################################
### code chunk number 21: counts
###################################################
counts <- summarizeOverlaps(chr4, bfl, ignore.strand=TRUE)
head(assay(counts))


###################################################
### code chunk number 22: countsplot
###################################################
plot(asinh(assay(counts)), asp=1, main="asinh(counts), chr4")
abline(0, 1, lty=2)


###################################################
### code chunk number 23: gff (eval = FALSE)
###################################################
## library(rtracklayer)  # import gff
## fl <- paste0("ftp://ftp.ensembl.org/pub/release-62/",
##              "gtf/drosophila_melanogaster/",
##              "Drosophila_melanogaster.BDGP5.25.62.gtf.gz")
## gffFile <- file.path(tempdir(), basename(fl))
## download.file(fl, gffFile)
## gff0 <- import(gffFile, asRangedData=FALSE)
## idx <- values(gff0)$source == "protein_coding" & 
##            values(gff0)$type == "exon" &
##                seqnames(gff0) == "4"
## gff1 <- gff0[idx]
## chr4 <- split(gff1, values(gff1)$gene_id)