###################################################
### chunk number 1: biocLite eval=FALSE
###################################################
## source("http://bioconductor.org/biocLite.R")
## biocLite("GenomicRanges")


###################################################
### chunk number 2: initialize
###################################################
library(GenomicRanges)


###################################################
### chunk number 3: YeastData
###################################################
library(Rsamtools)
testFile <- system.file("bam", "isowt5_13e.bam", package="leeBamViews")
aligns <- readBamGappedAlignments(testFile)


###################################################
### chunk number 4: GetAnnoations
###################################################
library(GenomicFeatures)
txdb <- makeTranscriptDbFromUCSC(genome="sacCer2", tablename="sgdGene")


###################################################
### chunk number 5: exonsBy
###################################################
exonRanges <- exonsBy(txdb, "tx")


###################################################
### chunk number 6: ammendData
###################################################
rname(aligns) <- sub("^Sc", "", rname(aligns))
rname(aligns) <- sub("13", "XIII", rname(aligns))


###################################################
### chunk number 7: count
###################################################
counts <- countOverlaps(exonRanges, aligns)


###################################################
### chunk number 8: numBases
###################################################
numBases <- sum(width(exonRanges))
geneLengthsInKB <- (numBases/1000)


###################################################
### chunk number 9: millionsMapped
###################################################
millionsMapped <- sum(counts)/1000000


###################################################
### chunk number 10: RPM
###################################################
# counted reads / total reads in millions
rpm <- counts/millionsMapped


###################################################
### chunk number 11: RPKM
###################################################
# reads per million per geneLength in Kb
rpkm <- rpm/geneLengthsInKB


###################################################
### chunk number 12: SessionInfo
###################################################
sessionInfo()