### R code from vignette source 'A3_AlignmentsAndRanges.Rnw' ################################################### ### code chunk number 1: setup ################################################### library(GenomicFeatures) plotRanges <- function(x, xlim = x, main = deparse(substitute(x)), col = "black", sep = 0.5, ...) { height <- 1 if (is(xlim, "Ranges")) xlim <- c(min(start(xlim)), max(end(xlim))) bins <- disjointBins(IRanges(start(x), end(x) + 1)) plot.new() par(mai=c(0.6, 0.2, 0.6, 0.2)) plot.window(xlim, c(0, max(bins)*(height + sep))) ybottom <- bins * (sep + height) - height rect(start(x)-0.5, ybottom, end(x)+0.5, ybottom + height, col = col, ...) title(main, cex.main=2.8, font.main=1) axis(1) } ################################################### ### code chunk number 2: GRanges-genes ################################################### genes <- GRanges(seqnames=c("3R", "X"), ranges=IRanges( start=c(19967117, 18962306), end=c(19973212, 18962925)), strand=c("+", "-"), seqlengths=c(`3R`=27905053L, `X`=22422827L)) ################################################### ### code chunk number 3: GRanges-display ################################################### genes ################################################### ### code chunk number 4: GRanges-help (eval = FALSE) ################################################### ## ?GRanges ################################################### ### code chunk number 5: GRanges-vignettes (eval = FALSE) ################################################### ## vignette(package="GenomicRanges") ################################################### ### code chunk number 6: GRanges-basics ################################################### genes[2] strand(genes) width(genes) length(genes) names(genes) <- c("FBgn0039155", "FBgn0085359") genes # now with names ################################################### ### code chunk number 7: ranges-ir ################################################### ir <- IRanges(start=c(7, 9, 12, 14, 22:24), end=c(15, 11, 12, 18, 26, 27, 28)) ################################################### ### code chunk number 8: ranges-ir-plot ################################################### png("ranges-ir-plot.png", width=800, height=160) plotRanges(ir, xlim=c(5, 35), main="Original") dev.off() png("ranges-shift-ir-plot.png", width=800, height=160) plotRanges(shift(ir, 5), xlim=c(5, 35), main="Shift") dev.off() png("ranges-reduce-ir-plot.png", width=800, height=160) plotRanges(reduce(ir), xlim=c(5, 35), main="Reduce") dev.off() ################################################### ### code chunk number 9: ranges-shift-ir ################################################### shift(ir, 5) ################################################### ### code chunk number 10: ranges-shift-ir-plot ################################################### png("ranges-shift-ir-plot.png", width=800, height=160) plotRanges(shift(ir, 5), xlim=c(5, 35), main="Shift") dev.off() ################################################### ### code chunk number 11: ranges-reduce-ir ################################################### reduce(ir) ################################################### ### code chunk number 12: ranges-reduce-ir-plot ################################################### png("ranges-reduce-ir-plot.png", width=800, height=160) plotRanges(reduce(ir), xlim=c(5, 35), main="Reduce") dev.off() ################################################### ### code chunk number 13: coverage ################################################### coverage(ir) ################################################### ### code chunk number 14: GRanges-mcols ################################################### mcols(genes) <- DataFrame(EntrezId=c("42865", "2768869"), Symbol=c("kal-1", "CG34330")) ################################################### ### code chunk number 15: GRanges-metadata ################################################### metadata(genes) <- list(CreatedBy="A. User", Date=date()) ################################################### ### code chunk number 16: GRangesList-eg-setup ################################################### library(TxDb.Dmelanogaster.UCSC.dm3.ensGene) txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene # alias fbgn <- exonsBy(txdb, "gene")["FBgn0039155"] seqlevels(fbgn) <- "chr3R" ################################################### ### code chunk number 17: GRangesList-eg ################################################### fbgn ################################################### ### code chunk number 18: txdb ################################################### library(TxDb.Dmelanogaster.UCSC.dm3.ensGene) txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene # alias ex0 <- exonsBy(txdb, "gene") head(table(elementLengths(ex0))) ids <- c("FBgn0002183", "FBgn0003360", "FBgn0025111", "FBgn0036449") ex <- reduce(ex0[ids]) ################################################### ### code chunk number 19: SAM ################################################### fl <- system.file("extdata", "ex1.sam", package="Rsamtools") strsplit(readLines(fl, 1), "\t")[[1]] ################################################### ### code chunk number 20: readGappedAlignments ################################################### alnFile <- system.file("extdata", "ex1.bam", package="Rsamtools") aln <- readGappedAlignments(alnFile) head(aln, 3) ################################################### ### code chunk number 21: GappedAlignments-accessors ################################################### table(strand(aln)) table(width(aln)) head(sort(table(cigar(aln)), decreasing=TRUE)) ################################################### ### code chunk number 22: bam-ex-fls ################################################### bigdata <- "~/BigData" fls <- dir(file.path(bigdata, "bam"), ".bam$", full=TRUE) #$ names(fls) <- sub("_subset.bam", "", basename(fls)) ################################################### ### code chunk number 23: bam-ex-input ################################################### ## input aln <- readGappedAlignments(fls[1]) xtabs( ~ seqnames + strand, as.data.frame(aln)) ################################################### ### code chunk number 24: bam-ex-roi ################################################### data(ex) # from a later exercise ################################################### ### code chunk number 25: bam-ex-strand ################################################### strand(aln) <- "*" # protocol not strand-aware ################################################### ### code chunk number 26: bam-ex-hits ################################################### hits <- countOverlaps(aln, ex) table(hits) ################################################### ### code chunk number 27: bam-ex-cnt ################################################### cnt <- countOverlaps(ex, aln[hits==1]) ################################################### ### code chunk number 28: bam-count-fun ################################################### counter <- function(filePath, range) { aln <- readGappedAlignments(filePath) strand(aln) <- "*" hits <- countOverlaps(aln, range) countOverlaps(range, aln[hits==1]) } ################################################### ### code chunk number 29: bam-count-all ################################################### counts <- sapply(fls, counter, ex) ################################################### ### code chunk number 30: bam-count-mclapply (eval = FALSE) ################################################### ## if (require(parallel)) ## simplify2array(mclapply(fls, counter, ex))