###################################################
### chunk number 1: 
###################################################
options(width=70)


###################################################
### chunk number 2: initialize
###################################################
require(Genominator)
require(ShortRead)
require(yeastRNASeq)


###################################################
### chunk number 3: datalist
###################################################
data(package = "yeastRNASeq")[[3]][, c("Item", "Title")]


###################################################
### chunk number 4: yeastAligned
###################################################
data("yeastAligned")
names(yeastAligned)
yeastAligned[[1]]


###################################################
### chunk number 5: eData
###################################################
eData <- importFromAlignedReads(yeastAligned, chrMap = levels(chromosome(yeastAligned[[1]])),
                                filename = "my.db", tablename = "raw", overwrite = TRUE)
head(eData)


###################################################
### chunk number 6: chrMap
###################################################
levels(chromosome(yeastAligned[[1]]))


###################################################
### chunk number 7: loadingYeastAnno
###################################################
data(yeastAnno)
head(yeastAnno)


###################################################
### chunk number 8: yeastAnnoFix
###################################################
yeastAnno$chr <- match(yeastAnno$chr, c(as.character(as.roman(1:16)), "MT", "2-micron"))
yeastAnno$start <- yeastAnno$start_position
yeastAnno$end <- yeastAnno$end_position
rownames(yeastAnno) <- yeastAnno$ensembl
yeastAnno <- yeastAnno[, c("chr", "start", "end", "strand", "gene_biotype")]
head(yeastAnno)


###################################################
### chunk number 9: geneCounts
###################################################
geneCounts <- summarizeByAnnotation(eData, yeastAnno, ignoreStrand = TRUE)
head(geneCounts)


###################################################
### chunk number 10: gof
###################################################
plot(regionGoodnessOfFit(geneCounts, groups = gsub("_[0-9]_f", "", colnames(geneCounts))),
     chisq = TRUE)