### R code from vignette source 'ExploringSequences-lab.Rnw' ################################################### ### code chunk number 1: install (eval = FALSE) ################################################### ## pkg <- "http://10.0.1.2/ExploringSequences_0.1.0.tar.gz" ## install.packages(pkg, repos=NULL, type="source") ################################################### ### code chunk number 2: setup ################################################### library(ExploringSequences) ################################################### ### code chunk number 3: browseVignettese (eval = FALSE) ################################################### ## browseVignettes(package="ExploringSequences") ################################################### ### code chunk number 4: setup-filepath ################################################### fastqFile <- system.file("data", "fq.Rda", package="ExploringSequences") bamFiles <- system.file("extdata", "bam", package="ExploringSequences") barFile <- system.file("data", "bar.Rda", package="ExploringSequences") ################################################### ### code chunk number 5: sanity-check ################################################### stopifnot(file.exists(fastqFile)) stopifnot(file.exists(bamFiles)) stopifnot(file.exists(barFile)) ################################################### ### code chunk number 6: fastq-input ################################################### load(fastqFile) fq ################################################### ### code chunk number 7: fastq-input ################################################### class(fq) # class ShortReadQ sread(fq) # short reads quality(fq) # quality scores ################################################### ### code chunk number 8: alphabet-mono ################################################### mono <- alphabetFrequency(sread(fq), baseOnly=TRUE, collapse=TRUE) mono ################################################### ### code chunk number 9: alphabet-dinuc-fun ################################################### di <- colSums(dinucleotideFrequency(sread(fq))) ################################################### ### code chunk number 10: alphabet-gc-fun ################################################### gcContent <- function(x) { abc <- alphabetFrequency(sread(x), baseOnly=TRUE) gcPerRead <- rowSums(abc[,c("G", "C")]) wd <- unique(width(sread(x))) tabulate(1 + gcPerRead, 1 + wd) } ################################################### ### code chunk number 11: alphabet-gc ################################################### gc <- gcContent(fq) gc ################################################### ### code chunk number 12: alphabet-dinuc-longdf ################################################### bins <- seq_along(gc) - 1 df <- data.frame(Count=unlist(gc, use.names=FALSE), GC=bins / bins[length(bins)]) ################################################### ### code chunk number 13: alphabet-dinuc-figure ################################################### xyplot(Count~GC, df, type=c("g", "b"), pch=20, file="gc") ################################################### ### code chunk number 14: dinuc-exp ################################################### ## outer product of mononucleotides frequencies f <- mono[-5] / sum(mono[-5]) exp0 <- as.vector(outer(f, f)) ## expected values, as counts n <- sum(di) exp <- exp0 * n mode(exp) <- "integer" head(exp, 3) ## Chi-squared test chisq.test(cbind(di, exp)) ################################################### ### code chunk number 15: dinuc-display ################################################### ## display df <- data.frame(Observed=as.vector(di), Expected=as.vector(exp), Label=names(di), stringsAsFactors=FALSE) xyplot(Observed - Expected ~ Expected, df, aspect="iso", panel=function(...) { panel.abline(h=0, col="gray") panel.text(..., labels=df$Label) }, file="dinuc") ################################################### ### code chunk number 16: abc-read ################################################### nuc <- c("A", "C", "G", "T", "N") abc <- alphabetByCycle(sread(fq))[nuc,] abc[,1:5] ################################################### ### code chunk number 17: abc-display ################################################### ## create a 'long' data frame ncycle <- ncol(abc) cycle <- rep(seq_len(ncycle), each=nrow(abc)) df <- data.frame(Count=as.vector(abc), Nucleotide=factor(nuc, levels=nuc), Cycle=cycle) ## plot Count as a function of cycle, with Nucleotide used to group lines ## within a panel. xyplot(Count ~ Cycle, group=Nucleotide, df, type="l", key=simpleKey(lines=TRUE, points=FALSE, columns=5, text=nuc), file="abc-read") ################################################### ### code chunk number 18: abc-quality ################################################### abc <- colMeans(as(quality(fq), "matrix")) ################################################### ### code chunk number 19: abc-quality-display ################################################### df <- data.frame(Mean=abc, Cycle=seq_along(abc)) xyplot(Mean ~ Cycle, df, type=c("g", "b"), pch=20, file="abc-quality") ################################################### ### code chunk number 20: quality-by-read ################################################### qual <- rowMeans(as(quality(fq), "matrix")) df <- data.frame(Quality=qual) histogram(~Quality, df, file="quality-by-read") ################################################### ### code chunk number 21: freqseq ################################################### t <- tables(sread(fq), n=0)[["distribution"]] df <- data.frame(nOccurrences=t$nOccurrences, CummulativeReads=cumsum(t$nOccurrences * t$nReads)) ## display xyplot(CummulativeReads ~ log(nOccurrences), df, type=c("b", "g"), pch=20, file="freqseq") ################################################### ### code chunk number 22: qa-path ################################################### qaFile <- system.file("data", "qa_GSM461176_81.rda", package="ExploringSequences") ################################################### ### code chunk number 23: qa-solution (eval = FALSE) ################################################### ## load(qaFile) ## rpt <- report(qa_GSM461176_81) ## browseURL(rpt) ################################################### ### code chunk number 24: qa-report (eval = FALSE) ################################################### ## rpt <- system.file("GSM461176_81_qa_report", "index.html", ## package="ExploringSequences") ## browseURL(rpt) ################################################### ### code chunk number 25: bam-open ################################################### fls <- list.files(bamFiles, "bam$", full=TRUE) stopifnot(length(fls) == 7) bam <- open(BamFileList(fls)) names(bam) <- sub("_subset.bam", "", basename(fls)) ################################################### ### code chunk number 26: bam-header (eval = FALSE) ################################################### ## seqinfo(bam[[1]]) ## h <- scanBamHeader(bam[[1]])[["text"]] ## noquote(unname(sapply(h[["@PG"]], strwrap))) ################################################### ### code chunk number 27: count-bam ################################################### cnt <- countBam(bam) cnt ################################################### ### code chunk number 28: count-bam-regions ################################################### which <- GRanges(c("chr3L", "chrX", "chr3L"), IRanges(c(1871574, 10675019, 14769596), c(1876336, 10680978, 14779523))) param <- ScanBamParam(which=which) head(cnts <- countBam(bam, param=param)) ## create a label, e.g., chr3L:1871574-1876336 cnts$rgn <- with(cnts, sprintf("%s:%d-%d", space, start, end)) xyplot(file~records, group=rgn, cnts, type="b", pch=20, auto.key=list(lines=TRUE, points=FALSE, columns=2), file="bam-regions") ################################################### ### code chunk number 29: flags ################################################### param1 <- ScanBamParam(which=which[3]) countFlags(bam, param=param1) bam2 <- bam[-c(1, 4, 5)] ################################################### ### code chunk number 30: insert-param ################################################### param2 <- ScanBamParam( what=c("pos", "mpos"), which=which[3], flag=scanBamFlag(isProperPair=TRUE, isFirstMateRead=TRUE)) ################################################### ### code chunk number 31: insert-scanBam ################################################### iwd0 <- lapply(bam2, scanBam, param=param2) str(iwd0) ################################################### ### code chunk number 32: insert-calc ################################################### fun <- function(elt) with(elt[[1]], abs(pos - mpos)) str(iwd <- lapply(iwd0, fun)) ################################################### ### code chunk number 33: insert-figure ################################################### df <- do.call(make.groups, iwd) densityplot(~data, group=which, df, plot.points=FALSE, auto.key=list(lines=TRUE, points=FALSE, columns=2), file="insert") ################################################### ### code chunk number 34: bar-input ################################################### load(barFile) bar summary(width(bar)) ################################################### ### code chunk number 35: bar-narrow ################################################### bar254 <- bar[width(bar) == 254] alphabetByCycle(narrow(sread(bar254), 101, 110))[1:4,] ################################################### ### code chunk number 36: bar-count ################################################### codes0 <- narrow(sread(bar), 1, 8) codes <- as.character(codes0) cnt <- sort(table(codes), decreasing=TRUE)[1:5] cnt aBar <- bar[codes==names(cnt)[[1]]] aBar ################################################### ### code chunk number 37: bar-trim ################################################### noBar <- narrow(aBar, 11, width(aBar)) noBar ################################################### ### code chunk number 38: bar-primer ################################################### pcrPrimer <- "GGACTACCVGGGTATCTAAT" ################################################### ### code chunk number 39: bar-trim ################################################### trimmed <- trimLRPatterns(pcrPrimer, subject=noBar, Lfixed=FALSE) trimmed table(width(noBar) - width(trimmed)) ################################################### ### code chunk number 40: qa-report (eval = FALSE) ################################################### ## fun <- function(fl) { ## id <- sub(".fastq$", "", basename(fl)) ## qa(readFastq(fl), id) ## } ## fls <- list.files(pattern=".fastq$", full=TRUE) ## doit <- ## if (suppressWarnings(require("multicore"))) { ## mclapply ## } else lapply ## qas <- doit(fls, fun) ## qa_GSM461176_81 <- do.call(rbind, qas)