### R code from vignette source 'GenomicFiles.Rnw' ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex() ################################################### ### code chunk number 2: biocLite (eval = FALSE) ################################################### ## source("http://bioconductor.org/biocLite.R") ## biocLite("GenomicFiles") ################################################### ### code chunk number 3: quick_start-load ################################################### library(GenomicFiles) ################################################### ### code chunk number 4: quick_start-ranges ################################################### gr <- GRanges("chr14", IRanges(c(19411500 + (1:5)*20), width=10)) ################################################### ### code chunk number 5: class-bam-data ################################################### library(RNAseqData.HNRNPC.bam.chr14) fls <- RNAseqData.HNRNPC.bam.chr14_BAMFILES ################################################### ### code chunk number 6: quick_start-MAP ################################################### MAP <- function(range, file, ...) { requireNamespace("Rsamtools") Rsamtools::pileup(file, scanBamParam=Rsamtools::ScanBamParam(which=range)) } ################################################### ### code chunk number 7: quick_start-reduceByFile ################################################### se <- reduceByFile(gr, fls, MAP, summarize=TRUE) se ################################################### ### code chunk number 8: quick_start-assays ################################################### dim(assays(se)$data) ## ranges x files ################################################### ### code chunk number 9: quick_start-MAP-REDUCE-reduceByRange ################################################### REDUCE <- function(mapped, ...) { cmb = do.call(rbind, mapped) xtabs(count ~ pos + nucleotide, cmb) } lst <- reduceByRange(gr, fls, MAP, REDUCE, iterate=FALSE) ################################################### ### code chunk number 10: quick_start-result ################################################### head(lst[[1]], 3) ################################################### ### code chunk number 11: overview-GenomicFiles ################################################### GenomicFiles(gr, fls) ################################################### ### code chunk number 12: pileups-ranges ################################################### gr <- GRanges("chr14", IRanges(c(19411677, 19659063, 105421963, 105613740), width=20)) ################################################### ### code chunk number 13: pileups-MAP ################################################### MAP <- function(range, file, ...) { requireNamespace("deepSNV") ct = deepSNV::bam2R(file, GenomeInfoDb::seqlevels(range), GenomicRanges::start(range), GenomicRanges::end(range), q=0) ct[, c("A", "T", "C", "G", "a", "t", "c", "g")] } ################################################### ### code chunk number 14: pileups-REDUCE ################################################### REDUCE <- function(mapped, ...) Reduce("+", mapped) ################################################### ### code chunk number 15: pileups-reduceByRange ################################################### pile2 <- reduceByRange(gr, fls, MAP, REDUCE) length(pile2) elementNROWS(pile2) ################################################### ### code chunk number 16: pileups-res ################################################### head(pile2[[1]]) ################################################### ### code chunk number 17: ttest-ranges ################################################### roi <- GRanges("chr14", IRanges(c(19411677, 19659063, 105421963, 105613740), width=20)) ################################################### ### code chunk number 18: ttest-group ################################################### grp <- factor(rep(c("A","B"), each=length(fls)/2)) ################################################### ### code chunk number 19: ttest-MAP ################################################### MAP <- function(range, file, ...) { requireNamespace("GenomicAlignments") param <- Rsamtools::ScanBamParam(which=range) as.numeric(unlist( GenomicAlignments::coverage(file, param=param)[range], use.names=FALSE)) } ################################################### ### code chunk number 20: ttest-REDUCE ################################################### REDUCE <- function(mapped, ..., grp) { mat = simplify2array(mapped) idx = which(rowSums(mat) != 0) df = genefilter::rowttests(mat[idx,], grp) cbind(offset = idx - 1, df) } ################################################### ### code chunk number 21: ttest-results (eval = FALSE) ################################################### ## ttest <- reduceByRange(roi, fls, MAP, REDUCE, iterate=FALSE, grp=grp) ################################################### ### code chunk number 22: junctions-ranges ################################################### gr <- GRanges("chr14", IRanges(c(19100000, 106000000), width=1e7)) ################################################### ### code chunk number 23: junctions-MAP ################################################### MAP <- function(range, file, ...) { requireNamespace("GenomicAlignments") ## for readGAlignments() ## ScanBamParam() param = Rsamtools::ScanBamParam(which=range) gal = GenomicAlignments::readGAlignments(file, param=param) table(GenomicAlignments::njunc(gal)) } ################################################### ### code chunk number 24: junctions-GenomicFiles ################################################### gf <- GenomicFiles(gr, fls) gf ################################################### ### code chunk number 25: junctions-counts1 ################################################### counts1 <- reduceByFile(gf[,1:3], MAP=MAP) length(counts1) ## 3 files elementNROWS(counts1) ## 2 ranges ################################################### ### code chunk number 26: junctions-counts1-show ################################################### counts1[[1]] ################################################### ### code chunk number 27: junctions-REDUCE ################################################### REDUCE <- function(mapped, ...) sum(sapply(mapped, "[", "1")) reduceByFile(gr, fls, MAP, REDUCE) ################################################### ### code chunk number 28: junctions-counts2 ################################################### counts2 <- reduceFiles(gf[,1:3], MAP=MAP) ################################################### ### code chunk number 29: junctions-counts2-show ################################################### ## reduceFiles returns counts for all ranges. counts2[[1]] ## reduceByFile returns counts for each range separately. counts1[[1]] ################################################### ### code chunk number 30: coverage1-tiles ################################################### chr14_seqlen <- seqlengths(seqinfo(BamFileList(fls))["chr14"]) tiles <- tileGenome(chr14_seqlen, ntile=5) ################################################### ### code chunk number 31: coverage1-tiles-show ################################################### tiles ################################################### ### code chunk number 32: coverage1-MAP ################################################### MAP = function(range, file, ...) { requireNamespace("GenomicAlignments") ## for ScanBamParam() and coverage() param = Rsamtools::ScanBamParam(which=range) rle = GenomicAlignments::coverage(file, param=param)[range] c(width = GenomicRanges::width(range), sum = sum(S4Vectors::runLength(rle) * S4Vectors::runValue(rle))) } ################################################### ### code chunk number 33: coverage1-REDUCE ################################################### REDUCE = function(mapped, ...) { Reduce(function(i, j) Map("+", i, j), mapped) } ################################################### ### code chunk number 34: coverage1-results (eval = FALSE) ################################################### ## cvg1 <- reduceByFile(tiles, fls, MAP, REDUCE, iterate=TRUE) ################################################### ### code chunk number 35: coverage2-MAP ################################################### MAP = function(range, file, ...) { requireNamespace("GenomicAlignments") ## for ScanBamParam() and coverage() GenomicAlignments::coverage( file, param=Rsamtools::ScanBamParam(which=range))[range] } ################################################### ### code chunk number 36: coverage2-REDUCE ################################################### REDUCE = function(mapped, ...) { sapply(mapped, Reduce, f = "+") } ################################################### ### code chunk number 37: coverage2-results ################################################### cvg2 <- reduceFiles(unlist(tiles), fls, MAP, REDUCE) ################################################### ### code chunk number 38: coverage2-show ################################################### cvg2[1] ################################################### ### code chunk number 39: coverage3-MAP ################################################### MAP = function(range, file, ...) { requireNamespace("BiocParallel") ## for bplapply() nranges = 2 idx = split(seq_along(range), ceiling(seq_along(range)/nranges)) BiocParallel::bplapply(idx, function(i, range, file) { requireNamespace("GenomicAlignments") ## ScanBamParam(), coverage() chunk = range[i] param = Rsamtools::ScanBamParam(which=chunk) cvg = GenomicAlignments::coverage(file, param=param)[chunk] Reduce("+", cvg) ## collapse coverage within chunks }, range, file) } ################################################### ### code chunk number 40: coverage3-REDUCE ################################################### REDUCE = function(mapped, ...) { sapply(mapped, Reduce, f = "+") } ################################################### ### code chunk number 41: coverage3-results (eval = FALSE) ################################################### ## cvg3 <- reduceFiles(unlist(tiles), fls, MAP, REDUCE) ################################################### ### code chunk number 42: reduceByYield-YIELD ################################################### library(GenomicAlignments) bf <- BamFile(fls[1], yieldSize=100000) YIELD <- function(x, ...) readGAlignments(x) ################################################### ### code chunk number 43: reduceByYield-MAP-REDUCE ################################################### gr <- unlist(tiles, use.names=FALSE) MAP <- function(value, gr, ...) { requireNamespace("GenomicRanges") ## for countOverlaps() GenomicRanges::countOverlaps(gr, value) } REDUCE <- `+` ################################################### ### code chunk number 44: reduceByYield-DONE ################################################### DONE <- function(value) length(value) == 0L ################################################### ### code chunk number 45: reduceByYield-bplapply ################################################### FUN <- function(file, gr, YIELD, MAP, REDUCE, tiles, ...) { requireNamespace("GenomicAlignments") ## for BamFile, readGAlignments() requireNamespace("GenomicFiles") ## for reduceByYield() gr <- unlist(tiles, use.names=FALSE) bf <- Rsamtools::BamFile(file, yieldSize=100000) YIELD <- function(x, ...) GenomicAlignments::readGAlignments(x) MAP <- function(value, gr, ...) { requireNamespace("GenomicRanges") ## for countOverlaps() GenomicRanges::countOverlaps(gr, value) } REDUCE <- `+` GenomicFiles::reduceByYield(bf, YIELD, MAP, REDUCE, gr=gr, parallel=FALSE) } ################################################### ### code chunk number 46: sessionInfo ################################################### toLatex(sessionInfo())