inst/doc/trackViewer.R

In trackViewer: A R/Bioconductor package with web interface for drawing elegant interactive tracks or lollipop plot to facilitate integrated analysis of multi-omics data

## ---- 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)
  library(httr)
})
knitr::opts_chunk$set(warning=FALSE, message=FALSE)

## ----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="-")

## ----Gviz,echo=FALSE,fig.cap='Plot data with **Gviz** and **trackViewer**. Please note that **trackViewer** can generate similar figure as **Gviz** with several lines of simple codes.',fig.width=8,fig.height=3----
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)

## ----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)

## ----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=8,fig.height=3----
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)

## ----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 any strand info, 
## we need to set it manually.
strand(repA$dat) <- "-"
strand(repA$dat2) <- "+"

## ----coverage-----------------------------------------------------------------
fox2 <- importScore(file.path(extdata, "fox2.bed"), format="BED",
                    ranges=GRanges("chr11", IRanges(122830799, 123116707)))
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. 
 
fox2$dat <- dat[strand(dat)=="+"]
fox2$dat2 <- dat[strand(dat)=="-"]

## ----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)

## ----geneTrack----------------------------------------------------------------
entrezIDforFMR1 <- get("FMR1", org.Hs.egSYMBOL2EG)
theTrack <- geneTrack(entrezIDforFMR1,TxDb.Hsapiens.UCSC.hg19.knownGene)[[1]]

## ----viewTracks,fig.cap='plot data and annotation information along genomic coordinates',fig.width=8,fig.height=3----
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 the bottom.
             label="label",
             col="blue",
             vp=vp)

## ----optSty-------------------------------------------------------------------
optSty <- optimizeStyle(trackList(repA, fox2, trs))
trackList <- optSty$tracks
viewerStyle <- optSty$style

## ----viewTracksXaxis,fig.cap='plot data with x-scale',fig.width=8,fig.height=3----
setTrackViewerStyleParam(viewerStyle, "xaxis", FALSE)
setTrackViewerStyleParam(viewerStyle, "margin", c(.01, .05, .01, .01))
setTrackXscaleParam(trackList[[1]], "draw", TRUE)
setTrackXscaleParam(trackList[[1]], "gp", list(cex=0.8))
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)
setTrackXscaleParam(trackList[[1]], attr="position", 
                    value=list(x=122929700, y=3, label=200))
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)

## ----viewTracksYaxis,fig.cap='plot data with y-axis in right side',fig.width=8,fig.height=3----
setTrackViewerStyleParam(viewerStyle, "margin", c(.01, .05, .01, .05))
for(i in 1:2){
    setTrackYaxisParam(trackList[[i]], "main", FALSE)
}
## Adjust the limit of y-axis
setTrackStyleParam(trackList[[1]], "ylim", c(0, 25))
setTrackStyleParam(trackList[[2]], "ylim", c(-25, 0))
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)

## ----viewTracksYlab,fig.cap='plot data with adjusted color and size of y-axis label',fig.width=8,fig.height=3----
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)

## ----viewTracksYlabTopBottom,fig.cap='plot data with adjusted y-axis label position',fig.width=8,fig.height=3----
setTrackStyleParam(trackList[[1]], "ylabpos", "bottomleft")
setTrackStyleParam(trackList[[2]], "ylabpos", "topright")
setTrackStyleParam(trackList[[2]], "marginTop", .2)
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)

## ----viewTracksYlabUpsDown,fig.cap='plot data with adjusted transcripts name position',fig.width=8,fig.height=3----
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)

## ----viewTracksCol,fig.cap='plot data with adjusted track color',fig.width=8,fig.height=3----
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)

## ----viewTracksHeight,fig.cap='plot data with adjusted track height',fig.width=8,fig.height=3----
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)

## ----viewTracksNames,fig.cap='change the track names',fig.width=8,fig.height=3----
names(trackList) <- c("cpsf160", "fox2", rep("HSPA8", 5))
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)

## ----viewTracksPaired,fig.cap='show two data in the same track',fig.width=8,fig.height=2.4----
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)

## ----viewTracksFlipped,fig.cap='show data in the flipped track',fig.width=8,fig.height=4----
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)

## ----themeBW,fig.cap='balck & white theme',fig.width=8,fig.height=4-----------
optSty <- optimizeStyle(trackList(repA, fox2, trs), theme="bw")
trackList <- optSty$tracks
viewerStyle <- optSty$style
vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)

## ----themeCol,fig.cap='colorful theme',fig.width=8,fig.height=4---------------
optSty <- optimizeStyle(trackList(repA, fox2, trs), theme="col")
trackList <- optSty$tracks
viewerStyle <- optSty$style
vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)

## ----themeSafe,fig.cap='safe theme',fig.width=8,fig.height=4------------------
optSty <- optimizeStyle(trackList(repA, fox2, trs), theme="safe")
trackList <- optSty$tracks
viewerStyle <- optSty$style
vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)

## ----axisBreak,fig.cap='axis with breaks',fig.width=8,fig.height=4------------
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)

## ----browseTrack,fig.cap='interactive tracks',fig.width=6,fig.height=4--------
browseTracks(trackList, gr=gr)

## ----plotgeneTrack,fig.cap='Plot multiple genes in one track',fig.width=6,fig.height=4----
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
grW <- parse2GRanges("chr11:122,830,799-123,116,707")
ids <- getGeneIDsFromTxDb(grW, TxDb.Hsapiens.UCSC.hg19.knownGene)
symbols <- mget(ids, org.Hs.egSYMBOL)
genes <- geneTrack(ids, TxDb.Hsapiens.UCSC.hg19.knownGene, 
                   symbols, asList=FALSE)
optSty <- optimizeStyle(trackList(repA, fox2, genes), theme="safe")
trackListW <- optSty$tracks
viewerStyleW <- optSty$style
viewTracks(trackListW, gr=grW, viewerStyle=viewerStyleW)

## ----viewTracksOperator1,fig.cap='show data with operator "+"',fig.width=8,fig.height=4----
newtrack <- repA
## Must keep the 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="+")

## ----viewTracksOperator2,fig.cap='show data with operator "-"',fig.width=8,fig.height=4----
viewTracks(trackList(newtrack, trs), gr=gr, viewerStyle=viewerStyle, operator="-")

## ----viewTracksOperator3,fig.cap='show data with operator "-"',fig.width=8,fig.height=4----
newtrack$dat <- GRoperator(newtrack$dat, newtrack$dat2, col="score", operator="-")
newtrack$dat2 <- GRanges()
viewTracks(trackList(newtrack, trs), gr=gr, viewerStyle=viewerStyle)

## ----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)))

## ----lolliplot2, fig.width=6, fig.height=3------------------------------------
features$fill <- c("#FF8833", "#51C6E6", "#DFA32D")
lolliplot(sample.gr, features)

## ----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)

## ----lolliplotOpacity, fig.width=6, fig.height=3------------------------------
sample.gr$alpha <- sample(100:255, length(SNP), replace = TRUE)/255
lolliplot(sample.gr, features)

## ----lolliplot.index, fig.width=6, fig.height=3-------------------------------
sample.gr$label <- as.character(1:length(sample.gr))
sample.gr$label.col <- ifelse(sample.gr$alpha>0.5, "white", "black")
lolliplot(sample.gr, features)

## ----lolliplot4, fig.width=6, fig.height=3------------------------------------
features$height <- c(0.02, 0.05, 0.08)
lolliplot(sample.gr, features)
## Specifying the height and its unit
features$height <- list(unit(1/16, "inches"), 
                     unit(3, "mm"), 
                     unit(12, "points"))
lolliplot(sample.gr, features)

## ----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 transcript
names(features.mul) <- features.mul$featureLayerID
lolliplot(sample.gr, features.mul)

## ----lolliplot4.1, fig.width=6, fig.height=3.5--------------------------------
#Note: the score value is an integer less than 10
sample.gr$score <- sample.int(5, length(sample.gr), replace = TRUE)
lolliplot(sample.gr, features)
##Remove y-axis
lolliplot(sample.gr, features, yaxis=FALSE)

## ----lolliplot4.2, fig.width=6, fig.height=4.5--------------------------------
#Try a score value greater than 10
sample.gr$score <- sample.int(20, length(sample.gr), replace=TRUE)
lolliplot(sample.gr, features)
#Try a float numeric score
sample.gr$score <- runif(length(sample.gr))*10
lolliplot(sample.gr, features)
# Score should not be smaller than 1
# remove the alpha for following samples
sample.gr$alpha <- NULL

## ----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)

## ----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] <- "y-axis" 
lolliplot(sample.gr, features, yaxis=yaxis)

## ----lolliplot.jitter, fig.width=6, fig.height=4.5----------------------------
sample.gr$dashline.col <- sample.gr$color
lolliplot(sample.gr, features, jitter="label")

## ----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, please try to set the legend by list.
## see more examples in the section [Plot multiple samples](#plot-multiple-samples)
legend <- list(legend)
lolliplot(sample.gr, features, legend=legend)

# from version 1.21.8, users can also try to set legend 
# as a column name in the metadata of GRanges.
sample.gr.newlegend <- sample.gr
sample.gr.newlegend$legend <- LETTERS[sample.gr$color]
lolliplot(sample.gr.newlegend, features, legend="legend")

## ----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)

## ----lolliplot.xlab.ylab.title, fig.width=6, fig.height=5.2-------------------
lolliplot(sample.gr.rot, features, legend=legend, ylab="y label here")
grid.text("label of x-axis here", x=.5, y=.01, just="bottom")
grid.text("title here", x=.5, y=.98, just="top", 
          gp=gpar(cex=1.5, fontface="bold"))

## ----lolliplot.labels.ctl.one.by.one, fig.width=6, fig.height=5---------------
label.parameter.gp.brown <- gpar(col="brown")
label.parameter.gp.blue <- gpar(col="blue")
label.parameter.gp.red <- gpar(col="red")
sample.gr$label.parameter.gp <- sample(list(label.parameter.gp.blue,
                                            label.parameter.gp.brown,
                                            label.parameter.gp.red),
                                       length(sample.gr), replace = TRUE)
lolliplot(sample.gr, features)

## ----lolliplotShape, fig.width=6, fig.height=4.5------------------------------
## shape must be "circle", "square", "diamond", "triangle_point_up", or "triangle_point_down"
available.shapes <- c("circle", "square", "diamond", 
                      "triangle_point_up", "triangle_point_down")
sample.gr$shape <- sample(available.shapes, size = length(sample.gr), replace = TRUE)
sample.gr$legend <- paste0("legend", as.numeric(factor(sample.gr$shape)))
lolliplot(sample.gr, features, type="circle", legend = "legend")

## ----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")

## ----lolliplotFlag, fig.width=6, fig.height=4---------------------------------
sample.gr.flag <- sample.gr
sample.gr.flag$label <- names(sample.gr) ## move the names to metadata:label
names(sample.gr.flag) <- NULL
lolliplot(sample.gr.flag, features, 
          ranges=GRanges("chr1", IRanges(0, 1600)), ## use ranges to leave more space on the right margin.
          type="flag")
## change the flag rotation angle
sample.gr.flag$label.rot <- 15
sample.gr.flag$label.rot[c(2, 5)] <- c(60, -15)
sample.gr.flag$label[7] <- "I have a long name"
lolliplot(sample.gr.flag, features, 
          ranges=GRanges("chr1", IRanges(0, 1600)), 
          type="flag")

## ----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", "alpha", "shape", "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 the color should be no less than that of value1 or value2
sample.gr$color <- rep(list(c("#87CEFA", '#98CE31')), length(SNP))
sample.gr$border <- "gray30"
lolliplot(sample.gr, features, type="pie")

## ----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)

## ----lolliplot.multiple.type, fig.width=6, fig.height=7.5---------------------
sample2.gr$score <- sample2.gr$value1 ## The circle layout needs the score column 
lolliplot(list(A=sample.gr, B=sample2.gr), 
          list(x=features, y=features2), 
          type=c("pie", "circle"), legend=legends)

## ----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")

## ----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]])

## ----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)

## ----ProteinsAPI, fig.width=6, fig.height=3-----------------------------------
library(httr) # load library to get data from REST API
APIurl <- "https://www.ebi.ac.uk/proteins/api/" # base URL of the API
taxid <- "9606" # human tax ID
gene <- "TP53" # target gene
orgDB <- "org.Hs.eg.db" # org database to get the uniprot accession id
eid <- mget("TP53", get(sub(".db", "SYMBOL2EG", orgDB)))[[1]]
chr <- mget(eid, get(sub(".db", "CHR", orgDB)))[[1]]
accession <- unlist(lapply(eid, function(.ele){
  mget(.ele, get(sub(".db", "UNIPROT", orgDB)))
}))
stopifnot(length(accession)<=20) # max number of accession is 20

tryCatch({ ## in case the internet connection does not work
  featureURL <- paste0(APIurl, 
                       "features?offset=0&size=-1&reviewed=true",
                       "&types=DNA_BIND%2CMOTIF%2CDOMAIN",
                       "&taxid=", taxid,
                       "&accession=", paste(accession, collapse = "%2C")
  )
  response <- GET(featureURL)
  if(!http_error(response)){
    content <- content(response)
    content <- content[[1]]
    acc <- content$accession
    sequence <- content$sequence
    gr <- GRanges(chr, IRanges(1, nchar(sequence)))
    domains <- do.call(rbind, content$features)
    domains <- GRanges(chr, IRanges(as.numeric(domains[, "begin"]),
                                     as.numeric(domains[, "end"]),
                                     names = domains[, "description"]))
    names(domains)[1] <- "DNA_BIND" ## this is hard coding.
    domains$fill <- 1+seq_along(domains)
    domains$height <- 0.04
    ## GET variations. This part can be replaced by user-defined data.
    variationURL <- paste0(APIurl,
                           "variation?offset=0&size=-1",
                           "&sourcetype=uniprot&dbtype=dbSNP",
                           "&taxid=", taxid,
                           "&accession=", acc)
    response <- GET(variationURL)
    if(!http_error(response)){
      content <- content(response)
      content <- content[[1]]
      keep <- sapply(content$features, function(.ele) length(.ele$evidences)>2 && # filter the data by at least 2 evidences
                       !grepl("Unclassified", .ele$clinicalSignificances)) # filter the data by classified clinical significances.
      nkeep <- c("wildType", "alternativeSequence", "begin", "end",
                 "somaticStatus", "consequenceType", "score")
      content$features <- lapply(content$features[keep], function(.ele){
        .ele$score <- length(.ele$evidences)
        unlist(.ele[nkeep]) 
      })
      variation <- do.call(rbind, content$features)
      variation <- 
        GRanges(chr, 
                IRanges(as.numeric(variation[, "begin"]),
                        width = 1,
                        names = paste0(variation[, "wildType"],
                                       variation[, "begin"],
                                       variation[, "alternativeSequence"])),
                score = as.numeric(variation[, "score"]),
                color = as.numeric(factor(variation[, "consequenceType"]))+1)
      variation$label.parameter.gp <- gpar(cex=.5)
      lolliplot(variation, domains, ranges = gr, ylab = "# evidences", yaxis = FALSE)
    }else{
      message("Can not get variations. http error")
    }
  }else{
    message("Can not get features. http error")
  }
},error=function(e){
  message(e)
},warning=function(w){
  message(w)
},interrupt=function(i){
  message(i)
})

## ----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"))
if(.Platform$OS.type!="windows"){# This line is for avoiding error from VariantAnnotation in the windows platform, which will be removed when VariantAnnotation's issue gets fixed.
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(mutation.frequency) <- "UCSC" 
}
seqlevelsStyle(gr) <- "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)
if(.Platform$OS.type!="windows"){
lolliplot(mutation.frequency, features, ranges=gr)
}

## ----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)

## ----methylation.hide, echo=FALSE---------------------------------------------
methy <- import(system.file("extdata", "methy.bed", package="trackViewer"), "BED")

## ----fig.width=6,fig.height=4-------------------------------------------------
lolliplot(methy, features, ranges=gr, type="pin")

## ----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")

## ----fig.width=6,fig.height=2.5-----------------------------------------------
methy$cex <- 1
lolliplot(methy, features, ranges=gr, rescale = TRUE)
## by set percentage for features and non-features segments
xaxis <- c(50968014, 50968514, 50968710, 50968838, 50970514)
rescale <- c(.3, .4, .3)
lolliplot(methy, features, ranges=gr, type="pin",
          rescale = rescale, xaxis = xaxis)
## by set data.frame to rescale
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)
)
lolliplot(methy, features, ranges=gr, type="pin",
          rescale = rescale, xaxis = xaxis)

## ----fig.width=6,fig.height=7.5-----------------------------------------------
grSplited <- tile(gr, n=2)
lolliplot(methy, features, ranges=grSplited, type="pin")

## ----fig.width=4.5,fig.height=3-----------------------------------------------
dandelion.plot(methy, features, ranges=gr, type="pin")

## ----fig.width=4.5,fig.height=3-----------------------------------------------
methy$color <- 3
methy$border <- "gray"
## Score info is required and 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")

## ----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)

## ----fig.width=4.5,fig.height=3.5---------------------------------------------
## Less dandelions
dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=1/10)

## ----fig.width=4.5,fig.height=4-----------------------------------------------
## More dandelions
dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=1/100)

## ----fig.width=4,fig.height=4-------------------------------------------------
maxgaps <- tile(gr, n = 10)[[1]]
dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=maxgaps)

## ----fig.width=4.5,fig.height=3-----------------------------------------------
dandelion.plot(methy, features, ranges=gr, type="pie", 
               maxgaps=1/100, yaxis = TRUE, heightMethod = mean,
               ylab='mean of methy scores')

## ----fig.width=4.5,fig.height=3-----------------------------------------------
yaxis = c(0, 0.5, 1)
dandelion.plot(methy, features, ranges=gr, type="pie", 
               maxgaps=1/100, yaxis = yaxis, heightMethod = mean,
               ylab='mean of methy scores')

## ----fig.width=8,fig.height=4-------------------------------------------------
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)

## ----eval=FALSE---------------------------------------------------------------
#  ideo <- loadIdeogram("hg38")

## ----echo=FALSE, results="hide"-----------------------------------------------
path <- system.file("extdata", "ideo.hg38.rds", package = "trackViewer")
ideo <- readRDS(path)

## -----------------------------------------------------------------------------
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")))

## -----------------------------------------------------------------------------
library(InteractionSet)
gi <- readRDS(system.file("extdata", "nij.chr6.51120000.53200000.gi.rds", package="trackViewer"))
head(gi)
## hicexplorer:hicConvertFormat tool can be used to convert other formats into GInteractions
## eg: hicConvertFormat -m mESC_rep.hic --inputFormat hic --outputFormat ginteractions -o mESC_rep.ginteractions --resolutions 5000
## please note that metadata:score is used for plot.
range <- GRanges("chr6", IRanges(51120000, 53200000))
tr <- gi2track(gi)
ctcf <- readRDS(system.file("extdata", "ctcf.sample.rds", package="trackViewer"))
viewTracks(trackList(ctcf, tr), gr=range, autoOptimizeStyle = TRUE)

## -----------------------------------------------------------------------------
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)

## ----sessionInfo, results='asis'----------------------------------------------
sessionInfo()