### R code from vignette source 'Ranges-slides.Rnw' ################################################### ### code chunk number 1: setup ################################################### options(width=84) 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) } ################################################### ### code chunk number 2: IRanges_example1 ################################################### library(IRanges) ir1 <- IRanges(start=c(12, -9, NA, 12), end=c(NA, 0, 15, NA), width=c(4, NA, 4, 3)) ir1 # "show" method not yet consistent with the other "show" methods (TODO) start(ir1) end(ir1) width(ir1) successiveIRanges(c(10, 5, 38), from=101) ################################################### ### code chunk number 3: GRanges_accessors3 ################################################### names(ir1) <- LETTERS[1:4] names(ir1) mcols(ir1) <- DataFrame(score=11:14, GC=seq(1, 0, length=4)) mcols(ir1) ir1 ################################################### ### code chunk number 4: IRanges_example2 ################################################### ir1[-2] ir2 <- c(ir1, IRanges(-10, 0)) ir2 ################################################### ### code chunk number 5: IRanges_example3 ################################################### duplicated(ir2) unique(ir2) ################################################### ### code chunk number 6: IRanges_sort ################################################### order(ir2) sort(ir2) ################################################### ### code chunk number 7: IRanges_example12 ################################################### ok <- c(FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE) ir4 <- as(ok, "IRanges") # from logical vector to IRanges ir4 ################################################### ### code chunk number 8: IRanges_example12 ################################################### as.data.frame(ir4) ################################################### ### code chunk number 9: ranges-ir0-plot ################################################### ir0 <- IRanges(start=c(7, 9, 12, 14, 22:24), end=c(15, 11, 12, 18, 26, 27, 28)) png("ranges-ir0-plot.png", width=800, height=170) plotRanges(ir0, xlim=c(5, 35), main="ir0", col="blue") dev.off() ################################################### ### code chunk number 10: ranges-shift-ir0-plot ################################################### png("ranges-shift-ir0-plot.png", width=800, height=170) plotRanges(shift(ir0, 5), xlim=c(5, 35), main="shift(ir0, 5)", col="blue") dev.off() ################################################### ### code chunk number 11: ranges-disjoin-ir0-plot ################################################### png("ranges-disjoin-ir0-plot.png", width=800, height=170) plotRanges(disjoin(ir0), xlim=c(5, 35), main="disjoin(ir0)", col="blue") dev.off() ################################################### ### code chunk number 12: ranges-reduce-ir0-plot ################################################### png("ranges-reduce-ir0-plot.png", width=800, height=170) plotRanges(reduce(ir0), xlim=c(5, 35), main="reduce(ir0)", col="blue") dev.off() ################################################### ### code chunk number 13: IRanges_example4 ################################################### ir1 ################################################### ### code chunk number 14: IRanges_example5 ################################################### shift(ir1, -start(ir1)) flank(ir1, 10, start=FALSE) ################################################### ### code chunk number 15: IRanges_example6 ################################################### ir1 ################################################### ### code chunk number 16: IRanges_example7 ################################################### range(ir1) reduce(ir1) ################################################### ### code chunk number 17: IRanges_example8 ################################################### union(ir1, IRanges(-2, 6)) intersect(ir1, IRanges(-2, 13)) setdiff(ir1, IRanges(-2, 13)) ################################################### ### code chunk number 18: IRanges_example9 ################################################### ir3 <- IRanges(5:1, width=12) ir3 ################################################### ### code chunk number 19: IRanges_example10 ################################################### ir2 ################################################### ### code chunk number 20: IRanges_example11 ################################################### pintersect(ir3, ir2, resolve.empty="max.start") ################################################### ### code chunk number 21: GRanges_constructor ################################################### library(GenomicRanges) gr1 <- GRanges(seqnames=Rle(c("ch1", "chMT"), lengths=c(2, 4)), ranges=IRanges(start=16:21, end=20), strand=rep(c("+", "-", "*"), 2)) gr1 ################################################### ### code chunk number 22: GRanges_accessors1 ################################################### length(gr1) seqnames(gr1) ranges(gr1) ################################################### ### code chunk number 23: GRanges_accessors2 ################################################### start(gr1) end(gr1) width(gr1) strand(gr1) strand(gr1) <- c("-", "-", "+") strand(gr1) ################################################### ### code chunk number 24: GRanges_accessors3 ################################################### names(gr1) <- LETTERS[1:6] names(gr1) mcols(gr1) <- DataFrame(score=11:16, GC=seq(1, 0, length=6)) mcols(gr1) gr1 ################################################### ### code chunk number 25: GRanges_accessors4 ################################################### seqinfo(gr1) seqlevels(gr1) seqlengths(gr1) seqlengths(gr1) <- c(50000, 800) seqlengths(gr1) ################################################### ### code chunk number 26: GRanges_Vector_ops1 ################################################### gr1[c("F", "A")] gr1[strand(gr1) == "+"] ################################################### ### code chunk number 27: GRanges_Vector_ops2 ################################################### gr1 <- gr1[-5] gr1 ################################################### ### code chunk number 28: GRanges_Vector_ops3 ################################################### gr2 <- GRanges(seqnames="ch2", ranges=IRanges(start=c(2:1,2), width=6), score=15:13, GC=seq(0, 0.4, length=3)) gr12 <- c(gr1, gr2) gr12 ################################################### ### code chunk number 29: GRanges_Vector_ops4 ################################################### gr12[length(gr12)] == gr12 duplicated(gr12) unique(gr12) ################################################### ### code chunk number 30: GRanges_Vector_ops5 ################################################### sort(gr12) ################################################### ### code chunk number 31: GRanges_Ranges_ops1 ################################################### gr2 shift(gr2, 50) narrow(gr2, start=2, end=-2) ################################################### ### code chunk number 32: GRanges_Ranges_ops2 ################################################### gr1 resize(gr1, 12) ################################################### ### code chunk number 33: GRanges_Ranges_ops3 ################################################### gr1 flank(gr1, 3) ################################################### ### code chunk number 34: GRanges_Ranges_ops4 ################################################### gr3 <- shift(gr1, c(35000, rep(0, 3), 100)) width(gr3)[c(3,5)] <- 117 gr3 range(gr3) ################################################### ### code chunk number 35: GRanges_Ranges_ops5 ################################################### gr3 disjoin(gr3) ################################################### ### code chunk number 36: GRanges_Ranges_ops6 ################################################### gr3 reduce(gr3) ################################################### ### code chunk number 37: GRanges_Ranges_ops7 ################################################### gr3 gaps(gr3) ################################################### ### code chunk number 38: GRanges_split ################################################### split(gr3, seqnames(gr3)) ################################################### ### code chunk number 39: GRangesList_constructor ################################################### grl <- GRangesList(gr3, gr2) grl ################################################### ### code chunk number 40: GRangesList_accessors1 ################################################### length(grl) ################################################### ### code chunk number 41: GRangesList_accessors2 ################################################### seqnames(grl) ################################################### ### code chunk number 42: GRangesList_accessors3 ################################################### strand(grl) ################################################### ### code chunk number 43: GRangesList_accessors4 ################################################### ranges(grl) ################################################### ### code chunk number 44: GRangesList_accessors5 ################################################### start(grl) end(grl) width(grl) ################################################### ### code chunk number 45: GRangesList_accessors6 ################################################### names(grl) <- c("TX1", "TX2") grl ################################################### ### code chunk number 46: GRangesList_accessors7 ################################################### mcols(grl)$geneid <- c("GENE1", "GENE2") mcols(grl) grl ################################################### ### code chunk number 47: GRangesList_accessors8 ################################################### seqinfo(grl) ################################################### ### code chunk number 48: GRangesList_Vector_ops1 ################################################### grl[c("TX2", "TX1")] ################################################### ### code chunk number 49: GRangesList_Vector_ops2 ################################################### c(grl, GRangesList(gr3)) ################################################### ### code chunk number 50: GRangesList_List_ops1 ################################################### grl[[2]] elementLengths(grl) unlisted <- unlist(grl, use.names=FALSE) # same as c(grl[[1]], grl[[2]]) unlisted ################################################### ### code chunk number 51: GRangesList_List_ops2 ################################################### grl100 <- relist(shift(unlisted, 100), grl) grl100 ################################################### ### code chunk number 52: GRangesList_List_ops3 ################################################### grl100b <- endoapply(grl, shift, 100) grl100b mcols(grl100) mcols(grl100b) ################################################### ### code chunk number 53: GRangesList_Ranges_ops1 ################################################### grl ################################################### ### code chunk number 54: GRangesList_Ranges_ops2 ################################################### shift(grl, 100) ################################################### ### code chunk number 55: GRangesList_Ranges_ops3 ################################################### grl ################################################### ### code chunk number 56: GRangesList_Ranges_ops4 ################################################### flank(grl, 10) ################################################### ### code chunk number 57: GRangesList_Ranges_ops5 ################################################### grl ################################################### ### code chunk number 58: GRangesList_Ranges_ops6 ################################################### range(grl) ################################################### ### code chunk number 59: GRangesList_Ranges_ops7 ################################################### grl ################################################### ### code chunk number 60: GRangesList_Ranges_ops8 ################################################### reduce(grl) ################################################### ### code chunk number 61: GRangesList_Ranges_ops9 ################################################### grl2 <- grl grl2[[1]] <- grl2[[1]][3]; grl2[[2]] <- grl2[[2]][1] grl3 <- unname(grl2) grl3[[1]] <- narrow(unname(grl3[[1]]), start=5, end=-5) ################################################### ### code chunk number 62: GRangesList_Ranges_ops10 ################################################### grl2 grl3 ################################################### ### code chunk number 63: GRangesList_Ranges_ops11 ################################################### psetdiff(grl2, grl3) ################################################### ### code chunk number 64: GAlignments_constructor ################################################### gal0 <- GAlignments(seqnames=Rle(c("ch1", "ch2"), c(3, 1)), pos=1L + 10L*0:3, cigar=c("36M", "20M3D16M", "20M703N16M", "14M2I20M"), strand=strand(c("+", "-", "-", "+"))) gal0 ################################################### ### code chunk number 65: readGAlignments ################################################### library(pasillaBamSubset) U1gal <- readGAlignments(untreated1_chr4()) length(U1gal) head(U1gal) ################################################### ### code chunk number 66: GAlignments_accessors1 ################################################### seqnames(U1gal) table(as.factor(seqnames(U1gal))) strand(U1gal) table(as.factor(strand(U1gal))) head(cigar(U1gal)) head(qwidth(U1gal)) table(qwidth(U1gal)) ################################################### ### code chunk number 67: GAlignments_accessors2 ################################################### head(start(U1gal)) head(end(U1gal)) head(width(U1gal)) head(ngap(U1gal)) table(ngap(U1gal)) ################################################### ### code chunk number 68: GAlignments_accessors3 ################################################### mcols(U1gal) seqinfo(U1gal) ################################################### ### code chunk number 69: readGAlignments2 ################################################### param <- ScanBamParam(what=c("flag", "mapq"), tag=c("NH", "NM")) U1gal <- readGAlignments(untreated1_chr4(), use.names=TRUE, param=param) U1gal[1:5] any(duplicated(names(U1gal))) ################################################### ### code chunk number 70: from_GAlignments_to_GRanges ################################################### as(U1gal, "GRanges") ################################################### ### code chunk number 71: from_GAlignments_to_GRangesList1 ################################################### U1grl <- as(U1gal, "GRangesList") U1grl ################################################### ### code chunk number 72: from_GAlignments_to_GRangesList2 ################################################### all(elementLengths(U1grl) == ngap(U1gal) + 1) ################################################### ### code chunk number 73: readGAlignmentPairs ################################################### library(pasillaBamSubset) U3galp <- readGAlignmentPairs(untreated3_chr4()) length(U3galp) head(U3galp) ################################################### ### code chunk number 74: GAlignmentPairs_accessors1 ################################################### head(first(U3galp)) head(last(U3galp)) ################################################### ### code chunk number 75: GAlignmentPairs_accessors2 ################################################### seqnames(U3galp) strand(U3galp) head(ngap(U3galp)) table(ngap(U3galp)) ################################################### ### code chunk number 76: from_GAlignmentPairs_to_GRangesList1 ################################################### U3grl <- as(U3galp, "GRangesList") U3grl ################################################### ### code chunk number 77: from_GAlignmentPairs_to_GRangesList2 ################################################### U3grl[ngap(U3galp) != 0] ################################################### ### code chunk number 78: coverage1 ################################################### U1cvg <- coverage(U1grl) U1cvg ################################################### ### code chunk number 79: coverage2 ################################################### mean(U1cvg) max(U1cvg) ################################################### ### code chunk number 80: slice1 ################################################### U1sl <- slice(U1cvg, lower=10) U1sl elementLengths(U1sl) head(U1sl$chr4) head(mean(U1sl$chr4)) head(max(U1sl$chr4)) ################################################### ### code chunk number 81: exbytx ################################################### library(TxDb.Dmelanogaster.UCSC.dm3.ensGene) txdb <- TxDb.Dmelanogaster.UCSC.dm3.ensGene exbytx <- exonsBy(txdb, by="tx", use.names=TRUE) exbytx ################################################### ### code chunk number 82: single-end-overlaps ################################################### U1txhits <- countOverlaps(exbytx, U1grl) length(U1txhits) head(U1txhits) sum(U1txhits) # total nb of hits head(sort(U1txhits, decreasing=TRUE)) ################################################### ### code chunk number 83: paired-end-overlaps ################################################### U3txhits <- countOverlaps(exbytx, U3grl) length(U3txhits) head(U3txhits) sum(U3txhits) # total nb of hits head(sort(U3txhits, decreasing=TRUE))