### R code from vignette source 'A2_ReadsAndStrings.Rnw'

###################################################
### code chunk number 1: gc-genome
###################################################
library(BSgenome.Dmelanogaster.UCSC.dm3)
Dmelanogaster
library(Morgan2013)
data(ex)
ex[1]
nm <- "chr3L"
chr <- Dmelanogaster[[nm]]
v <- Views(chr, start=start(ex[[1]]), end=end(ex[[1]]))


###################################################
### code chunk number 2: gcFunction-definition
###################################################
gcFunction


###################################################
### code chunk number 3: gcFunction-genome (eval = FALSE)
###################################################
## subjectGC <- gcFunction(v)


###################################################
### code chunk number 4: fastq-format
###################################################
bigdata <- "~/BigData"
fls <- dir(file.path(bigdata, "fastq"), full=TRUE)
cat(noquote(readLines(fls[[1]], 4)), sep="\n")


###################################################
### code chunk number 5: ascii
###################################################
cat(rawToChar(as.raw(32+1:75)),
    rawToChar(as.raw(32+76:94)), sep="\n")


###################################################
### code chunk number 6: readFastq
###################################################
library(ShortRead)
bigdata <- "~/BigData"  ## e.g., "C:/Users/mtmorgan/"
fastqDir <- file.path(bigdata, "fastq")
fastqFiles <- dir(fastqDir, full=TRUE)
fq <- readFastq(fastqFiles[1])
fq


###################################################
### code chunk number 7: sread
###################################################
head(sread(fq), 3)
head(quality(fq), 3)
head(id(fq), 3)


###################################################
### code chunk number 8: getClasse
###################################################
getClass("ShortReadQ")


###################################################
### code chunk number 9: showMethods-ShortRead (eval = FALSE)
###################################################
## showMethods(class="ShortRead", where=search())


###################################################
### code chunk number 10: width-ShortRead
###################################################
table(width(fq))


###################################################
### code chunk number 11: width-ShortReadQ
###################################################
abc <- alphabetByCycle(sread(fq))
abc[1:4, 1:8]


###################################################
### code chunk number 12: FastqSampler
###################################################
sampler <- FastqSampler(fastqFiles[1], 1000000)
yield(sampler) # sample of 1000000 reads


###################################################
### code chunk number 13: qa (eval = FALSE)
###################################################
## qas0 <- Map(function(fl, nm) {
##     fq <- FastqSampler(fl)
##     qa(yield(fq), nm)
## }, fastqFiles,
##    sub("_subset.fastq", "", basename(fastqFiles)))
## qas <- do.call(rbind, qas0)
## rpt <- report(qas, dest=tempfile())
## browseURL(rpt)


###################################################
### code chunk number 14: report (eval = FALSE)
###################################################
## rpt <- system.file("GSM461176_81_qa_report", "index.html",
##                    package="Morgan2013")
## browseURL(rpt)


###################################################
### code chunk number 15: fastq-discovery
###################################################
dir(bigdata)
fls <- dir(file.path(bigdata, "fastq"), full=TRUE)


###################################################
### code chunk number 16: fastq-input-gc
###################################################
rm(fq); invisible(gc())


###################################################
### code chunk number 17: fastq-input
###################################################
fq <- readFastq(fls[1])


###################################################
### code chunk number 18: gcC
###################################################
alf0 <- alphabetFrequency(sread(fq), as.prob=TRUE, collapse=TRUE)
sum(alf0[c("G", "C")])


###################################################
### code chunk number 19: gc-reads
###################################################
gc <- gcFunction(sread(fq))
hist(gc)


###################################################
### code chunk number 20: abc
###################################################
abc <- alphabetByCycle(sread(fq))
matplot(t(abc[c("A", "C", "G", "T"),]), type="l")


###################################################
### code chunk number 21: abc-mclapply (eval = FALSE)
###################################################
## library(parallel)
## gc0 <- mclapply(fls, function(fl) {
##   fq <- readFastq(fl)
##   gc <- gcFunction(sread(fq))
##   table(cut(gc, seq(0, 1, .05)))
## })
## ## simplify list of length 2 to 2-D array
## gc <- simplify2array(gc0)
## matplot(gc, type="s")


###################################################
### code chunk number 22: quality
###################################################
head(quality(fq))
qual <- as(quality(fq), "matrix")
dim(qual)
plot(colMeans(qual), type="b")