--- title: "trackViewer Vignette" author: "Jianhong Ou, Lihua Julie Zhu" date: "`r doc_date()`" package: "`r pkg_ver('trackViewer')`" abstract: > Visualize mapped reads along with annotation as track layers for NGS dataset such as ChIP-seq, RNA-seq, miRNA-seq, DNA-seq, SNPs and methylation data. vignette: > %\VignetteIndexEntry{trackViewer Vignette} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} output: BiocStyle::html_document --- ```{r, echo=FALSE, results="hide", warning=FALSE} suppressPackageStartupMessages({ library(trackViewer) library(rtracklayer) library(Gviz) library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(org.Hs.eg.db) library(VariantAnnotation) }) knitr::opts_chunk$set(warning=FALSE, message=FALSE) ``` # Introduction There are two packages available in Bioconductor for visualizing genomic data: **rtracklayer** and **Gviz**. **rtracklayer** provides an interface to genome browsers and associated annotation tracks. **Gviz** plots data and annotation information along genomic coordinates. **TrackViewer** is a light-weighted visualization tool for generating neat figures for publication. It utilizes **Gviz**, is easy to use, and has a low memory and cpu consumption. ```{r plotComp,echo=TRUE,fig.keep='none'} library(Gviz) library(rtracklayer) library(trackViewer) extdata <- system.file("extdata", package="trackViewer", mustWork=TRUE) gr <- GRanges("chr11", IRanges(122929275, 122930122), strand="-") fox2 <- importScore(file.path(extdata, "fox2.bed"), format="BED", ranges=gr) fox2$dat <- coverageGR(fox2$dat) viewTracks(trackList(fox2), gr=gr, autoOptimizeStyle=TRUE, newpage=FALSE) dt <- DataTrack(range=fox2$dat[strand(fox2$dat)=="-"] , genome="hg19", type="hist", name="fox2", window=-1, chromosome="chr11", fill.histogram="black", col.histogram="NA", background.title="white", col.frame="white", col.axis="black", col="black", col.title="black") plotTracks(dt, from=122929275, to=122930122, strand="-") ``` ```{r Gviz,echo=FALSE,fig.cap='Plot data with **Gviz** and **trackViewer**. Note that **trackViewer** can generate similar figure as **Gviz** with several lines of simple codes.',fig.width=6,fig.height=1.5} viewerStyle <- trackViewerStyle() setTrackViewerStyleParam(viewerStyle, "margin", c(.01, .13, .02, .02)) empty <- DataTrack(showAxis=FALSE, showTitle=FALSE, background.title="white") plotTracks(list(empty, dt), from=122929275, to=122930122, strand="-") pushViewport(viewport(0, .5, 1, .5, just=c(0, 0))) viewTracks(trackList(fox2), viewerStyle=viewerStyle, gr=gr, autoOptimizeStyle=TRUE, newpage=FALSE) popViewport() grid.text(label="Gviz track", x=.3, y=.4) grid.text(label="trackViewer track", x=.3, y=.9) ``` **TrackViewer** not only has the functionalities to plot the figures generated by **Gviz**, as shown in Figure above, but also provides additional plotting styles as shown in Figure below. The mimimalist design requires minimum input from users while retaining the flexibility to change output style easily. ```{r lostcode,echo=TRUE,fig.keep='none'} gr <- GRanges("chr1", IRanges(c(1, 6, 10), c(3, 6, 12)), score=c(3, 4, 1)) dt <- DataTrack(range=gr, data="score", type="hist") plotTracks(dt, from=2, to=11) tr <- new("track", dat=gr, type="data", format="BED") viewTracks(trackList(tr), chromosome="chr1", start=2, end=11) ``` ```{r GvizLost,echo=FALSE,fig.cap='Plot data with **Gviz** and **trackViewer**. Note that **trackViewer** is not only including more details but also showing all the data involved in the given range.',fig.width=6,fig.height=1.5} plotTracks(list(empty, dt), from=2, to=11) pushViewport(viewport(0, .5, 1, .5, just=c(0, 0))) viewTracks(trackList(tr), viewerStyle=viewerStyle, chromosome="chr1", start=2, end=11, autoOptimizeStyle=TRUE, newpage=FALSE) popViewport() grid.text(label="Gviz track", x=.3, y=.4) grid.text(label="trackViewer track", x=.3, y=.9) ``` It requires huge memory space to handle big wig files. To solve this problem, **trackViewer** rewrote the import function to import whole file first and parse it later when plot. **trackViewer** provides higher import speed (21 min vs. over 180 min) and acceptable memory cost (5.32G vs. over 10G) for a half giga wig file (GSM917672) comparing to **Gviz**. # Browse tracks ## Steps of using \Biocpkg{trackViewer} ### step1 import data Function **importScore** is used to import BED, WIG, bedGraph or BigWig files. Function **importBam** is employed to import bam file. Here is the example. ```{r importData} library(trackViewer) extdata <- system.file("extdata", package="trackViewer", mustWork=TRUE) repA <- importScore(file.path(extdata, "cpsf160.repA_-.wig"), file.path(extdata, "cpsf160.repA_+.wig"), format="WIG") ## because the wig file does not contain strand info, ## we need to set it manually strand(repA$dat) <- "-" strand(repA$dat2) <- "+" ``` Function **coverageGR** could be used to calculate coverage after importing if needed. ```{r coverage} fox2 <- importScore(file.path(extdata, "fox2.bed"), format="BED", ranges=GRanges("chr11", IRanges(122929000, 122931000))) dat <- coverageGR(fox2$dat) ## we can split the data by strand into two different track channels ## here we set the dat2 slot to save the negative strand info, ## reverse order as previous. fox2$dat <- dat[strand(dat)=="+"] fox2$dat2 <- dat[strand(dat)=="-"] ``` ### step2 build gene model The gene model can be built for a given genomic range using **geneModelFromTxdb** function which uses **TranscriptDb** object as input. ```{r geneModel} library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(org.Hs.eg.db) gr <- GRanges("chr11", IRanges(122929275, 122930122), strand="-") trs <- geneModelFromTxdb(TxDb.Hsapiens.UCSC.hg19.knownGene, org.Hs.eg.db, gr=gr) ``` ### step3 view the tracks Use **viewTracks** function to plot data and annotation information along genomic coordinates. **addGuideLine** or **addArrowMark** can be used to highlight the peaks. ```{r viewTracks,fig.cap='plot data and annotation information along genomic coordinates',fig.width=6,fig.height=1.5} viewerStyle <- trackViewerStyle() setTrackViewerStyleParam(viewerStyle, "margin", c(.1, .05, .02, .02)) vp <- viewTracks(trackList(repA, fox2, trs), gr=gr, viewerStyle=viewerStyle, autoOptimizeStyle=TRUE) addGuideLine(c(122929767, 122929969), vp=vp) addArrowMark(list(x=122929650, y=2), # 2 means track 2 from bottom. label="label", col="blue", vp=vp) ``` ## Adjust the styles ### adjust x axis or x-scale In most cases, researchers are interested in the relative position of peaks in the gene. Sometimes, margin needs to be adjusted to be able to show the entire gene model. Figure below shows how to add an x-scale and remove x-axis using **addGuideLine** Function . ```{r optSty} optSty <- optimizeStyle(trackList(repA, fox2, trs)) trackList <- optSty$tracks viewerStyle <- optSty$style ``` ```{r viewTracksXaxis,fig.cap='plot data with x-scale',fig.width=6,fig.height=1.5} setTrackViewerStyleParam(viewerStyle, "xaxis", FALSE) setTrackViewerStyleParam(viewerStyle, "margin", c(.01, .05, .01, .01)) setTrackXscaleParam(trackList[[1]], "draw", TRUE) setTrackXscaleParam(trackList[[1]], "gp", list(cex=.5)) viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ``` ### adjust y axis y-axis can be put to right side of the track by setting main slot to FALSE in y-axis slot of each track. And ylim can be set by **setTrackStyleParam**. ```{r viewTracksYaxis,fig.cap='plot data with y-axis in right side',fig.width=6,fig.height=1.5} setTrackViewerStyleParam(viewerStyle, "margin", c(.01, .05, .01, .05)) for(i in 1:2){ setTrackYaxisParam(trackList[[i]], "main", FALSE) } ## adjust y scale setTrackStyleParam(trackList[[1]], "ylim", c(0, 25)) setTrackStyleParam(trackList[[2]], "ylim", c(-25, 0)) viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ``` ### adjust y label Y label style can be changed by setting the ylabgp slot in style of each track. ```{r viewTracksYlab,fig.cap='plot data with adjusted color and size of y label',fig.width=6,fig.height=1.5} setTrackStyleParam(trackList[[1]], "ylabgp", list(cex=.8, col="green")) ## set cex to avoid automatic adjust setTrackStyleParam(trackList[[2]], "ylabgp", list(cex=.8, col="blue")) setTrackStyleParam(trackList[[2]], "marginBottom", .2) viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ``` Y label can be also put to top or bottom of each track. ```{r viewTracksYlabTopBottom,fig.cap='plot data with adjusted y label position',fig.width=6,fig.height=1.5} setTrackStyleParam(trackList[[1]], "ylabpos", "bottomleft") setTrackStyleParam(trackList[[2]], "ylabpos", "topright") setTrackStyleParam(trackList[[2]], "marginTop", .2) viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ``` For each transcript, the transcript name can be put to upstream or downstream of the transcript. ```{r viewTracksYlabUpsDown,fig.cap='plot data with adjusted transcripts name position',fig.width=6,fig.height=1.5} trackListN <- trackList setTrackStyleParam(trackListN[[3]], "ylabpos", "upstream") setTrackStyleParam(trackListN[[4]], "ylabpos", "downstream") ## set cex to avoid automatic adjust setTrackStyleParam(trackListN[[3]], "ylabgp", list(cex=.6)) setTrackStyleParam(trackListN[[4]], "ylabgp", list(cex=.6)) gr1 <- range(unlist(GRangesList(sapply(trs, function(.ele) .ele$dat)))) start(gr1) <- start(gr1) - 2000 end(gr1) <- end(gr1) + 2000 viewTracks(trackListN, gr=gr1, viewerStyle=viewerStyle) ``` ### adjust track color The track color can be changed by setting the color slot in style of each track. The first color is for dat slot of **track** and seconde color is for dat2 slot. ```{r viewTracksCol,fig.cap='plot data with adjusted track color',fig.width=6,fig.height=1.5} setTrackStyleParam(trackList[[1]], "color", c("green", "black")) setTrackStyleParam(trackList[[2]], "color", c("black", "blue")) for(i in 3:length(trackList)) setTrackStyleParam(trackList[[i]], "color", "black") viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ``` ### adjust track height The track height can be changed by setting the height slot in style of each track. However, the total height for all the tracks should be 1. ```{r viewTracksHeight,fig.cap='plot data with adjusted track height',fig.width=6,fig.height=1.5} trackListH <- trackList setTrackStyleParam(trackListH[[1]], "height", .1) setTrackStyleParam(trackListH[[2]], "height", .44) for(i in 3:length(trackListH)){ setTrackStyleParam(trackListH[[i]], "height", (1-(0.1+0.44))/(length(trackListH)-2)) } viewTracks(trackListH, gr=gr, viewerStyle=viewerStyle) ``` ### change track names The track names such as gene model names can also be edited easily by changing the names of **trackList**. ```{r viewTracksNames,fig.cap='change the track names',fig.width=6,fig.height=1.5} names(trackList) <- c("cpsf160", "fox2", rep("HSPA8", 5)) viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ``` ### show paired data in the same track **trackViewer** can be used to show to-be-compared data in the same track side by side. ```{r viewTracksPaired,fig.cap='show two data in the same track',fig.width=6,fig.height=1.2} cpsf160 <- importScore(file.path(extdata, "cpsf160.repA_-.wig"), file.path(extdata, "cpsf160.repB_-.wig"), format="WIG") strand(cpsf160$dat) <- strand(cpsf160$dat2) <- "-" setTrackStyleParam(cpsf160, "color", c("black", "red")) viewTracks(trackList(trs, cpsf160), gr=gr, viewerStyle=viewerStyle) ``` ### flip the x-axis The x-axis can be horizotally flipped for the genes in negative strand. ```{r viewTracksFlipped,fig.cap='show data in the flipped track',fig.width=6,fig.height=2} viewerStyleF <- viewerStyle setTrackViewerStyleParam(viewerStyleF, "flip", TRUE) setTrackViewerStyleParam(viewerStyleF, "xaxis", TRUE) setTrackViewerStyleParam(viewerStyleF, "margin", c(.1, .05, .01, .01)) vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyleF) addGuideLine(c(122929767, 122929969), vp=vp) addArrowMark(list(x=122929650, y=2), label="label", col="blue", vp=vp) ``` ### optimize with theme We support two themes now: bw and col. ```{r theme_bw,fig.cap='theme_bw',fig.width=6,fig.height=2} optSty <- optimizeStyle(trackList(repA, fox2, trs), theme="bw") trackList <- optSty$tracks viewerStyle <- optSty$style vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ``` ```{r theme_col,fig.cap='theme_col',fig.width=6,fig.height=2} optSty <- optimizeStyle(trackList(repA, fox2, trs), theme="col") trackList <- optSty$tracks viewerStyle <- optSty$style vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ``` ### plot with breaks We could plot the tracks with breaks by set multiple genomic ranges. ```{r axis.break,fig.cap='axis.break',fig.width=6,fig.height=2} gr.breaks <- GRanges("chr11", IRanges(c(122929275, 122929575, 122929775), c(122929555, 122929725, 122930122)), strand="-", percentage=c(.4, .2, .4)) vp <- viewTracks(trackList, gr=gr.breaks, viewerStyle=viewerStyle) ``` ## browse tracks by htmlwidgets To change the parameter may be complicated for users. Users may want to have a try function `browseTracks`. ```{r browseTrack,fig.cap='browseTrack',fig.width=8,fig.height=6} browseTracks(trackList, gr=gr) ``` # Operators If there are two tracks and we want to draw the two track by adding or substract one from another, we can try operators. ```{r viewTracksOperator1,fig.cap='show data with operator "+"',fig.width=6,fig.height=2} newtrack <- repA ## must keep same format for dat and dat2 newtrack <- parseWIG(newtrack, "chr11", 122929275, 122930122) newtrack$dat2 <- newtrack$dat newtrack$dat <- fox2$dat2 setTrackStyleParam(newtrack, "color", c("blue", "red")) viewTracks(trackList(newtrack, trs), gr=gr, viewerStyle=viewerStyle, operator="+") ``` ```{r viewTracksOperator2,fig.cap='show data with operator "-"',fig.width=6,fig.height=2} viewTracks(trackList(newtrack, trs), gr=gr, viewerStyle=viewerStyle, operator="-") ``` Or try **GRoperator** before view tracks. ```{r viewTracksOperator3,fig.cap='show data with operator "-"',fig.width=6,fig.height=2} newtrack$dat <- GRoperator(newtrack$dat, newtrack$dat2, col="score", operator="-") newtrack$dat2 <- GRanges() viewTracks(trackList(newtrack, trs), gr=gr, viewerStyle=viewerStyle) ``` # Lolliplot Lolliplot is a mutation distribution graphics tool. ```{r lolliplot1, fig.width=6, fig.height=3} SNP <- c(10, 12, 1400, 1402) sample.gr <- GRanges("chr1", IRanges(SNP, width=1, names=paste0("snp", SNP))) features <- GRanges("chr1", IRanges(c(1, 501, 1001), width=c(120, 400, 405), names=paste0("block", 1:3))) lolliplot(sample.gr, features) ## more SNPs SNP <- c(10, 100, 105, 108, 400, 410, 420, 600, 700, 805, 840, 1400, 1402) sample.gr <- GRanges("chr1", IRanges(SNP, width=1, names=paste0("snp", SNP))) lolliplot(sample.gr, features) ## define the range lolliplot(sample.gr, features, ranges = GRanges("chr1", IRanges(104, 109))) ``` ## Change lolliplot colors ### Change the color of features. ```{r lolliplot2, fig.width=6, fig.height=3} features$fill <- c("#FF8833", "#51C6E6", "#DFA32D") lolliplot(sample.gr, features) ``` ### Change the color of lollipop. ```{r lolliplot3, fig.width=6, fig.height=3} sample.gr$color <- sample.int(6, length(SNP), replace=TRUE) sample.gr$border <- sample(c("gray80", "gray30"), length(SNP), replace=TRUE) lolliplot(sample.gr, features) ``` ## Add index labels in the node ```{r lolliplot.index, fig.width=6, fig.height=3} sample.gr$label <- as.character(1:length(sample.gr)) sample.gr$label.col <- "white" lolliplot(sample.gr, features) ``` ## Change the height of features ```{r lolliplot4, fig.width=6, fig.height=3} features$height <- c(0.02, 0.05, 0.08) lolliplot(sample.gr, features) ## keep the height by giving the unit features$height <- list(unit(1/16, "inches"), unit(3, "mm"), unit(12, "points")) lolliplot(sample.gr, features) ``` ## Multiple transcripts in the features The metadata 'featureLayerID' are used for drawing features in different layers. ```{r lolliplot.mul.features, fig.width=6, fig.height=3} features.mul <- rep(features, 2) features.mul$height[4:6] <- list(unit(1/8, "inches"), unit(0.5, "lines"), unit(.2, "char")) features.mul$fill <- c("#FF8833", "#F9712A", "#DFA32D", "#51C6E6", "#009DDA", "#4B9CDF") end(features.mul)[5] <- end(features.mul[5])+50 features.mul$featureLayerID <- paste("tx", rep(1:2, each=length(features)), sep="_") names(features.mul) <- paste(features.mul$featureLayerID, rep(1:length(features), 2), sep="_") lolliplot(sample.gr, features.mul) ## one name per transcripts names(features.mul) <- features.mul$featureLayerID lolliplot(sample.gr, features.mul) ``` ## Change the height of lolliplot ```{r lolliplot4.1, fig.width=6, fig.height=3.5} #Note: the score value is integer less than 10 sample.gr$score <- sample.int(5, length(sample.gr), replace = TRUE) lolliplot(sample.gr, features) ##remove yaxis lolliplot(sample.gr, features, yaxis=FALSE) ``` ```{r lolliplot4.2, fig.width=6, fig.height=4.5} #Try score value greater than 10 sample.gr$score <- sample.int(20, length(sample.gr), replace=TRUE) lolliplot(sample.gr, features) #Try float numeric score sample.gr$score <- runif(length(sample.gr))*10 lolliplot(sample.gr, features) # score should not be smaller than 1 ``` ## Customize xaxis label positions ```{r lolliplot.xaxis, fig.width=6, fig.height=4.5} xaxis <- c(1, 200, 400, 701, 1000, 1200, 1402) ## define the position lolliplot(sample.gr, features, xaxis=xaxis) names(xaxis) <- xaxis # define the labels names(xaxis)[4] <- "center" lolliplot(sample.gr, features, xaxis=xaxis) ``` ## Customize yaxis label positions ```{r lolliplot.yaxis, fig.width=6, fig.height=4.5} yaxis <- c(0, 5) ## define the position lolliplot(sample.gr, features, yaxis=yaxis) yaxis <- c(0, 5, 10, 15) names(yaxis) <- yaxis # define the labels names(yaxis)[3] <- "yaxis" lolliplot(sample.gr, features, yaxis=yaxis) ``` ## Jitter the label only ```{r lolliplot.jitter, fig.width=6, fig.height=4.5} sample.gr$dashline.col <- sample.gr$color lolliplot(sample.gr, features, jitter="label") ``` ## Add legend ```{r lolliplot.legend, fig.width=6, fig.height=5} legend <- 1:6 ## legend fill color names(legend) <- paste0("legend", letters[1:6]) ## legend labels lolliplot(sample.gr, features, legend=legend) ## use list to define more attributes. see ?grid::gpar to get more details. legend <- list(labels=paste0("legend", LETTERS[1:6]), col=palette()[6:1], fill=palette()[legend]) lolliplot(sample.gr, features, legend=legend) ## if you have multiple tracks, try to set the legend by list. ## see more in section [Plot multiple samples](#plot-multiple-samples) legend <- list(legend) lolliplot(sample.gr, features, legend=legend) ``` ## Control the labels Use can control the paramters of labels by name the metadata start as label.parameter. such as label.parameter.rot, label.parameter.gp. The parameter is used for \link[grid]{grid.text}. ```{r lolliplot.labels.control, fig.width=6, fig.height=5} sample.gr.rot <- sample.gr sample.gr.rot$label.parameter.rot <- 45 lolliplot(sample.gr.rot, features, legend=legend) sample.gr.rot$label.parameter.rot <- 60 sample.gr.rot$label.parameter.gp <- gpar(col="brown") lolliplot(sample.gr.rot, features, legend=legend) ``` If you want to change the text in ylab, please try to set the labels in ylab. **lolliplot** does not support the parameters for title and xlab. If you want to add title and xlab, please try to add them by \link[grid]{grid.text}. ```{r lolliplot.xlab.ylab.title, fig.width=6, fig.height=5.2} lolliplot(sample.gr.rot, features, legend=legend, ylab="y label here") grid.text("x label here", x=.5, y=.01, just="bottom") grid.text("title here", x=.5, y=.98, just="top", gp=gpar(cex=1.5, fontface="bold")) ``` ## Change the type of lolliplot ### Google pin ```{r lolliplot5, fig.width=6, fig.height=4.5} lolliplot(sample.gr, features, type="pin") sample.gr$color <- lapply(sample.gr$color, function(.ele) c(.ele, sample.int(6, 1))) sample.gr$border <- sample.int(6, length(SNP), replace=TRUE) lolliplot(sample.gr, features, type="pin") ``` ### Pie plot ```{r lolliplot6, fig.width=6, fig.height=3} sample.gr$score <- NULL ## must be removed, because pie will consider all the numeric columns except column "color", "fill", "lwd", "id" and "id.col". sample.gr$label <- NULL sample.gr$label.col <- NULL x <- sample.int(100, length(SNP)) sample.gr$value1 <- x sample.gr$value2 <- 100 - x ## the length of color should be no less than the values number sample.gr$color <- rep(list(c("#87CEFA", '#98CE31')), length(SNP)) sample.gr$border <- "gray30" lolliplot(sample.gr, features, type="pie") ``` ## Plot multiple samples ### multiple layers ```{r lolliplot7, fig.width=6, fig.height=5.5} SNP2 <- sample(4000:8000, 30) x2 <- sample.int(100, length(SNP2), replace=TRUE) sample2.gr <- GRanges("chr3", IRanges(SNP2, width=1, names=paste0("snp", SNP2)), value1=x2, value2=100-x2) sample2.gr$color <- rep(list(c('#DB7575', '#FFD700')), length(SNP2)) sample2.gr$border <- "gray30" features2 <- GRanges("chr3", IRanges(c(5001, 5801, 7001), width=c(500, 500, 405), names=paste0("block", 4:6)), fill=c("orange", "gray30", "lightblue"), height=unit(c(0.5, 0.3, 0.8), "cm")) legends <- list(list(labels=c("WT", "MUT"), fill=c("#87CEFA", '#98CE31')), list(labels=c("WT", "MUT"), fill=c('#DB7575', '#FFD700'))) lolliplot(list(A=sample.gr, B=sample2.gr), list(x=features, y=features2), type="pie", legend=legends) ``` Different layouts are also be possible. ```{r lolliplot.multiple.type, fig.width=6, fig.height=7.5} sample2.gr$score <- sample2.gr$value1 ## circle layout need score column lolliplot(list(A=sample.gr, B=sample2.gr), list(x=features, y=features2), type=c("pie", "circle"), legend=legends) ``` ### pie.stack layout ```{r lolliplot.pie.stack, fig.width=6, fig.height=5} rand.id <- sample.int(length(sample.gr), 3*length(sample.gr), replace=TRUE) rand.id <- sort(rand.id) sample.gr.mul.patient <- sample.gr[rand.id] ## pie.stack require metadata "stack.factor", and the metadata can not be ## stack.factor.order or stack.factor.first len.max <- max(table(rand.id)) stack.factors <- paste0("patient", formatC(1:len.max, width=nchar(as.character(len.max)), flag="0")) sample.gr.mul.patient$stack.factor <- unlist(lapply(table(rand.id), sample, x=stack.factors)) sample.gr.mul.patient$value1 <- sample.int(100, length(sample.gr.mul.patient), replace=TRUE) sample.gr.mul.patient$value2 <- 100 - sample.gr.mul.patient$value1 patient.color.set <- as.list(as.data.frame(rbind(rainbow(length(stack.factors)), "#FFFFFFFF"), stringsAsFactors=FALSE)) names(patient.color.set) <- stack.factors sample.gr.mul.patient$color <- patient.color.set[sample.gr.mul.patient$stack.factor] legend <- list(labels=stack.factors, col="gray80", fill=sapply(patient.color.set, `[`, 1)) lolliplot(sample.gr.mul.patient, features, type="pie.stack", legend=legend, dashline.col="gray") ``` ### caterpillar layout Metadata SNPsideID is used to trigger caterpillar layout. SNPsideID must be 'top' or 'bottom'. ```{r lolliplot.caterpillar, fig.width=6, fig.height=4} sample.gr$SNPsideID <- sample(c("top", "bottom"), length(sample.gr), replace=TRUE) lolliplot(sample.gr, features, type="pie", legend=legends[[1]]) ``` ```{r lolliplot.caterpillar2, fig.width=6, fig.height=12} ## two layers sample2.gr$SNPsideID <- "top" idx <- sample.int(length(sample2.gr), 15) sample2.gr$SNPsideID[idx] <- "bottom" sample2.gr$color[idx] <- '#FFD700' lolliplot(list(A=sample.gr, B=sample2.gr), list(x=features.mul, y=features2), type=c("pie", "circle"), legend=legends) ``` ## Variant Call Format (VCF) data ```{r VCF, fig.width=6, fig.height=5} library(VariantAnnotation) library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(org.Hs.eg.db) fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation") gr <- GRanges("22", IRanges(50968014, 50970514, names="TYMP")) tab <- TabixFile(fl) vcf <- readVcf(fl, "hg19", param=gr) mutation.frequency <- rowRanges(vcf) mcols(mutation.frequency) <- cbind(mcols(mutation.frequency), VariantAnnotation::info(vcf)) mutation.frequency$border <- "gray30" mutation.frequency$color <- ifelse(grepl("^rs", names(mutation.frequency)), "lightcyan", "lavender") ## plot Global Allele Frequency based on AC/AN mutation.frequency$score <- mutation.frequency$AF*100 seqlevelsStyle(gr) <- seqlevelsStyle(mutation.frequency) <- "UCSC" trs <- geneModelFromTxdb(TxDb.Hsapiens.UCSC.hg19.knownGene, org.Hs.eg.db, gr=gr) features <- c(range(trs[[1]]$dat), range(trs[[5]]$dat)) names(features) <- c(trs[[1]]$name, trs[[5]]$name) features$fill <- c("lightblue", "mistyrose") features$height <- c(.02, .04) lolliplot(mutation.frequency, features, ranges=gr) ``` ## Methylation data ```{r methylation, eval=FALSE, echo=TRUE} library(rtracklayer) session <- browserSession() query <- ucscTableQuery(session, track="HAIB Methyl RRBS", range=GRangesForUCSCGenome("hg19", seqnames(gr), ranges(gr))) tableName(query) <- tableNames(query)[1] methy <- track(query) methy <- GRanges(methy) ``` ```{r methylation.hide, echo=FALSE} methy <- import(system.file("extdata", "methy.bed", package="trackViewer"), "BED") ``` ```{r fig.width=6,fig.height=4} lolliplot(methy, features, ranges=gr, type="pin") ``` ## Change node size In above sample, some of the nodes overlap each other. To change the node size, cex prameter could be used. ```{r fig.width=6,fig.height=2.5} methy$lwd <- .5 lolliplot(methy, features, ranges=gr, type="pin", cex=.5) lolliplot(methy, features, ranges=gr, type="circle", cex=.5) methy$score2 <- max(methy$score) - methy$score lolliplot(methy, features, ranges=gr, type="pie", cex=.5) ## we can change it one by one methy$cex <- runif(length(methy)) lolliplot(methy, features, ranges=gr, type="pin") lolliplot(methy, features, ranges=gr, type="circle") ``` ## Change the xscale In above samples, some of the nodes are moved too far from the original coordinates. To rescale the xaxis could be used. ```{r fig.width=6,fig.height=2.5} methy$cex <- 1 lolliplot(methy, features, ranges=gr, rescale = TRUE) rescale <- data.frame( from.start = c(50968014, 50968515, 50968838), from.end = c(50968514, 50968837, 50970514), to.start = c(50968014, 50968838, 50969501), to.end = c(50968837, 50969500, 50970514) ) xaxis <- c(50968014, 50968514, 50968710, 50968838, 50970514) lolliplot(methy, features, ranges=gr, type="pin", rescale = rescale, xaxis = xaxis) ``` # Dandelion Plot Sometimes, there are too many SNPs invoved in one gene. If you want to plot such as more than handred SNPs, Dandelion plot may be the good chioce. Please note, the height of the dandelion means the desity of the SNPs. ```{r fig.width=4.5,fig.height=3} dandelion.plot(methy, features, ranges=gr, type="pin") ``` ## change the type of Dandelion plot There are one more type for dandelion plot, type "fan". The area of the fan indicate the percentage of methylation or rate of mutation. ```{r fig.width=4.5,fig.height=3} methy$color <- 3 methy$border <- "gray" ## score info is required, the score must be a number in [0, 1] m <- max(methy$score) methy$score <- methy$score/m dandelion.plot(methy, features, ranges=gr, type="fan") ``` ```{r fig.width=4.5,fig.height=4} methy$color <- rep(list(c(3, 5)), length(methy)) methy$score2 <- methy$score2/m legends <- list(list(labels=c("s1", "s2"), fill=c(3, 5))) dandelion.plot(methy, features, ranges=gr, type="pie", legend=legends) ``` ## change the number of dandelions ```{r fig.width=4.5,fig.height=2.5} ## want less dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=1/10) ``` ```{r fig.width=4.5,fig.height=4.5} ## want more dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=1/100) ``` # Plot lollipop plot with tracks ```{r fig.width=6,fig.height=2} gene <- geneTrack(get("HSPA8", org.Hs.egSYMBOL2EG), TxDb.Hsapiens.UCSC.hg19.knownGene)[[1]] SNPs <- GRanges("chr11", IRanges(sample(122929275:122930122, size = 20), width = 1), strand="-") SNPs$score <- sample.int(5, length(SNPs), replace = TRUE) SNPs$color <- sample.int(6, length(SNPs), replace=TRUE) SNPs$border <- "gray80" SNPs$feature.height = .1 SNPs$cex <- .5 gene$dat2 <- SNPs optSty <- optimizeStyle(trackList(repA, fox2, gene), theme="col") trackList <- optSty$tracks viewerStyle <- optSty$style gr <- GRanges("chr11", IRanges(122929275, 122930122)) setTrackStyleParam(trackList[[3]], "ylabgp", list(cex=.8)) vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) ## lollipopData track SNPs2 <- GRanges("chr11", IRanges(sample(122929275:122930122, size = 30), width = 1), strand="-") SNPs2 <- c(SNPs2, promoters(gene$dat, upstream = 0, downstream = 1)) SNPs2$score <- sample.int(3, length(SNPs2), replace = TRUE) SNPs2$color <- sample.int(6, length(SNPs2), replace=TRUE) SNPs2$border <- "gray30" SNPs2$feature.height = .1 SNPs2$cex <- .5 SNPs$cex <- .5 lollipopData <- new("track", dat=SNPs, dat2=SNPs2, type="lollipopData") gene <- geneTrack(get("HSPA8", org.Hs.egSYMBOL2EG), TxDb.Hsapiens.UCSC.hg19.knownGene)[[1]] optSty <- optimizeStyle(trackList(repA, lollipopData, gene, heightDist = c(3, 3, 1)), theme="col") trackList <- optSty$tracks viewerStyle <- optSty$style gr <- GRanges("chr11", IRanges(122929275, 122930122)) setTrackStyleParam(trackList[[2]], "ylabgp", list(cex=.8)) vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle) addGuideLine(122929538, vp=vp) ``` # Ideogram Plot Plot ideograms with data. ```{r eval=FALSE} ideo <- loadIdeogram("hg38") ``` ```{r echo=FALSE, results="hide"} path <- system.file("extdata", "ideo.hg38.rds", package = "trackViewer") ideo <- readRDS(path) ``` ```{r} dataList <- ideo dataList$score <- as.numeric(dataList$gieStain) dataList <- dataList[dataList$gieStain!="gneg"] dataList <- GRangesList(dataList) ideogramPlot(ideo, dataList, layout=list("chr1", c("chr3", "chr22"), c("chr4", "chr21"))) ``` # Plot genomic interactions data Plot GInteractions. Different from most of the available tools, plotGInteractions try to plot the data with 2D structure. The nodes indicate the region with interactions and the edges indicates the interactions. The size of the nodes are relative to the width of the region. The features could be enhancer, promoter or gene. The enhancer and promoter are shown as points with symbol 11 and 13. ```{r} library(TxDb.Hsapiens.UCSC.hg19.knownGene) library(InteractionSet) gi <- readRDS(system.file("extdata", "gi.rds", package="trackViewer")) range <- GRanges("chr2", IRanges(234500000, 235000000)) feature.gr <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene) feature.gr <- subsetByOverlaps(feature.gr, regions(gi)) feature.gr$col <- sample(1:7, length(feature.gr), replace=TRUE) feature.gr$type <- sample(c("promoter", "enhancer", "gene"), length(feature.gr), replace=TRUE, prob=c(0.1, 0.2, 0.7)) plotGInteractions(gi, range, feature.gr) ``` # Session Info ```{r sessionInfo, results='asis'} sessionInfo() ```