### R code from vignette source 'vignettes/ggbio/inst/doc/intro.Rnw' ################################################### ### code chunk number 1: options ################################################### options(width=72) ################################################### ### code chunk number 2: load ################################################### library(ggbio) ################################################### ### code chunk number 3: simul ################################################### set.seed(1) N <- 1000 library(GenomicRanges) gr <- GRanges(seqnames = sample(c("chr1", "chr2", "chr3"), size = N, replace = TRUE), IRanges( start = sample(1:300, size = N, replace = TRUE), width = sample(70:75, size = N,replace = TRUE)), strand = sample(c("+", "-", "*"), size = N, replace = TRUE), value = rnorm(N, 10, 3), score = rnorm(N, 100, 30), sample = sample(c("Normal", "Tumor"), size = N, replace = TRUE), pair = sample(letters, size = N, replace = TRUE)) idx <- sample(1:length(gr), size = 50) ################################################### ### code chunk number 4: default ################################################### autoplot(gr[idx]) ################################################### ### code chunk number 5: bar-default-pre ################################################### set.seed(123) gr.b <- GRanges(seqnames = "chr1", IRanges(start = seq(1, 100, by = 10), width = sample(4:9, size = 10, replace = TRUE)), score = rnorm(10, 10, 3), value = runif(10, 1, 100)) gr.b2 <- GRanges(seqnames = "chr2", IRanges(start = seq(1, 100, by = 10), width = sample(4:9, size = 10, replace = TRUE)), score = rnorm(10, 10, 3), value = runif(10, 1, 100)) gr.b <- c(gr.b, gr.b2) head(gr.b) ################################################### ### code chunk number 6: bar-default ################################################### p1 <- autoplot(gr.b, geom = "bar") ## use value to fill the bar p2 <- autoplot(gr.b, geom = "bar", aes(fill = value)) tracks(default = p1, fill = p2) ################################################### ### code chunk number 7: intro.Rnw:393-394 ################################################### autoplot(gr[idx], geom = "arch", aes(color = value), facets = sample ~ seqnames) ################################################### ### code chunk number 8: gr-group ################################################### gra <- GRanges("chr1", IRanges(c(1,7,20), end = c(4,9,30)), group = c("a", "a", "b")) ## if you desn't specify group, then group based on stepping levels, and gaps are computed without ## considering extra group method p1 <- autoplot(gra, aes(fill = group), geom = "alignment") ## when use group method, gaps only computed for grouped intervals. ## default is group.selfish = TRUE, each group keep one row. ## in this way, group labels could be shown as y axis. p2 <- autoplot(gra, aes(fill = group, group = group), geom = "alignment") ## group.selfish = FALSE, save space p3 <- autoplot(gra, aes(fill = group, group = group), geom = "alignment", group.selfish = FALSE) tracks('non-group' = p1,'group.selfish = TRUE' = p2 , 'group.selfish = FALSE' = p3) ################################################### ### code chunk number 9: gr-facet-strand ################################################### autoplot(gr, stat = "coverage", geom = "area", facets = strand ~ seqnames, aes(fill = strand)) ################################################### ### code chunk number 10: gr-autoplot-circle ################################################### autoplot(gr[idx], layout = 'circle') ################################################### ### code chunk number 11: gr-circle ################################################### seqlengths(gr) <- c(400, 500, 700) values(gr)$to.gr <- gr[sample(1:length(gr), size = length(gr))] idx <- sample(1:length(gr), size = 50) gr <- gr[idx] ggplot() + layout_circle(gr, geom = "ideo", fill = "gray70", radius = 7, trackWidth = 3) + layout_circle(gr, geom = "bar", radius = 10, trackWidth = 4, aes(fill = score, y = score)) + layout_circle(gr, geom = "point", color = "red", radius = 14, trackWidth = 3, grid = TRUE, aes(y = score)) + layout_circle(gr, geom = "link", linked.to = "to.gr", radius = 6, trackWidth = 1) ################################################### ### code chunk number 12: seqinfo-src ################################################### data(hg19Ideogram, package = "biovizBase") sq <- seqinfo(hg19Ideogram) sq ################################################### ### code chunk number 13: seqinfo ################################################### autoplot(sq[paste0("chr", c(1:22, "X"))]) ################################################### ### code chunk number 14: ir-load ################################################### set.seed(1) N <- 100 ir <- IRanges(start = sample(1:300, size = N, replace = TRUE), width = sample(70:75, size = N,replace = TRUE)) ## add meta data df <- DataFrame(value = rnorm(N, 10, 3), score = rnorm(N, 100, 30), sample = sample(c("Normal", "Tumor"), size = N, replace = TRUE), pair = sample(letters, size = N, replace = TRUE)) values(ir) <- df ir ################################################### ### code chunk number 15: ir-exp ################################################### p1 <- autoplot(ir) p2 <- autoplot(ir, aes(fill = pair)) + theme(legend.position = "none") p3 <- autoplot(ir, stat = "coverage", geom = "line", facets = sample ~. ) p4 <- autoplot(ir, stat = "reduce") tracks(p1, p2, p3, p4) ################################################### ### code chunk number 16: grl-simul ################################################### set.seed(1) N <- 100 ## ====================================================================== ## simmulated GRanges ## ====================================================================== gr <- GRanges(seqnames = sample(c("chr1", "chr2", "chr3"), size = N, replace = TRUE), IRanges( start = sample(1:300, size = N, replace = TRUE), width = sample(30:40, size = N,replace = TRUE)), strand = sample(c("+", "-", "*"), size = N, replace = TRUE), value = rnorm(N, 10, 3), score = rnorm(N, 100, 30), sample = sample(c("Normal", "Tumor"), size = N, replace = TRUE), pair = sample(letters, size = N, replace = TRUE)) grl <- split(gr, values(gr)$pair) ################################################### ### code chunk number 17: grl-exp ################################################### ## default gap.geom is 'chevron' p1 <- autoplot(grl, group.selfish = TRUE) p2 <- autoplot(grl, group.selfish = TRUE, main.geom = "arrowrect", gap.geom = "segment") tracks(p1, p2) ################################################### ### code chunk number 18: grl-name ################################################### autoplot(grl, aes(fill = ..grl_name..)) ## equal to ## autoplot(grl, aes(fill = grl_name)) ################################################### ### code chunk number 19: rle-simul ################################################### library(IRanges) library(ggbio) set.seed(1) lambda <- c(rep(0.001, 4500), seq(0.001, 10, length = 500), seq(10, 0.001, length = 500)) ## @knitr create xVector <- rpois(1e4, lambda) xRle <- Rle(xVector) xRle ################################################### ### code chunk number 20: rle-bin ################################################### p1 <- autoplot(xRle) p2 <- autoplot(xRle, nbin = 80) p3 <- autoplot(xRle, geom = "heatmap", nbin = 200) tracks('nbin = 30' = p1, "nbin = 80" = p2, "nbin = 200(heatmap)" = p3) ################################################### ### code chunk number 21: rle-id ################################################### p1 <- autoplot(xRle, stat = "identity") p2 <- autoplot(xRle, stat = "identity", geom = "point", color = "red") tracks('line' = p1, "point" = p2) ################################################### ### code chunk number 22: rle-slice ################################################### p1 <- autoplot(xRle, type = "viewMaxs", stat = "slice", lower = 5) p2 <- autoplot(xRle, type = "viewMaxs", stat = "slice", lower = 5, geom = "heatmap") tracks('bar' = p1, "heatmap" = p2) ################################################### ### code chunk number 23: rlel-simul ################################################### xRleList <- RleList(xRle, 2L * xRle) xRleList ################################################### ### code chunk number 24: rlel-bin ################################################### p1 <- autoplot(xRleList) p2 <- autoplot(xRleList, nbin = 80) p3 <- autoplot(xRleList, geom = "heatmap", nbin = 200) tracks('nbin = 30' = p1, "nbin = 80" = p2, "nbin = 200(heatmap)" = p3) ################################################### ### code chunk number 25: rlel-id ################################################### p1 <- autoplot(xRleList, stat = "identity") p2 <- autoplot(xRleList, stat = "identity", geom = "point", color = "red") tracks('line' = p1, "point" = p2) ################################################### ### code chunk number 26: rlel-slice ################################################### p1 <- autoplot(xRleList, type = "viewMaxs", stat = "slice", lower = 5) p2 <- autoplot(xRleList, type = "viewMaxs", stat = "slice", lower = 5, geom = "heatmap") tracks('bar' = p1, "heatmap" = p2) ################################################### ### code chunk number 27: txdb ################################################### library(TxDb.Hsapiens.UCSC.hg19.knownGene) data(genesymbol, package = "biovizBase") txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene ################################################### ### code chunk number 28: txdb-visual ################################################### p1 <- autoplot(txdb, which = genesymbol["ALDOA"], names.expr = "tx_name:::gene_id") p2 <- autoplot(txdb, which = genesymbol["ALDOA"], stat = "reduce", color = "brown", fill = "brown") tracks(full = p1, reduce = p2, heights = c(5, 1)) + ylab("") ################################################### ### code chunk number 29: ga-load ################################################### library(Rsamtools) data("genesymbol", package = "biovizBase") bamfile <- system.file("extdata", "SRR027894subRBM17.bam", package="biovizBase") ## need to set use.names = TRUE ga <- readBamGappedAlignments(bamfile, param = ScanBamParam(which = genesymbol["RBM17"]), use.names = TRUE) ################################################### ### code chunk number 30: ga-exp ################################################### p1 <- autoplot(ga) p2 <- autoplot(ga, geom = "rect") p3 <- autoplot(ga, geom = "line", stat = "coverage") tracks(default = p1, rect = p2, coverage = p3) ################################################### ### code chunk number 31: bf-load (eval = FALSE) ################################################### ## library(Rsamtools) ## bamfile <- "./wgEncodeCaltechRnaSeqK562R1x75dAlignsRep1V2.bam" ## bf <- BamFile(bamfile) ################################################### ### code chunk number 32: bf-est-cov (eval = FALSE) ################################################### ## autoplot(bamfile) ## autoplot(bamfile, which = c("chr1", "chr2")) ## autoplot(bf) ## autoplot(bf, which = c("chr1", "chr2")) ## ## data(genesymbol, package = "biovizBase") ## autoplot(bamfile, method = "raw", which = genesymbol["ALDOA"]) ## ## library(BSgenome.Hsapiens.UCSC.hg19) ## autoplot(bf, stat = "mismatch", which = genesymbol["ALDOA"], bsgenome = Hsapiens) ################################################### ### code chunk number 33: char-bam (eval = FALSE) ################################################### ## bamfile <- "./wgEncodeCaltechRnaSeqK562R1x75dAlignsRep1V2.bam" ## autoplot(bamfile) ################################################### ### code chunk number 34: char-gr ################################################### library(rtracklayer) test_path <- system.file("tests", package = "rtracklayer") test_bed <- file.path(test_path, "test.bed") autoplot(test_bed, aes(fill = name)) ################################################### ### code chunk number 35: matrix-default ################################################### autoplot(volcano) ################################################### ### code chunk number 36: matrix-default-scale ################################################### autoplot(volcano-150)+scale_fill_fold_change() ################################################### ### code chunk number 37: matrix-default-gene ################################################### colnames(volcano) <- sort(sample(1:300, size = ncol(volcano), replace = FALSE)) autoplot(volcano-150, genomic.pos = TRUE)+scale_fill_fold_change() ################################################### ### code chunk number 38: intro.Rnw:989-995 ################################################### library(Biobase) data(sample.ExpressionSet) sample.ExpressionSet set.seed(1) idx <- sample(seq_len(dim(sample.ExpressionSet)[1]), size = 50) eset <- sample.ExpressionSet[idx,] ################################################### ### code chunk number 39: eset-default ################################################### p1 <- autoplot(eset) p1 ################################################### ### code chunk number 40: eset-default-scale ################################################### p2 <- p1 + scale_fill_fold_change() p2 ################################################### ### code chunk number 41: eset-default-pcp ################################################### autoplot(eset, type = "pcp") ################################################### ### code chunk number 42: eset-default-boxplot ################################################### autoplot(eset, type = "boxplot") ################################################### ### code chunk number 43: eset-default-sm ################################################### autoplot(eset[, 1:7], type = "scatterplot.matrix") ################################################### ### code chunk number 44: eset-default-ms ################################################### autoplot(eset, type = "mean-sd") ################################################### ### code chunk number 45: eset-default-volcano ################################################### autoplot(eset, type = "volcano", fac = pData(sample.ExpressionSet)$type) ################################################### ### code chunk number 46: sset ################################################### nrows <- 200; ncols <- 6 counts <- matrix(runif(nrows * ncols, 1, 1e4), nrows) rowData <- GRanges(rep(c("chr1", "chr2"), c(50, 150)), IRanges(floor(runif(200, 1e5, 1e6)), width=100), strand=sample(c("+", "-"), 200, TRUE)) colData <- DataFrame(Treatment=rep(c("ChIP", "Input"), 3), row.names=LETTERS[1:6]) sset <- SummarizedExperiment(assays=SimpleList(counts=counts), rowData=rowData, colData=colData) ################################################### ### code chunk number 47: sset-heatmap ################################################### autoplot(sset) + scale_fill_fold_change() ################################################### ### code chunk number 48: sset-pcp ################################################### autoplot(sset, type = "pcp") ################################################### ### code chunk number 49: sset-boxplot ################################################### autoplot(sset, type = "boxplot") ################################################### ### code chunk number 50: sset-sm ################################################### autoplot(sset, type = "scatterplot.matrix") ################################################### ### code chunk number 51: load ################################################### library(biovizBase) data(hg19IdeogramCyto) ## data structure hg19IdeogramCyto ## return TRUE, if the object could be visualized by ggbio biovizBase::isIdeogram(hg19IdeogramCyto) ################################################### ### code chunk number 52: ideo-ori ################################################### library(ggbio) p <- plotIdeogram(hg19IdeogramCyto, "chr1") print(p) ## big, need to be resized. ################################################### ### code chunk number 53: ideo-ori-viewport ################################################### ## to use function viewport, load grid library(grid) print(p, vp = viewport(height = 0.15, width = 1)) ################################################### ### code chunk number 54: ideo-ori-zoom ################################################### plotIdeogram(hg19IdeogramCyto, "chr1", zoom.region = c(1e8, 1.5e8)) ################################################### ### code chunk number 55: ideo-ori-xlabel ################################################### plotIdeogram(hg19IdeogramCyto, "chr1", xlabel = TRUE) ################################################### ### code chunk number 56: ideo-ori-nocyto ################################################### library(GenomicRanges) ## there are no seqlengths seqlengths(hg19IdeogramCyto) ## so directly plot will try to aggregate and estimate lengths of chromosomes, ## this is not accurate data(hg19IdeogramCyto) p1 <- plotIdeogram(hg19IdeogramCyto, "chr1", cytoband = FALSE, xlabel = TRUE) ## let's assign a short length to this object hg19_fake_chr1 <- hg19IdeogramCyto seqlengths(hg19_fake_chr1)[1] <- 1e8 ## this will use it's "seqlengths" information to visualize the chromosome. p2 <- plotIdeogram(hg19_fake_chr1, "chr1", cytoband = FALSE, xlabel = TRUE) ## see the difference tracks(p1, p2) ################################################### ### code chunk number 57: hg19-nocyto ################################################### data(hg19Ideogram) head(hg19Ideogram) ################################################### ### code chunk number 58: getIdeogram (eval = FALSE) ################################################### ## library(rtracklayer) ## ## need UCSC connection ## head(ucscGenomes()) ################################################### ### code chunk number 59: getideo-no (eval = FALSE) ################################################### ## library(biovizBase) ## obj <- getIdeogram() ################################################### ### code chunk number 60: mm9 ################################################### library(biovizBase) ## just need information about chromosome lengths mm9 <- getIdeogram("mm9", cytoband = FALSE) ## have head(mm9) ## need information about cytoband mm9 <- getIdeogram("mm9") head(mm9) ## with extra column 'name' and 'gieStain'. ################################################### ### code chunk number 61: cytoband ################################################### cyto.def <- getOption("biovizBase")$cytobandColor cyto.def setdiff(unique(values(mm9)$gieStain), names(cyto.def)) ################################################### ### code chunk number 62: cyto-color ################################################### p1 <- plotIdeogram(mm9, "chr1") cyto.def cyto.new <- c(cyto.def, c(gpos33 = "grey80", gpos66 = "grey60")) ## method 1: optlist <- getOption("biovizBase") optlist$cytobandColor <- cyto.new options(biovizBase = optlist) p2 <- plotIdeogram(mm9, "chr1") p3 <- plotIdeogram(mm9, "chr1") + scale_fill_manual(values = cyto.new) tracks(p1, p2, p3) ################################################### ### code chunk number 63: seqinfo ################################################### data(hg19Ideogram) seqs <- seqinfo(hg19Ideogram) class(seqs) p1 <- autoplot(seqs["chr1"]) p2 <- autoplot(seqs["chr1"], FALSE) tracks(type1 = p1, type2 = p2) ################################################### ### code chunk number 64: gr-simul ################################################### set.seed(1) library(GenomicRanges) ## simulate short reads gr1 <- GRanges("chr1", IRanges(start = sample(1:10, size = 50, replace = TRUE), width = 4)) ## simulate exons grl <- GRangesList(GRanges("chr1", IRanges(start = c(1, 10, 20), width = 5)), GRanges("chr1", IRanges(start = c(1, 20), width = 5))) ################################################### ### code chunk number 65: shortread-plot ################################################### library(ggbio) p1 <- autoplot(gr1) print(p1) ################################################### ### code chunk number 66: exons-plot ################################################### p2 <- autoplot(grl) print(p2) ################################################### ### code chunk number 67: grid_arrange ################################################### library(gridExtra) grid.arrange(p1, p2) ################################################### ### code chunk number 68: tracks-sr-ex ################################################### tracks(p1, p2) ################################################### ### code chunk number 69: tracks-sr-ex-more ################################################### tracks("Short Reads" = p1, "Exons" = p2, heights = c(3, 1)) ## this is equivalent to ## tracks(list("Short Reads" = p1, "Exons" = p2)) ################################################### ### code chunk number 70: tracks-sr-ex-more2 ################################################### ## to use unit(), load grid. library(grid) tracks("Short Reads" = p1, "Exons" = p2, heights = c(3, 1), label.text.cex = 2, label.text.color = "white", label.bg.fill = "brown", label.width = unit(5, "line")) ################################################### ### code chunk number 71: plus ################################################### ## change track plot background color p.track <- tracks("Short Reads" = p1, "Exons" = p2, heights = c(4, 1), label.text.color = "white", label.bg.fill = "brown", track.plot.color = c("white", "#FFFEDB")) ## apply theme_null() to all track p.track + theme_null() ################################################### ### code chunk number 72: specific-plus ################################################### ## only apply change to one track p.track@grobs[[1]] <- p.track@grobs[[1]] + theme_null() p.track ################################################### ### code chunk number 73: zoom-in (eval = FALSE) ################################################### ## p.track <- tracks("Short Reads" = p1, "Exons" = p2, heights = c(4, 1)) ## ## ## remove data outside c(1, 6) ## p.track + xlim(c(1, 6)) ## ## ## keep original data, simply zooming to [1, 10] ## p.track + coord_cartesian(xlim = c(1, 10)) ## ## ## more easy with ggbio's 'xlim<-', you need to print it to see the effect. ## xlim(p.track) <- c(1, 6) ## p.track ## ## ## or update, if a graphic device is open, you can see the effect, immediately. ## update(p.track, xlim = c(1, 4)) ################################################### ### code chunk number 74: backup-and-reset ################################################### p.track xlim(p.track) <- c(1, 4) p.track ## back it up p.track <- backup(p.track) ## change to another view xlim(p.track) <- c(1, 10) p.track reset(p.track) ################################################### ### code chunk number 75: ggplot ################################################### p.gg <- qplot(x = seq(from = 1, to = 20, length.out = 500), y = rnorm(500), geom = "area") p.track <- tracks("density" = p.gg, "Short Reads" = p1, "Exons" = p2, heights = c(2, 5, 1)) p.track ################################################### ### code chunk number 76: ggplot-facet ################################################### p.gg <- qplot(data = mtcars, x = mpg, y = wt, facets = cyl ~.) p.track <- tracks("Random" = p.gg, "Short Reads" = p1, "Exons" = p2, heights = c(4, 5, 2)) p.track ################################################### ### code chunk number 77: loading-txdb ################################################### library(TxDb.Hsapiens.UCSC.hg19.knownGene) txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene ## suppose you already know the region you want to visualize ## or for human genome, you can try following commented code ## data(genesymbol, package = "biovizBase") ## genesymbol["ALDOA"] aldoa.gr <- GRanges("chr16", IRanges(30064491, 30081734)) aldoa.gr ################################################### ### code chunk number 78: txdb-full ################################################### library(ggbio) p1 <- autoplot(txdb, which = aldoa.gr) print(p1) ################################################### ### code chunk number 79: txdb-full-aes ################################################### library(ggbio) p1 <- autoplot(txdb, which = aldoa.gr, fill = "brown", color = "brown") print(p1) ################################################### ### code chunk number 80: txdb-reduce ################################################### p2 <- autoplot(txdb, which = aldoa.gr, stat = "reduce") print(p2) ################################################### ### code chunk number 81: tracks ################################################### tracks(full = p1, reduced = p2, heights = c(4,1)) + theme_alignment(grid=FALSE, border = FALSE) ################################################### ### code chunk number 82: change-intron-geom ################################################### autoplot(txdb, which = aldoa.gr, gap.geom = "chevron") ################################################### ### code chunk number 83: change-intron-geom-arrow ################################################### library(grid) autoplot(txdb, which = aldoa.gr, arrow.rate = 0.001, length = unit(0.35, "cm")) ################################################### ### code chunk number 84: parsing-expression ################################################### p <- autoplot(txdb, which = aldoa.gr, names.expr = "gene_id:::tx_name") p ################################################### ### code chunk number 85: change-unit ################################################### p + scale_x_sequnit("kb") ################################################### ### code chunk number 86: stat_gene (eval = FALSE) ################################################### ## p1 <- ggplot() + stat_gene(txdb, which = aldoa.gr) ################################################### ### code chunk number 87: loading ################################################### library(ggbio) data(hg19IdeogramCyto, package = "biovizBase") head(hg19IdeogramCyto) ## default pre-set color stored in getOption("biovizBase")$cytobandColor ################################################### ### code chunk number 88: default ################################################### autoplot(hg19IdeogramCyto, layout = "karyogram", cytoband = TRUE) ################################################### ### code chunk number 89: change-order ################################################### library(GenomicRanges) hg19 <- keepSeqlevels(hg19IdeogramCyto, paste0("chr", c(1:22, "X", "Y"))) head(hg19) autoplot(hg19, layout = "karyogram", cytoband = TRUE) ################################################### ### code chunk number 90: cyto-normal ################################################### library(GenomicRanges) ## it's a 'ideogram' biovizBase::isIdeogram(hg19) ## set to FALSE autoplot(hg19, layout = "karyogram", cytoband = FALSE, aes(fill = gieStain)) + scale_fill_giemsa() ################################################### ### code chunk number 91: load-RNAediting ################################################### data(darned_hg19_subset500, package = "biovizBase") dn <- darned_hg19_subset500 head(dn) ## add seqlengths ## we have seqlegnths information in another data set data(hg19Ideogram, package = "biovizBase") seqlengths(dn) <- seqlengths(hg19Ideogram)[names(seqlengths(dn))] ## now we have seqlengths head(dn) ## then we change order dn <- keepSeqlevels(dn, paste0("chr", c(1:22, "X"))) autoplot(dn, layout = "karyogram") ## this equivalent to ## autoplot(seqinfo(dn)) ################################################### ### code chunk number 92: load-RNAediting-color ################################################### ## since default is geom rectangle, even though it's looks like segment ## we still use both fill/color to map colors autoplot(dn, layout = "karyogram", aes(color = exReg, fill = exReg)) ################################################### ### code chunk number 93: load-RNAediting-color-NA ################################################### ## since default is geom rectangle, even though it's looks like segment ## we still use both fill/color to map colors autoplot(dn, layout = "karyogram", aes(color = exReg, fill = exReg)) + scale_color_discrete(na.value = "brown") ################################################### ### code chunk number 94: load-RNAediting-color-fake ################################################### dn2 <- dn seqlengths(dn2) <- rep(max(seqlengths(dn2)), length(seqlengths(dn2)) ) autoplot(dn2, layout = "karyogram", aes(color = exReg, fill = exReg)) ################################################### ### code chunk number 95: plotKaryogram (eval = FALSE) ################################################### ## plotKaryogram(dn) ## plotKaryogram(dn, aes(color = exReg, fill = exReg)) ################################################### ### code chunk number 96: low-default ################################################### ## plot ideogram p <- ggplot() + layout_karyogram(hg19, cytoband = TRUE) p ## eqevelant autoplot(hg19, layout = "karyogram", cytoband = TRUE) ################################################### ### code chunk number 97: low-default-addon ################################################### p <- p + layout_karyogram(dn, geom = "rect", ylim = c(11, 21), color = "red") ## commented line below won't work ## the cytoband fill color has been used already. ## p <- p + layout_karyogram(dn, aes(fill = exReg, color = exReg), geom = "rect") p ################################################### ### code chunk number 98: edit-space ################################################### ## plot chromosome space p <- autoplot(seqinfo(dn)) ## make sure you pass rect as geom ## otherwise you just get background p <- p + layout_karyogram(dn, aes(fill = exReg, color = exReg), geom = "rect") values(dn)$pvalue <- rnorm(length(dn)) p + layout_karyogram(dn, aes(x = start, y = pvalue), ylim = c(10, 30), geom = "line", color = "red") p ################################################### ### code chunk number 99: processing ################################################### crc1 <- system.file("extdata", "crc1-missense.csv", package = "biovizBase") crc1 <- read.csv(crc1) library(GenomicRanges) mut.gr <- with(crc1,GRanges(Chromosome, IRanges(Start_position, End_position), strand = Strand)) values(mut.gr) <- subset(crc1, select = -c(Start_position, End_position, Chromosome)) data("hg19Ideogram", package = "biovizBase") seqs <- seqlengths(hg19Ideogram) ## subset_chr chr.sub <- paste("chr", 1:22, sep = "") ## levels tweak seqlevels(mut.gr) <- c(chr.sub, "chrX") mut.gr <- keepSeqlevels(mut.gr, chr.sub) seqs.sub <- seqs[chr.sub] ## remove wrong position bidx <- end(mut.gr) <= seqs.sub[match(as.character(seqnames(mut.gr)), names(seqs.sub))] mut.gr <- mut.gr[which(bidx)] ## assign_seqlengths seqlengths(mut.gr) <- seqs.sub ## reanme to shorter names new.names <- as.character(1:22) names(new.names) <- paste("chr", new.names, sep = "") new.names mut.gr.new <- renameSeqlevels(mut.gr, new.names) head(mut.gr.new) ################################################### ### code chunk number 100: ideo ################################################### hg19Ideo <- hg19Ideogram hg19Ideo <- keepSeqlevels(hg19Ideogram, chr.sub) hg19Ideo <- renameSeqlevels(hg19Ideo, new.names) head(hg19Ideo) ################################################### ### code chunk number 101: lower-ideo-track ################################################### library(ggbio) p <- ggplot() + layout_circle(hg19Ideo, geom = "ideo", fill = "gray70", radius = 30, trackWidth = 4) p ################################################### ### code chunk number 102: lower-scale-track ################################################### p <- p + layout_circle(hg19Ideo, geom = "scale", size = 2, radius = 35, trackWidth = 2) p ################################################### ### code chunk number 103: lower-text-track ################################################### p <- p + layout_circle(hg19Ideo, geom = "text", aes(label = seqnames), vjust = 0, radius = 38, trackWidth = 7) p ################################################### ### code chunk number 104: lower-mut-track ################################################### p <- p + layout_circle(mut.gr, geom = "rect", color = "steelblue", radius = 23 ,trackWidth = 6) p ################################################### ### code chunk number 105: links ################################################### rearr <- read.csv(system.file("extdata", "crc-rearrangment.csv", package = "biovizBase")) ## start position gr1 <- with(rearr, GRanges(chr1, IRanges(pos1, width = 1))) ## end position gr2 <- with(rearr, GRanges(chr2, IRanges(pos2, width = 1))) ## add extra column nms <- colnames(rearr) .extra.nms <- setdiff(nms, c("chr1", "chr2", "pos1", "pos2")) values(gr1) <- rearr[,.extra.nms] ## remove out-of-limits data seqs <- as.character(seqnames(gr1)) .mx <- seqlengths(hg19Ideo)[seqs] idx1 <- start(gr1) > .mx seqs <- as.character(seqnames(gr2)) .mx <- seqlengths(hg19Ideo)[seqs] idx2 <- start(gr2) > .mx idx <- !idx1 & !idx2 gr1 <- gr1[idx] seqlengths(gr1) <- seqlengths(hg19Ideo) gr2 <- gr2[idx] seqlengths(gr2) <- seqlengths(hg19Ideo) ################################################### ### code chunk number 106: link-data ################################################### values(gr1)$to.gr <- gr2 ## rename to gr gr <- gr1 ################################################### ### code chunk number 107: rearr ################################################### values(gr)$rearrangements <- ifelse(as.character(seqnames(gr)) == as.character(seqnames((values(gr)$to.gr))), "intrachromosomal", "interchromosomal") ################################################### ### code chunk number 108: subset-crc-1 ################################################### gr.crc1 <- gr[values(gr)$individual == "CRC-1"] ################################################### ### code chunk number 109: lower-point-track ################################################### p <- p + layout_circle(gr.crc1, geom = "point", aes(y = score, size = tumreads), color = "red", radius = 12 ,trackWidth = 10, grid = TRUE) + scale_size(range = c(1, 2.5)) p ################################################### ### code chunk number 110: lower-link-track ################################################### p <- p + layout_circle(gr.crc1, geom = "link", linked.to = "to.gr", aes(color = rearrangements), radius = 10 ,trackWidth = 1) p ################################################### ### code chunk number 111: single-arr ################################################### cols <- RColorBrewer::brewer.pal(3, "Set2")[2:1] names(cols) <- c("interchromosomal", "intrachromosomal") p0 <- ggplot() + layout_circle(gr.crc1, geom = "link", linked.to = "to.gr", aes(color = rearrangements), radius = 7.1) + layout_circle(hg19Ideo, geom = "ideo", trackWidth = 1.5, color = "gray70", fill = "gray70") + scale_color_manual(values = cols) p0 ################################################### ### code chunk number 112: legend ################################################### library(gridExtra) g <- ggplotGrob(p0) idx <- which(g$layout$name == "guide-box") ## gg <- editGrob(getGrob(g, gPath("guide-box"), ## grep=TRUE), vp=viewport()) gg <- g$grobs[[idx]] ################################################### ### code chunk number 113: arrangement ################################################### grl <- split(gr, values(gr)$individual) ## need "unit", load grid library(grid) lst <- lapply(grl, function(gr.cur){ print(unique(as.character(values(gr.cur)$individual))) cols <- RColorBrewer::brewer.pal(3, "Set2")[2:1] names(cols) <- c("interchromosomal", "intrachromosomal") p <- ggplot() + layout_circle(gr.cur, geom = "link", linked.to = "to.gr", aes(color = rearrangements), radius = 7.1) + layout_circle(hg19Ideo, geom = "ideo", trackWidth = 1.5, color = "gray70", fill = "gray70") + scale_color_manual(values = cols) + labs(title = (unique(values(gr.cur)$individual))) + theme(plot.margin = unit(rep(0, 4), "lines")) }) lst.nolegend <- lapply(lst, function(p) p + theme(legend.position = "none")) l.g <- lapply(lst.nolegend, ggplotGrob) ################################################### ### code chunk number 114: 9-circle ################################################### grid.arrange(do.call(arrangeGrob, l.g), gg, ncol = 2, widths = c(4/5, 1/5)) ################################################### ### code chunk number 115: simple-wrapper (eval = FALSE) ################################################### ## ## noitce this code doesn't have theme(legend.position = "none"), which means ## ## we keep the legend for each plot. ## arrangeGrobByParsingLegend(lst, widths = c(4, 1), legend.idx = 1, ncol = 2) ################################################### ### code chunk number 116: load ################################################### library(ggbio) data(hg19IdeogramCyto, package = "biovizBase") data(hg19Ideogram, package = "biovizBase") library(GenomicRanges) ################################################### ### code chunk number 117: simul_gr ################################################### library(biovizBase) gr <- GRanges(rep(c("chr1", "chr2"), each = 5), IRanges(start = rep(seq(1, 100, length = 5), times = 2), width = 50)) autoplot(gr, aes(fill = seqnames)) ################################################### ### code chunk number 118: coord-genome ################################################### autoplot(gr, coord = "genome", aes(fill = seqnames)) ################################################### ### code chunk number 119: is ################################################### gr.t <- transformToGenome(gr) head(gr.t) is_coord_genome(gr.t) metadata(gr.t)$coord ################################################### ### code chunk number 120: simul_snp ################################################### chrs <- as.character(levels(seqnames(hg19IdeogramCyto))) seqlths <- seqlengths(hg19Ideogram)[chrs] set.seed(1) nchr <- length(chrs) nsnps <- 100 gr.snp <- GRanges(rep(chrs,each=nsnps), IRanges(start = do.call(c, lapply(chrs, function(chr){ N <- seqlths[chr] runif(nsnps,1,N) })), width = 1), SNP=sapply(1:(nchr*nsnps), function(x) paste("rs",x,sep='')), pvalue = -log10(runif(nchr*nsnps)), group = sample(c("Normal", "Tumor"), size = nchr*nsnps, replace = TRUE) ) genome(gr.snp) <- "hg19" gr.snp ################################################### ### code chunk number 121: shorter ################################################### seqlengths(gr.snp) nms <- seqnames(seqinfo(gr.snp)) nms.new <- gsub("chr", "", nms) names(nms.new) <- nms gr.snp <- renameSeqlevels(gr.snp, nms.new) seqlengths(gr.snp) ################################################### ### code chunk number 122: unorder ################################################### autoplot(gr.snp, coord = "genome", geom = "point", aes(y = pvalue), space.skip = 0.01) ################################################### ### code chunk number 123: sort ################################################### gr.snp <- keepSeqlevels(gr.snp, c(1:22, "X", "Y")) autoplot(gr.snp, coord = "genome", geom = "point", aes(y = pvalue), space.skip = 0.01) ################################################### ### code chunk number 124: with_seql ################################################### names(seqlths) <- gsub("chr", "", names(seqlths)) seqlengths(gr.snp) <- seqlths[names(seqlengths(gr.snp))] autoplot(gr.snp, coord = "genome", geom = "point", aes(y = pvalue), space.skip = 0.01) ################################################### ### code chunk number 125: line ################################################### autoplot(gr.snp, coord = "genome", geom = "line", aes(y = pvalue, group = seqnames, color = seqnames)) ################################################### ### code chunk number 126: plotGrandLinear ################################################### plotGrandLinear(gr.snp, aes(y = pvalue)) ################################################### ### code chunk number 127: morecolor1 ################################################### plotGrandLinear(gr.snp, aes(y = pvalue, color = seqnames)) ################################################### ### code chunk number 128: morecolor2 ################################################### plotGrandLinear(gr.snp, aes(y = pvalue), color = c("green", "deepskyblue")) ################################################### ### code chunk number 129: morecolor3 ################################################### plotGrandLinear(gr.snp, aes(y = pvalue), color = c("green", "deepskyblue", "red")) ################################################### ### code chunk number 130: morecolor4 ################################################### plotGrandLinear(gr.snp, aes(y = pvalue), color = "red") ################################################### ### code chunk number 131: cutoff ################################################### plotGrandLinear(gr.snp, aes(y = pvalue), cutoff = 3, cutoff.color = "blue", cutoff.size = 4) ################################################### ### code chunk number 132: cutoff-low (eval = FALSE) ################################################### ## plotGrandLinear(gr.snp, aes(y = pvalue)) + geom_hline(yintercept = 3, color = "blue", size = 4) ################################################### ### code chunk number 133: longer ################################################### ## let's make a long name nms <- seqnames(seqinfo(gr.snp)) nms.new <- paste("chr00000", nms, sep = "") names(nms.new) <- nms gr.snp <- renameSeqlevels(gr.snp, nms.new) seqlengths(gr.snp) ################################################### ### code chunk number 134: rotate ################################################### plotGrandLinear(gr.snp, aes(y = pvalue)) + theme(axis.text.x=element_text(angle=-90, hjust=0)) ################################################### ### code chunk number 135: load ################################################### library(ggbio) library(BSgenome.Hsapiens.UCSC.hg19) data("genesymbol", package = "biovizBase") ################################################### ### code chunk number 136: load_bam ################################################### bamfile <- system.file("extdata", "SRR027894subRBM17.bam", package="biovizBase") library(Rsamtools) bf <- BamFile(bamfile) ################################################### ### code chunk number 137: BamFile ################################################### ggplot() + stat_mismatch(bf, which = genesymbol["RBM17"], bsgenome = Hsapiens,show.coverage = TRUE) + coord_cartesian(xlim = c(6134000, 6135000)) + theme_bw() ################################################### ### code chunk number 138: pag ################################################### library(biovizBase) pgr <- pileupAsGRanges(bamfile, region = genesymbol["RBM17"]) pgr.match <- pileupGRangesAsVariantTable(pgr, genome = Hsapiens) ################################################### ### code chunk number 139: pag_v ################################################### ggplot() + stat_mismatch(pgr.match) ################################################### ### code chunk number 140: pag_v2 ################################################### p1 <- ggplot() + stat_mismatch(pgr.match, show.coverage = FALSE, geom = "bar") + xlim(6132060,6132120) + ylim(0, 10) p2<- ggplot() + stat_mismatch(pgr.match, geom = "segment") + xlim(6132060,6132120) + ylim(0, 10) tracks(segment = p2, bar = p1) +scale_x_sequnit('Mb') ################################################### ### code chunk number 141: autoplot ################################################### autoplot(bf, which = genesymbol["RBM17"], bsgenome = Hsapiens,show.coverage = TRUE, stat = "mismatch", geom = "bar") + xlim(6132060,6132120) + ylim(0, 10) ################################################### ### code chunk number 142: make ################################################### library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(ggbio) data(genesymbol, package = "biovizBase") txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene model <- exonsBy(txdb, by = "tx") model17 <- subsetByOverlaps(model, genesymbol["RBM17"]) exons <- exons(txdb) exon17 <- subsetByOverlaps(exons, genesymbol["RBM17"]) ## reduce to make sure there is no overlap ## just for example exon.new <- reduce(exon17) ## suppose set.seed(1) values(exon.new)$sample1 <- rnorm(length(exon.new), 10, 3) values(exon.new)$sample2 <- rnorm(length(exon.new), 10, 10) values(exon.new)$significant <- c(TRUE, rep(FALSE,length(exon.new)-1)) head(exon.new) ################################################### ### code chunk number 143: link1 ################################################### plotRangesLinkedToData(exon.new, stat.col = 1:2) ## equivilent to ## plotRangesLinkedToData(exon.new, stat.col = c("sample1", "sample2")) ################################################### ### code chunk number 144: link3 ################################################### plotRangesLinkedToData(exon.new, stat.col = 1:2, size = 3, linetype = 4) ################################################### ### code chunk number 145: link4 ################################################### plotRangesLinkedToData(exon.new, stat.col = 1:2, size = 3, linetype = 4, sig = "significant", sig.col = c("gray70", "red")) ################################################### ### code chunk number 146: load ################################################### library(chipseq) library(GenomicFeatures) data(cstest) unique(seqnames(cstest)) ## subset chrs <- c("chr10", "chr11", "chr12") cstest <- keepSeqlevels(cstest, chrs) ## estimate fragment length fraglen <- estimate.mean.fraglen(cstest$ctcf) fraglen[!is.na(fraglen)] ## that's around 200 ## cstest.gr <- stack(cstest) ## head(cstest.gr) ## cstest.ext <- resize(cstest.gr, width = 200) ## extending them ctcf.ext <- resize(cstest$ctcf, width = 200) cov.ctcf <- coverage(ctcf.ext) gfp.ext <- resize(cstest$gfp, width = 200) cov.gfp <- coverage(gfp.ext) ## estimate peak cutoff peakCutoff(cov.ctcf, fdr = 0.0001) ## we can use 8 ################################################### ### code chunk number 147: find-close ################################################### c.p <- cstest$ctcf[seqnames(cstest$ctcf) == "chr10" & strand(cstest$ctcf) == "+",] c.n <- cstest$ctcf[seqnames(cstest$ctcf) == "chr10" & strand(cstest$ctcf) == "-",] cov.p <- coverage(c.p) cov.n <- coverage(c.n) v1 <- viewWhichMaxs(slice(cov.p, lower = 8))$chr10 v2 <- viewWhichMaxs(slice(cov.n, lower = 8))$chr10 res <- expand.grid(v1, v2) wh <- as.numeric(res[order(abs(res[,1] - res[, 2]))[1], ]) gr.wh <- GRanges("chr10", IRanges(wh[1], wh[2])) gr.wh <- resize(gr.wh, width(gr.wh) + 200, fix = "center") ################################################### ### code chunk number 148: reads-close ################################################### library(ggbio) ctcf.sub <- subsetByOverlaps(cstest$ctcf, gr.wh) p1 <- autoplot(ctcf.sub, aes(fill = strand), facets = strand ~ .) ctcf.ext.sub <- subsetByOverlaps(ctcf.ext, gr.wh) p2 <- autoplot(ctcf.ext.sub, aes(fill = strand), facets = strand ~ .) tracks("original" = p1, "extending" = p2, heights = c(3, 5)) ################################################### ### code chunk number 149: coverage-close ################################################### ctcf.sub <- subsetByOverlaps(cstest$ctcf, gr.wh) p1 <- autoplot(ctcf.sub, aes(fill = strand), facets = strand ~ ., stat = "coverage", geom = "area") ctcf.ext.sub <- subsetByOverlaps(ctcf.ext, gr.wh) p2 <- autoplot(ctcf.ext.sub, aes(fill = strand), facets = strand ~ ., stat = "coverage", geom = "area") tracks("original" = p1, "extending" = p2) ################################################### ### code chunk number 150: genome-bin-no-ylim ################################################### library(ggbio) p1 <- autoplot(cov.ctcf) p2 <- autoplot(cov.gfp) tracks(ctcf = p1, gfp = p2) ################################################### ### code chunk number 151: genome-bin-ylim ################################################### library(ggbio) p1 <- autoplot(cov.ctcf) p2 <- autoplot(cov.gfp) tracks(ctcf = p1, gfp = p2) + coord_cartesian(ylim = c(0, 2e6)) ################################################### ### code chunk number 152: genome-bin-maxs ################################################### p1 <- autoplot(cov.ctcf, type = "viewMaxs") p2 <- autoplot(cov.gfp,type = "viewMaxs") tracks(ctcf = p1, gfp = p2) + coord_cartesian(ylim = c(0, 2e6)) ################################################### ### code chunk number 153: genome-bin-100 ################################################### p1 <- autoplot(cov.ctcf, type = "viewMaxs", nbin = 100) p2 <- autoplot(cov.gfp,type = "viewMaxs", nbin = 100) tracks("ctcf" = p1, "gfp" = p2) + coord_cartesian(ylim = c(0, 1e6)) ################################################### ### code chunk number 154: genome-bin-heat ################################################### p1 <- autoplot(cov.ctcf, type = "viewMeans", nbin = 100, geom = "heatmap") p2 <- autoplot(cov.gfp,type = "viewMeans", nbin = 100, geom = "heatmap") tracks(ctcf = p1, gfp = p2) + scale_fill_continuous(limits = c(0, 8e05)) + scale_color_continuous(limits = c(0, 8e05)) ################################################### ### code chunk number 155: genome-slice ################################################### p1 <- autoplot(cov.ctcf, type = "viewMaxs", stat = "slice", lower = 8) p2 <- autoplot(cov.gfp,type = "viewMaxs", stat = "slice", lower = 8) tracks(ctcf = p1, gfp = p2) + coord_cartesian(ylim = c(0, 15000)) ################################################### ### code chunk number 156: genome-slice-drop ################################################### p1 <- autoplot(cov.ctcf, type = "viewMaxs", stat = "slice", lower = 8) p2 <- autoplot(cov.gfp,type = "viewMaxs", stat = "slice", lower = 8, drop = FALSE) tracks(ctcf = p1, gfp = p2) + coord_cartesian(ylim = c(0, 15000)) ################################################### ### code chunk number 157: genome-slice-rect ################################################### p1 <- autoplot(cov.ctcf, stat = "slice", lower = 5, geom = "rect") p2 <- autoplot(cov.gfp, stat = "slice", lower = 5, geom = "rect") tracks(ctcf = p1, gfp = p2) ################################################### ### code chunk number 158: genome-slice-heatmap ################################################### p1 <- autoplot(cov.ctcf, type = "viewMaxs", stat = "slice", lower = 8, geom = "heatmap") p2 <- autoplot(cov.gfp,type = "viewMaxs", stat = "slice", lower = 8, geom = "heatmap", drop = FALSE) tracks("ctcf" = p1, "gfp" = p2) + scale_fill_continuous(limits = c(1000, 6000)) + scale_color_continuous(limits = c(1000, 6000)) ################################################### ### code chunk number 159: region ################################################### peaks.ctcf <- slice(cov.ctcf, lower = 8) peaks.gfp <- slice(cov.gfp, lower = 8) ## this function is from vignette of chipseq peakSummary <- diffPeakSummary(peaks.gfp, peaks.ctcf) peakSummary <-within(peakSummary, { diffs <- asinh(sums2) - asinh(sums1) resids <- (diffs - median(diffs)) / mad(diffs) up <- resids > 2 down <- resids < -2 change <- ifelse(up, "up", ifelse(down, "down", "flat")) }) ps.down <- peakSummary[peakSummary$change == "down" & peakSummary$space == "chr11", ] pk.down <- ps.down[order(ps.down$diffs),] pk.down ## library(TxDb.Mmusculus.UCSC.mm9.knownGene) txdb <- TxDb.Mmusculus.UCSC.mm9.knownGene tx <- transcripts(txdb) gn <- transcriptsBy(txdb, by = "gene") fu <- fiveUTRsByTranscript(txdb) idx <- which(countOverlaps(as(pk.down, "GRanges"), flank(fu, width = 100)) == 1) wh.p <- as(pk.down[idx[2], ], "GRanges") wh.pw <- resize(wh.p, width = 30000, fix = "center") ################################################### ### code chunk number 160: ideogram ################################################### library(biovizBase) mm9 <- getIdeogram("mm9") cyto.def <- getOption("biovizBase")$cytobandColor cyto.new <- c(cyto.def, c(gpos33 = "grey80", gpos66 = "grey60")) p.ideo <- plotIdeogram(mm9, "chr10", zoom = c(start(wh.pw),end(wh.pw))) + scale_fill_manual(values = cyto.new) print(p.ideo) ################################################### ### code chunk number 161: txdb ################################################### p.gene <- autoplot(txdb, which = wh.pw) ################################################### ### code chunk number 162: ideo-features ################################################### ## coverage cstest.s <- stack(cstest) cstest.s <- resize(cstest.s, width = 200) cstest.sub <- subsetByOverlaps(cstest.s, wh.pw) p.cov <- autoplot(cstest.sub, stat = "coverage", facets = sample ~ ., geom = "area") ## ideogram tracks("ideogram" = p.ideo, "coverage" = p.cov, "gene" = p.gene, xlim = as(wh.pw, "GRanges"), xlim.change = c(FALSE, TRUE, TRUE), heights = c(1, 5, 5))