################################################### ### chunk number 1: setup ################################################### options(width=60) olocale=Sys.setlocale(locale="C") ################################################### ### chunk number 2: ShortRead ################################################### suppressMessages(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()