################################################### ### chunk number 1: setup ################################################### options(width=60) olocale=Sys.setlocale(locale="C") library(ShortReadTutorial) ################################################### ### chunk number 2: ShortReadVersion ################################################### library(ShortReadTutorial) packageDescription("ShortRead")$Version ################################################### ### chunk number 3: setwd eval=FALSE ################################################### ## setwd("c:/Documents and Settings/Desktop/Course/labs") ################################################### ### chunk number 4: extdataDir ################################################### extdataDir <- system.file("extdata", package="ShortReadTutorial") file.exists(extdataDir) ################################################### ### chunk number 5: dirsExist ################################################### stopifnot(file.exists(extdataDir)) ################################################### ### chunk number 6: filePattern ################################################### list.files(extdataDir) ################################################### ### chunk number 7: pattern ################################################### pattern <- "s_1_1_export_head.txt" list.files(extdataDir, pattern) ################################################### ### chunk number 8: readAligned ################################################### aln <- readAligned(extdataDir, pattern, "SolexaExport") ################################################### ### chunk number 9: aln ################################################### aln ################################################### ### chunk number 10: aln-head ################################################### head(sread(aln)) table(strand(aln), useNA="ifany") sum(is.na(position(aln))) ################################################### ### chunk number 11: quality ################################################### head(quality(aln)) ################################################### ### chunk number 12: quality-scores ################################################### alf <- alphabet(quality(aln)) m <- as(quality(aln), "matrix") colMeans(m) ################################################### ### chunk number 13: alignQuality ################################################### alignQuality(aln) q <- quality(alignQuality(aln)) sum(q==0) print(densityplot(q[q>1], plot.points=FALSE, xlab="Alignment quality")) ################################################### ### chunk number 14: alignData ################################################### alignData(aln) table(alignData(aln)$filtering) ################################################### ### chunk number 15: select ################################################### filtIdx <- alignData(aln)$filtering=="Y" alignedIdx <- !is.na(strand(aln)) aln[filtIdx & alignedIdx] ################################################### ### chunk number 16: filter ################################################### filt1 <- alignDataFilter(expression(filtering=="Y")) filt2 <- chromosomeFilter("chr[0-9XYM]+.fa") filt <- compose(filt1, filt2) caln <- aln[filt(aln)] caln ################################################### ### chunk number 17: readAliged-filter eval=FALSE ################################################### ## readAligned(extdatDir, pattern, filter=filt) ################################################### ### chunk number 18: 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 19: 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 20: fastq ################################################### reads <- readFastq(extdataDir, "s_1_1_sequence_head.txt$") head(id(reads)) ################################################### ### chunk number 21: readXStringColumns ################################################### fl <- list.files(extdataDir, 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(extdataDir, pattern, colClasses=colClasses) head(strings[[2]]) ################################################### ### chunk number 22: lowlevel-int-columns ################################################### fl <- list.files(extdataDir, ".*_int.*", full=TRUE) strsplit(readLines(gzfile(fl, open="rb"), 1), "\t")[[1]][1:7] ################################################### ### chunk number 23: lowlevel-int ################################################### int <- readIntensities(extdataDir, type="SolexaIntensity") arr <- as(intensity(int), "array") ################################################### ### chunk number 24: intensities-plot-early ################################################### print(splom(arr[,,5], pch=".")) ################################################### ### chunk number 25: qa-input ################################################### data("qa_080828_081110", package="ShortReadTutorial") qa ################################################### ### chunk number 26: qa-report ################################################### rpt <- report(qa, dest=tempfile()) ################################################### ### chunk number 27: qa-browse eval=FALSE ################################################### ## browseURL(rpt) ################################################### ### chunk number 28: qa-readCount ################################################### qa[["readCounts"]] ################################################### ### chunk number 29: qa-template-source ################################################### ShortRead:::.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, as.prob=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(extdataDir, ".*_export.txt", type="SolexaExport") ## mpi.close.Rslaves() ################################################### ### chunk number 37: workflow-qa eval=FALSE ################################################### ## rpt <- report(qa(extdataDir, "_export.txt"), ## 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(extdataDir, pattern, filter=filt) ################################################### ### chunk number 39: filter-MAQ eval=FALSE ################################################### ## maqDir <- "path/to/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())