## ----style, eval=TRUE, echo=FALSE, results='asis'------------------------ options(max.print=1000) stopifnot(BiocInstaller::biocVersion() == "2.13") BiocStyle::latex() library(knitr) opts_chunk$set(cache=TRUE, tidy=FALSE) 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.5, 0.2, 1.2, 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) } ## ----echo=FALSE---------------------------------------------------------- suppressPackageStartupMessages(library(GenomicRanges)) gr1 <- GRanges(Rle(c("chr1", "chr19", "chrX", "chrM", "chr19"), c(2, 1, 2, 1, 1)), IRanges(c(10001, 5508, 25644, 685, 142501, 33077, 8001), width=c(400, 285, 1620, 333, 1, 2000, 421)), strand=c("-", "*", "+", "+", "-", "*", "+")) ## ------------------------------------------------------------------------ gr1 ## ----echo=FALSE---------------------------------------------------------- suppressPackageStartupMessages(library(BSgenome.Scerevisiae.UCSC.sacCer2)) tiles <- tileGenome(seqinfo(Scerevisiae), tilewidth=50000, cut.last.tile.in.chrom=TRUE) ## ------------------------------------------------------------------------ tiles ## ----echo=FALSE---------------------------------------------------------- gr2 <- gr1 mcols(gr2) <- DataFrame(desc=c("exon", "repeat region", "gene", "exon", "exon", "snp", "repeat region"), GC_content=c(0.55, 0.495, 0.562, 0.53, 0.552, NA, 0.51)) ## ------------------------------------------------------------------------ gr2 ## ----echo=FALSE---------------------------------------------------------- suppressPackageStartupMessages(library(TxDb.Hsapiens.UCSC.hg19.knownGene)) hg19_cds <- cds(TxDb.Hsapiens.UCSC.hg19.knownGene, columns="gene_id") ## ------------------------------------------------------------------------ hg19_cds ## ----echo=FALSE---------------------------------------------------------- while (search()[2L] != "package:parallel") detach(search()[2L], character.only=TRUE) detach("package:parallel") ## ----results='asis'------------------------------------------------------ library(GenomicRanges) ## ----eval=FALSE---------------------------------------------------------- ## class?GRanges ## ----results='hide'------------------------------------------------------ example(GRanges) ## ------------------------------------------------------------------------ gr ## ------------------------------------------------------------------------ seqnames(gr) ## ------------------------------------------------------------------------ class(seqnames(gr)) ## ------------------------------------------------------------------------ width(gr) ## ------------------------------------------------------------------------ all(width(gr) == end(gr) - start(gr) + 1) ## ------------------------------------------------------------------------ names(gr) ## ------------------------------------------------------------------------ names(gr) <- toupper(names(gr)) ## ------------------------------------------------------------------------ mcols(gr) ## ------------------------------------------------------------------------ mcols(gr)$GC ## ------------------------------------------------------------------------ gr[mcols(gr)$GC >= 0.3 & mcols(gr)$GC <= 0.6] ## ------------------------------------------------------------------------ subset(gr, GC >= 0.3 & GC <= 0.6) ## ------------------------------------------------------------------------ mcols(gr)$id <- sprintf("ID%03d", seq_along(gr)) gr ## ----ranges-gr0-plot,results='hide',echo=FALSE,fig.keep='none'----------- gr0 <- IRanges(start=c(7, 9, 12, 14, 22:24), end=c(15, 11, 12, 18, 26, 27, 28)) png("ranges-gr0-plot.png", width=800, height=170) plotRanges(gr0, xlim=c(5, 35), main="gr0", col="blue") dev.off() ## ----ranges-shift-gr0-plot,results='hide',echo=FALSE,fig.keep='none'----- png("ranges-shift-gr0-plot.png", width=800, height=170) plotRanges(shift(gr0, 5), xlim=c(5, 35), main="shift(gr0, 5)", col="blue") dev.off() ## ----ranges-disjoin-gr0-plot,results='hide',echo=FALSE,fig.keep='none'---- png("ranges-disjoin-gr0-plot.png", width=800, height=170) plotRanges(disjoin(gr0), xlim=c(5, 35), main="disjoin(gr0)", col="blue") dev.off() ## ----ranges-reduce-gr0-plot,results='hide',echo=FALSE,fig.keep='none'---- png("ranges-reduce-gr0-plot.png", width=800, height=170) plotRanges(reduce(gr0), xlim=c(5, 35), main="reduce(gr0)", col="blue") dev.off() ## ------------------------------------------------------------------------ gr2 <- shift(gr, 500) ## ------------------------------------------------------------------------ flank(gr2, width=100) ## ------------------------------------------------------------------------ shift ## ------------------------------------------------------------------------ showMethods(shift) ## ------------------------------------------------------------------------ selectMethod(shift, "GRanges") ## ------------------------------------------------------------------------ is(gr, "GenomicRanges") class(gr) ## ----eval=FALSE---------------------------------------------------------- ## ?`shift,GenomicRanges-method` ## ------------------------------------------------------------------------ library(TxDb.Hsapiens.UCSC.hg19.knownGene) TxDb.Hsapiens.UCSC.hg19.knownGene ## ------------------------------------------------------------------------ ex_by_gn <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, by="gene") class(ex_by_gn) ## ------------------------------------------------------------------------ length(ex_by_gn) ## ------------------------------------------------------------------------ "1008" %in% names(ex_by_gn) ## ------------------------------------------------------------------------ ex_by_gn[["1008"]] ## ------------------------------------------------------------------------ seqnames(ex_by_gn[["1008"]]) strand(ex_by_gn[["1008"]]) ## ------------------------------------------------------------------------ nex_per_gn <- elementLengths(ex_by_gn) ## ------------------------------------------------------------------------ max(nex_per_gn) ## ------------------------------------------------------------------------ which(nex_per_gn == max(nex_per_gn)) ## ------------------------------------------------------------------------ which.max(nex_per_gn) ## ------------------------------------------------------------------------ library(parathyroidSE) bamdir <- system.file("extdata", package="parathyroidSE") bampaths <- list.files(bamdir, pattern="bam$", full.names=TRUE) bampaths ## ------------------------------------------------------------------------ library(Rsamtools) reads <- readGAlignmentsFromBam(bampaths[1]) reads ## ------------------------------------------------------------------------ library(parathyroidSE) bamdir <- system.file("extdata", package="parathyroidSE") bampaths <- list.files(bamdir, pattern="bam$", full.names=TRUE) reads <- readGAlignmentsFromBam(bampaths[3]) class(reads) reads ## ------------------------------------------------------------------------ length(grep("N", cigar(reads))) ## ------------------------------------------------------------------------ max(ngap(reads)) table(ngap(reads)) ## ------------------------------------------------------------------------ read_ranges <- as(reads, "GRangesList") ## ------------------------------------------------------------------------ table(elementLengths(read_ranges)) ## ------------------------------------------------------------------------ as(reads, "GRanges") ## ------------------------------------------------------------------------ cvg <- coverage(read_ranges) class(cvg) cvg ## ------------------------------------------------------------------------ cvg[["12"]] max(cvg[["12"]]) ## ------------------------------------------------------------------------ max(cvg)