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


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


###################################################
### chunk number 3: setwd eval=FALSE
###################################################
## setwd("c:/Documents and Settings/mtmorgan/Desktop/CourseData")


###################################################
### chunk number 4: setwd
###################################################
setwd("../")


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


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


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


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


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


###################################################
### chunk number 10: aln
###################################################
aln


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


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


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


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


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


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


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


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


###################################################
### chunk number 19: 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 20: sampleFunction
###################################################
samplingFilter <- function(sampleSize) {
    srFilter(function(x) {
        idx <- seq_len(length(x))
        idx %in% sample(idx, sampleSize)
    }, name="Martin's demo sampling filter")
}
sample100 <- samplingFilter(100)
caln[sample100(caln)]


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


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


###################################################
### chunk number 24: lowlevel-int
###################################################
int <- matrix(as.numeric(scan(gzfile(fl, open="rb"))),
              ncol=4*(1 + 2*36), byrow=TRUE)


###################################################
### chunk number 25: intensities-plot-early
###################################################
intplot <- function(intensities, cycle, ...) {
    m <- intensities[,4 + (cycle - 1) * 4 + 1:4]
    colnames(m) <- c("A", "C", "G", "T")
    splom(m, ...)
}
print(intplot(int, 5, pch=".", log="xy"))


###################################################
### chunk number 26: qa eval=FALSE
###################################################
## qa <- qa(sp)


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


###################################################
### chunk number 28: report eval=FALSE
###################################################
## rpt <- report(qa, dest=tempfile())


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


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


###################################################
### chunk number 31: 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 32: 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 33: freq-seqs
###################################################
df <- qa[["frequentSequences"]]
head(df[df$lane=="s_5_1_export.txt" & df$type=="read",1:2])


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


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


###################################################
### chunk number 36: 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 37: Rmpi eval=FALSE
###################################################
## library(Rmpi)
## mpi.spawn.Rslaves(nsl=8)
## qa <- qa(sp)
## mpi.close.Rslaves()


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


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


###################################################
### chunk number 40: 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 41: sessionInfo
###################################################
sessionInfo()