###################################################
### chunk number 1: setup
###################################################
options(width=60)
olocale=Sys.setlocale(locale="C")


###################################################
### chunk number 2: ShortRead
###################################################
library("ShortRead")


###################################################
### chunk number 3: ShortReadVersion
###################################################
library("ShortRead")
packageDescription("ShortRead")$Version


###################################################
### chunk number 4: setwd eval=FALSE
###################################################
## setwd("c:/Documents and Settings/Desktop/Course/labs")


###################################################
### chunk number 5: setwd
###################################################
oldwd = setwd("..")


###################################################
### chunk number 6: SolexaPath
###################################################
sp <- SolexaPath("extdata/ELAND/080828_HWI-EAS88_0003")


###################################################
### chunk number 7: sp-display
###################################################
sp
analysisPath(sp)


###################################################
### chunk number 8: filePattern
###################################################
list.files(analysisPath(sp)[[1]])


###################################################
### chunk number 9: pattern
###################################################
pattern <- "s_1_1_export_head.txt"
list.files(analysisPath(sp), pattern)


###################################################
### chunk number 10: readAligned
###################################################
aln <- readAligned(sp, pattern)


###################################################
### chunk number 11: aln
###################################################
aln


###################################################
### chunk number 12: aln-head
###################################################
head(sread(aln))
table(strand(aln), useNA="ifany")
sum(is.na(position(aln)))


###################################################
### chunk number 13: quality
###################################################
head(quality(aln))


###################################################
### chunk number 14: quality-scores
###################################################
alf <- alphabet(quality(aln))
m <- as(quality(aln), "matrix")
colMeans(m)


###################################################
### chunk number 15: alignQuality
###################################################
alignQuality(aln)
q <- quality(alignQuality(aln))
sum(q==0)
print(densityplot(q[q>1], plot.points=FALSE, 
                  xlab="Alignment quality"))


###################################################
### chunk number 16: alignData
###################################################
alignData(aln)
table(alignData(aln)$filtering)


###################################################
### chunk number 17: select
###################################################
filtIdx <- alignData(aln)$filtering=="Y"
alignedIdx <- !is.na(strand(aln))
aln[filtIdx & alignedIdx]


###################################################
### chunk number 18: filter
###################################################
filt1 <- alignDataFilter(expression(filtering=="Y"))
filt2 <- chromosomeFilter("chr[0-9XYM]+.fa")
filt <- compose(filt1, filt2)
caln <- aln[filt(aln)]
caln


###################################################
### chunk number 19: readAliged-filter eval=FALSE
###################################################
## readAligned(sp, pattern, filter=filt)


###################################################
### chunk number 20: ualignFilter
###################################################
ualignFilter <- srFilter(function(x) {
    ## create a numerical index of reads. Divide the index, position,
    ## and strand information between chromosomes. Select the index of
    ## a single read at each unique position and strand. Return the
    ## selected index as a logical vector with the same length as x
    oindex <- seq_len(length(x))
    index <- tapply(oindex, chromosome(x), c)
    pdup <- tapply(position(x), chromosome(x), duplicated)
    sdup <- tapply(strand(x), chromosome(x), duplicated)
    keep <- oindex  %in% unlist(mapply(function(i, p, s) {
        i[!(p & s)]
    }, index, pdup, sdup))
}, name="select only one read per position & strand ")
caln[ualignFilter(caln)]


###################################################
### chunk number 21: sampleFunction
###################################################
samplingFilter <- function(sampleSize) {
    srFilter(function(x) {
        idx <- seq_len(length(x))
        idx %in% sample(idx, sampleSize)
    }, name="Demo sampling filter")
}
sample100 <- samplingFilter(100)
caln[sample100(caln)]


###################################################
### chunk number 22: fastq
###################################################
ap <- analysisPath(sp)[[1]]
reads <- readFastq(ap, "s_1_1_sequence_head.txt$")
head(id(reads))


###################################################
### chunk number 23: readXStringColumns
###################################################
fl <- list.files(ap, pattern, full=TRUE)
cols <- strsplit(readLines(fl, 1), "\t")[[1]]
length(cols)
cols[9:10]

colClasses <- rep(list(NULL), 22)
colClasses[9:10] <- c("DNAString", "BString")
strings <- readXStringColumns(ap, pattern, colClasses=colClasses)
head(strings[[2]])


###################################################
### chunk number 24: lowlevel-int-columns
###################################################
fl <- list.files(imageAnalysisPath(sp),".*_int.*", full=TRUE)
strsplit(readLines(gzfile(fl, open="rb"), 1), "\t")[[1]][1:7]


###################################################
### chunk number 25: lowlevel-int
###################################################
int <- readIntensities(sp)
arr <- as(intensity(int), "array")


###################################################
### chunk number 26: intensities-plot-early
###################################################
print(splom(arr[,,5], pch="."))


###################################################
### chunk number 27: qa-input
###################################################
load(file.path("data", "qa_080828_081110.rda"))
qa


###################################################
### chunk number 28: qa-readCount
###################################################
qa[["readCounts"]]


###################################################
### chunk number 29: qa-template-source
###################################################
source(file.path("scripts", "qa_solexa.R"))
ppnCount(qa[["readCounts"]])


###################################################
### chunk number 30: hoover
###################################################
df <- qa[["sequenceDistribution"]]
df5raw <- df[df$lane=="s_5_1_export.txt" & df$type=="read",]
head(df5raw)
print(xyplot(log10(nReads)~log10(nOccurrences), df5raw,
             xlab="Copies per read (log 10)",
             ylab="Unique reads (log 10)"))


###################################################
### chunk number 31: choover
###################################################
csum <- with(df5raw, cumsum(nReads * nOccurrences))
csum <- csum / csum[length(csum)]
print(xyplot(csum ~log10(nOccurrences), df5raw,
             xlab="Copies per read (log 10)",
             ylab="Cumulative proportion of reads",
             type="l"))


###################################################
### chunk number 32: freq-seqs
###################################################
df <- qa[["frequentSequences"]]
head(df[df$lane=="s_5_1_export.txt" & df$type=="read",1:2])


###################################################
### chunk number 33: srdistance
###################################################
seq <- "CGGTTCAGCAGGAATGCCGAGATCGGAAGAGCGGT"
dist <- srdistance(clean(aln), seq)[[1]]
head(table(dist))


###################################################
### chunk number 34: alphabetFreqeuncy
###################################################
alphabetFrequency(sread(aln), collapse=TRUE, baseOnly=TRUE, freq=TRUE)


###################################################
### chunk number 35: abc
###################################################
abc <- alphabetByCycle(sread(aln))
dim(abc)
abc[1:4,1:5]
abc <- abc[rowSums(abc)!=0,]
df <- as.data.frame(t(abc[1:4,]))
print(xyplot(A+C+G+T~1:nrow(df), df, type="l",
             auto.key=list(x=.75, y=.95, points=FALSE, lines=TRUE),
             xlab="Cycle", ylab="Count"))


###################################################
### chunk number 36: Rmpi eval=FALSE
###################################################
## library("Rmpi")
## mpi.spawn.Rslaves(nsl=8)
## qa <- qa(sp)
## mpi.close.Rslaves()


###################################################
### chunk number 37: workflow-qa eval=FALSE
###################################################
## sp <- SolexaPath("extdata/ELAND/080828_HWI-EAS88_0003")
## rpt <- report(qa(sp), dest="reports/my_report.pdf")


###################################################
### chunk number 38: filter eval=FALSE
###################################################
## filt1 <- nFilter()
## filt2 <- alignDataFilter(expression(filtering=="Y"))
## filt3 <- alignQualityFilter(threshold=1)
## filt4 <- srdistanceFilter("CGGTTCAGCAGGAATGCCGAGATCGGAAGAGCGGT", 4)
## filt5 <- chromosomeFilter("chr[0-9]+.fa")
## 
## filt <- compose(filt1, filt2, filt3, filt4, filt5)
## aln <- readAligned(sp, "s_1_1_export.txt$", filter=filt)


###################################################
### chunk number 39: filter-MAQ eval=FALSE
###################################################
## maqDir <- file.path("extdata", "MAQ")
## filt5 <- chromosomeFilter("chr[0-9XY]+$")
## filt <- compose(filt1, filt3, filt4, filt5)
## maq <- readAligned(maqDir, "s_8.map", "MAQMap", filter=filt)


###################################################
### chunk number 40: sessionInfo
###################################################
toLatex(sessionInfo())