## ----style, echo = FALSE, results = 'asis'-------------------------------------------------------- options(width=100) knitr::opts_chunk$set( eval=as.logical(Sys.getenv("KNITR_EVAL", "TRUE")), cache=as.logical(Sys.getenv("KNITR_CACHE", "TRUE"))) ## ----setup, echo=FALSE---------------------------------------------------------------------------- suppressPackageStartupMessages({ library(Biostrings) library(GenomicRanges) library(SummarizedExperiment) library(airway) library(rtracklayer) library(ShortRead) library(GenomicAlignments) library(RNAseqData.HNRNPC.bam.chr14) library(VariantAnnotation) }) ## ----install, eval=FALSE-------------------------------------------------------------------------- # source("https://bioconductor.org/biocLite.R") # biocLite(c("DESeq2", "org.Hs.eg.db")) ## ----require-------------------------------------------------------------------------------------- library(GenomicRanges) ## ----help-bioc, eval=FALSE------------------------------------------------------------------------ # help(package="GenomicRanges") # vignette(package="GenomicRanges") # vignette(package="GenomicRanges", "GenomicRangesHOWTOs") # ?GRanges ## ----five-lines----------------------------------------------------------------------------------- x <- rnorm(1000) y <- x + rnorm(1000) df <- data.frame(X=x, Y=y) plot(Y ~ X, df) fit <- lm(Y ~ X, df) anova(fit) abline(fit) ## ----help-r, eval=FALSE--------------------------------------------------------------------------- # class(fit) # methods(class=class(fit)) # methods(plot) # ?"plot" # ?"plot.formula" ## ----classes-and-methods-------------------------------------------------------------------------- library(Biostrings) dna <- DNAStringSet(c("AACAT", "GGCGCCT")) reverseComplement(dna) ## ----classes-and-methods-discovery, eval=FALSE---------------------------------------------------- # class(dna) # ?"DNAStringSet-class" # ?"reverseComplement,DNAStringSet-method" ## ------------------------------------------------------------------------------------------------- library(Biostrings) data(phiX174Phage) phiX174Phage letterFrequency(phiX174Phage, c("A", "C", "G", "T")) letterFrequency(phiX174Phage, "GC", as.prob=TRUE) ## ----ranges, message=FALSE------------------------------------------------------------------------ library(GenomicRanges) gr <- GRanges("A", IRanges(c(10, 20, 22), width=5), "+") shift(gr, 1) # intra-range range(gr) # inter-range reduce(gr) # inter-range snps <- GRanges("A", IRanges(c(11, 17, 24), width=1)) findOverlaps(snps, gr) # between-range setdiff(range(gr), gr) # 'introns' ## ------------------------------------------------------------------------------------------------- library(airway) data(airway) airway ## ------------------------------------------------------------------------------------------------- x <- assay(airway) class(x) dim(x) head(x) colData(airway) rowRanges(airway) ## ------------------------------------------------------------------------------------------------- cidx <- colData(airway)$dex %in% "trt" airway[, cidx] ## shortcut airway[, airway$dex %in% "trt"] ## ------------------------------------------------------------------------------------------------- chr14 <- as(seqinfo(rowRanges(airway)), "GRanges")["14"] ridx <- rowRanges(airway) %over% chr14 airway[ridx,] ## shortcut chr14 <- as(seqinfo(airway), "GRanges")["14"] airway[airway %over% chr14,] ## ------------------------------------------------------------------------------------------------- library(RNAseqData.HNRNPC.bam.chr14) fls = RNAseqData.HNRNPC.bam.chr14_BAMFILES basename(fls) ## ------------------------------------------------------------------------------------------------- library(GenomicAlignments) bfls = BamFileList(fls) bfl = bfls[[1]] ## ------------------------------------------------------------------------------------------------- ga = readGAlignments(bfl) ga table(strand(ga)) ## ------------------------------------------------------------------------------------------------- tail(sort(table(cigar(ga)))) ga[cigar(ga) != "72M"] ## ------------------------------------------------------------------------------------------------- summarizeJunctions(ga) junctions <- summarizeJunctions(ga, with.revmap=TRUE) ga[ junctions$revmap[[1]] ] ## ------------------------------------------------------------------------------------------------- coverage(bfl)$chr14