### R code from vignette source 'S4ContainersForHTS_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: Rle_example1 ################################################### library(IRanges) set.seed(2012) rle1 <- Rle(sample(c(-0.9, 0), 20, replace=TRUE)) rle1 runLength(rle1) runValue(rle1) as.vector(rle1) rle1[c(TRUE, FALSE)] ################################################### ### code chunk number 3: Rle_example2 ################################################### sort(rle1) rle1 * 50.1 sum(rle1) cumsum(rle1) cumsum(rle1) <= -4.2 rle1[cumsum(rle1) <= -4.2] ################################################### ### code chunk number 4: Rle_example3 ################################################### rle2 <- Rle(c("ch1", "chMT", "ch1", "ch2", "chMT"), c(4, 2, 1, 5, 1)) rle2 as.vector(rle2) c(rle2, c("chMT", "chX")) ################################################### ### code chunk number 5: Rle_example4 ################################################### runValue(rle2) <- factor(runValue(rle2)) rle2 runValue(rle2) as.vector(rle2) as.factor(rle2) ################################################### ### code chunk number 6: Rle_example5 ################################################### rle1 == 0 as(rle1 == 0, "IRanges") ################################################### ### code chunk number 7: IRanges_example1 ################################################### 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 8: IRanges_example2 ################################################### ir1[-2] ir2 <- c(ir1, IRanges(-10, 0)) ir2 ################################################### ### code chunk number 9: IRanges_example3 ################################################### duplicated(ir2) unique(ir2) sort(ir2) ################################################### ### code chunk number 10: 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 11: 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 12: 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 13: 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 14: IRanges_example4 ################################################### ir1 ################################################### ### code chunk number 15: IRanges_example5 ################################################### shift(ir1, -start(ir1)) flank(ir1, 10, start=FALSE) ################################################### ### code chunk number 16: IRanges_example6 ################################################### ir1 ################################################### ### code chunk number 17: IRanges_example7 ################################################### range(ir1) reduce(ir1) ################################################### ### code chunk number 18: IRanges_example8 ################################################### union(ir1, IRanges(-2, 6)) intersect(ir1, IRanges(-2, 13)) setdiff(ir1, IRanges(-2, 13)) ################################################### ### code chunk number 19: IRanges_example9 ################################################### ir3 <- IRanges(5:1, width=12) ir3 ################################################### ### code chunk number 20: IRanges_example10 ################################################### ir2 ################################################### ### code chunk number 21: IRanges_example11 ################################################### pintersect(ir3, ir2, resolve.empty="max.start") ################################################### ### code chunk number 22: IRanges_example12 ################################################### ok <- c(FALSE, FALSE, TRUE, TRUE, TRUE, FALSE, FALSE, TRUE) ir4 <- as(ok, "IRanges") # from logical vector to IRanges ir4 as(which(ok), "IRanges") # from integer vector to IRanges rle2[ir4] # IRanges subscript ################################################### ### code chunk number 23: DataFrame_example1 ################################################### library(Biostrings) dna <- DNAStringSet(c("AAA", "TGGATT", "CATTNGAGC", "TAATAG")) af <- alphabetFrequency(dna, baseOnly=TRUE) df <- DataFrame(dna, af) df df$G ################################################### ### code chunk number 24: DataFrame_example2 ################################################### df$cds_id <- paste("CDS", 1:4, sep="") df$cds_range <- successiveIRanges(width(dna), from=51) df as.data.frame(df) ################################################### ### code chunk number 25: CharacterList_example1 ################################################### ccl <- CharacterList(one=c("aaa", "bb", "c"), two=c("dd", "e", "fff", "gggg")) ################################################### ### code chunk number 26: CharacterList_example2 ################################################### ccl length(ccl) names(ccl) ################################################### ### code chunk number 27: CharacterList_example3 ################################################### as.list(ccl) ccl[[2]] ################################################### ### code chunk number 28: CharacterList_example4 ################################################### toupper(ccl) elementLengths(ccl) # fast version of sapply(ccl, length) unlist(ccl) # insane! will be changed soon... unlist(ccl, use.names=FALSE) ################################################### ### code chunk number 29: IntegerList_example1 ################################################### cil <- IntegerList(6:-2, 5, integer(0), 14:21) cil cil * cil ################################################### ### code chunk number 30: IntegerList_example2 ################################################### cil * 100L - 2L relist(unlist(cil) * 100L - 2L, cil) ################################################### ### code chunk number 31: IntegerList_example3 ################################################### cumsum(cil) ################################################### ### code chunk number 32: GRanges_constructor ################################################### library(GenomicRanges) gr1 <- GRanges(seqnames=rep(c("ch1", "chMT"), c(2, 4)), ranges=IRanges(start=16:21, end=20), strand=rep(c("+", "-", "*"), 2)) gr1 ################################################### ### code chunk number 33: GRanges_accessors1 ################################################### length(gr1) seqnames(gr1) ranges(gr1) ################################################### ### code chunk number 34: GRanges_accessors2 ################################################### start(gr1) end(gr1) width(gr1) strand(gr1) strand(gr1) <- c("-", "-", "+") strand(gr1) ################################################### ### code chunk number 35: GRanges_accessors3 ################################################### names(gr1) <- LETTERS[1:6] names(gr1) elementMetadata(gr1) <- DataFrame(score=11:16, GC=seq(1, 0, length=6)) elementMetadata(gr1) gr1 ################################################### ### code chunk number 36: GRanges_accessors4 ################################################### seqinfo(gr1) seqlevels(gr1) seqlengths(gr1) seqlengths(gr1) <- c(50000, 800) seqlengths(gr1) ################################################### ### code chunk number 37: GRanges_Vector_ops1 ################################################### gr1[c("F", "A")] gr1[strand(gr1) == "+"] ################################################### ### code chunk number 38: GRanges_Vector_ops2 ################################################### gr1 <- gr1[-5] gr1 ################################################### ### code chunk number 39: 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 40: GRanges_Vector_ops4 ################################################### gr12[length(gr12)] == gr12 duplicated(gr12) unique(gr12) ################################################### ### code chunk number 41: GRanges_Vector_ops5 ################################################### sort(gr12) ################################################### ### code chunk number 42: GRanges_Ranges_ops1 ################################################### gr2 shift(gr2, 50) narrow(gr2, start=2, end=-2) ################################################### ### code chunk number 43: GRanges_Ranges_ops2 ################################################### gr1 resize(gr1, 12) ################################################### ### code chunk number 44: GRanges_Ranges_ops3 ################################################### gr1 flank(gr1, 3) ################################################### ### code chunk number 45: GRanges_Ranges_ops4 ################################################### gr3 <- shift(gr1, c(35000, rep(0, 3), 100)) width(gr3)[c(3,5)] <- 117 gr3 range(gr3) ################################################### ### code chunk number 46: GRanges_Ranges_ops5 ################################################### gr3 disjoin(gr3) ################################################### ### code chunk number 47: GRanges_Ranges_ops6 ################################################### gr3 reduce(gr3) ################################################### ### code chunk number 48: GRanges_Ranges_ops7 ################################################### gr3 gaps(gr3) ################################################### ### code chunk number 49: GRanges_split ################################################### split(gr3, seqnames(gr3)) ################################################### ### code chunk number 50: GRangesList_constructor ################################################### grl <- GRangesList(gr3, gr2) grl ################################################### ### code chunk number 51: GRangesList_accessors1 ################################################### length(grl) ################################################### ### code chunk number 52: GRangesList_accessors2 ################################################### seqnames(grl) ################################################### ### code chunk number 53: GRangesList_accessors3 ################################################### strand(grl) ################################################### ### code chunk number 54: GRangesList_accessors4 ################################################### ranges(grl) ################################################### ### code chunk number 55: GRangesList_accessors5 ################################################### start(grl) end(grl) width(grl) ################################################### ### code chunk number 56: GRangesList_accessors6 ################################################### names(grl) <- c("TX1", "TX2") grl ################################################### ### code chunk number 57: GRangesList_accessors7 ################################################### elementMetadata(grl)$geneid <- c("GENE1", "GENE2") elementMetadata(grl) grl ################################################### ### code chunk number 58: GRangesList_accessors8 ################################################### seqinfo(grl) ################################################### ### code chunk number 59: GRangesList_Vector_ops1 ################################################### grl[c("TX2", "TX1")] ################################################### ### code chunk number 60: GRangesList_Vector_ops2 ################################################### c(grl, GRangesList(gr3)) ################################################### ### code chunk number 61: GRangesList_List_ops1 ################################################### grl[[2]] elementLengths(grl) unlisted <- unlist(grl, use.names=FALSE) # same as c(grl[[1]], grl[[2]]) unlisted ################################################### ### code chunk number 62: GRangesList_List_ops2 ################################################### grl100 <- relist(shift(unlisted, 100), grl) grl100 ################################################### ### code chunk number 63: GRangesList_List_ops3 ################################################### grl100b <- endoapply(grl, shift, 100) grl100b elementMetadata(grl100) elementMetadata(grl100b) ################################################### ### code chunk number 64: GRangesList_Ranges_ops1 ################################################### grl ################################################### ### code chunk number 65: GRangesList_Ranges_ops2 ################################################### shift(grl, 100) ################################################### ### code chunk number 66: GRangesList_Ranges_ops3 ################################################### grl ################################################### ### code chunk number 67: GRangesList_Ranges_ops4 ################################################### flank(grl, 10) ################################################### ### code chunk number 68: GRangesList_Ranges_ops5 ################################################### grl ################################################### ### code chunk number 69: GRangesList_Ranges_ops6 ################################################### range(grl) ################################################### ### code chunk number 70: GRangesList_Ranges_ops7 ################################################### grl ################################################### ### code chunk number 71: GRangesList_Ranges_ops8 ################################################### reduce(grl) ################################################### ### code chunk number 72: 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 73: GRangesList_Ranges_ops10 ################################################### grl2 grl3 ################################################### ### code chunk number 74: GRangesList_Ranges_ops11 ################################################### psetdiff(grl2, grl3) ################################################### ### code chunk number 75: GappedAlignments_constructor ################################################### gal0 <- GappedAlignments(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 76: readGappedAlignments ################################################### library(pasillaBamSubset) U1gal <- readGappedAlignments(untreated1_chr4()) length(U1gal) head(U1gal) ################################################### ### code chunk number 77: GappedAlignments_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 78: GappedAlignments_accessors2 ################################################### head(start(U1gal)) head(end(U1gal)) head(width(U1gal)) head(ngap(U1gal)) table(ngap(U1gal)) ################################################### ### code chunk number 79: GappedAlignments_accessors3 ################################################### elementMetadata(U1gal) seqinfo(U1gal) ################################################### ### code chunk number 80: readGappedAlignments2 ################################################### param <- ScanBamParam(what=c("flag", "mapq"), tag=c("NH", "NM")) U1gal <- readGappedAlignments(untreated1_chr4(), use.names=TRUE, param=param) U1gal[1:5] any(duplicated(names(U1gal))) ################################################### ### code chunk number 81: from_GappedAlignments_to_GRanges ################################################### as(U1gal, "GRanges") ################################################### ### code chunk number 82: from_GappedAlignments_to_GRangesList1 ################################################### U1grl <- as(U1gal, "GRangesList") U1grl ################################################### ### code chunk number 83: from_GappedAlignments_to_GRangesList2 ################################################### all(elementLengths(U1grl) == ngap(U1gal) + 1) ################################################### ### code chunk number 84: readGappedAlignmentPairs ################################################### library(pasillaBamSubset) U3galp <- readGappedAlignmentPairs(untreated3_chr4()) length(U3galp) head(U3galp) ################################################### ### code chunk number 85: GappedAlignmentPairs_accessors1 ################################################### head(first(U3galp)) head(last(U3galp)) ################################################### ### code chunk number 86: GappedAlignmentPairs_accessors2 ################################################### seqnames(U3galp) strand(U3galp) head(ngap(U3galp)) table(ngap(U3galp)) ################################################### ### code chunk number 87: from_GappedAlignmentPairs_to_GRangesList1 ################################################### U3grl <- as(U3galp, "GRangesList") U3grl ################################################### ### code chunk number 88: from_GappedAlignmentPairs_to_GRangesList2 ################################################### U3grl[ngap(U3galp) != 0] ################################################### ### code chunk number 89: coverage1 ################################################### U1cvg <- coverage(U1grl) ################################################### ### code chunk number 90: coverage2 ################################################### mean(U1cvg) max(U1cvg) ################################################### ### code chunk number 91: slice1 ################################################### U1sl <- slice(U1cvg, lower=10) U1sl elementLengths(U1sl) head(U1sl$chr4) head(mean(U1sl$chr4)) head(max(U1sl$chr4)) ################################################### ### code chunk number 92: 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 93: 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 94: paired-end-overlaps ################################################### U3txhits <- countOverlaps(exbytx, U3grl) length(U3txhits) head(U3txhits) sum(U3txhits) # total nb of hits head(sort(U3txhits, decreasing=TRUE))