---
title: "trackViewer Vignette"
author: "Jianhong Ou, Lihua Julie Zhu"
date: "`r BiocStyle::doc_date()`"
package: "`r BiocStyle::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:
  html_document:
    theme: simplex
    toc: true
    toc_float: true
    toc_depth: 4
    fig_caption: true
---

```{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)
  library(httr)
})
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** can be used to plot coverage and annotation tracks. **TrackViewer** is a lightweight visualization tool for generating interactive figures for publication. Not only can **trackViewer** be used to visualize coverage and annotation tracks, but it can also be employed to generate lollipop/dandelion plots that depict dense methylation/mutation/variant data to facilitate an integrative analysis of these multi-omics data. It leverages **Gviz** and **rtracklayer**, is easy to use, and has a low memory and cpu consumption. In addition, we implemented a web application of **trackViewer** leveraging **Shiny** package. The web application of **trackViewer** is available at https://github.com/jianhong/trackViewer.documentation/tree/master/trackViewerShinyApp.

```{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**. 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)
```

 **trackViewer** not only has the functionalities to produce the figures generated by **Gviz**, as shown in the Figure above, but also provides additional plotting styles as shown in the Figure below. The mimimalist design requires minimum input from the users while retaining the flexibility to change the 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=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)
```


**Gviz** requires huge memory space to handle big wig files. To solve this problem, we rewrote the import function in **trackViewer** by importing the entire file first and parsing it later when plot. As a result, **trackViewer** decreases the import time from 180 min to 21 min and the memory cost from 10G to 5.32G for a half giga wig file (GSM917672). 

# Browse 
## Steps of using trackViewer \Biocpkg{trackViewer}
### Step 1. Import data
The function **importScore** is used to import BED, WIG, bedGraph or BigWig files. 
The function **importBam** is employed to import the bam files. Here is an 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 any strand info, 
## we need to set it manually.
strand(repA$dat) <- "-"
strand(repA$dat2) <- "+"
```

The function **coverageGR** could be used to calculate the coverage after the data is imported.

```{r 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)=="-"]
```

### Step 2. Build the gene model
The gene model can be built for a given genomic range using **geneModelFromTxdb** function which uses the **TranscriptDb** object as the 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)
```
Users can generate a track object with the geneTrack function by inputting a TxDb and a list of gene Entrez IDs. Entrez IDs can be obtained from other types of gene IDs such as gene symbol by using the ID mapping function. For example, to generate a track object given gene FMR1 and human TxDb, refer to the code below.

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

### Step 3. View the tracks
Use **viewTracks** function to plot data and annotation information along genomic coordinates.
**addGuideLine** or **addArrowMark** can be used to highlight a specific region.
    
```{r 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)
```


## Adjust the styles
### Adjust the x-axis or the X scale
In most cases, researchers are interested in the relative position of the peaks in the gene. 
Sometimes, margin needs to be adjusted to be able to show the entire  gene model. 
The Figure below shows how to add an X scale (x-scale) and 
remove the x-axis using the **setTrackXscaleParam** and **setTrackViewerStyleParam** functions. 

```{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=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)
```


### Adjust the y-axis
The y-axis can be put to the right side of the track by setting the main slot to FALSE in the y-axis slot of each track. 
In addition, the limit of y-axis (ylim) can be set by **setTrackStyleParam**. 
```{r 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)
```

### Adjust the label of y-axis
The style of y-axis can be changed by setting the ylabgp slot in the style of each track.
```{r 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)
```

The y-axis label can be put at the top or the bottom of each track.
```{r 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)
```

For each transcript, the transcript name can be put on the upstream or downstream of the transcript.

```{r 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)
```


### Adjust the track color
The track color can be changed by setting the color slot in the style of each track.
The first color is for the dat slot of 
**track** and the second color is for the dat2 slot.
```{r 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)
```

### Adjust the track height
The track height can be changed by setting the height slot in the 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=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)
```

### Change the track names
The track names such as gene model names can be edited easily by changing the names of **trackList**.

```{r 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)
```

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

### Flip the x-axis
The x-axis can be horizotally flipped for the genes in the negative strand.

```{r 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)
```

### Optimize the theme

Currently, we support two themes: bw (black and white) and col (colored).
```{r 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)
```

```{r 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)
```

```{r 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)
```

### Plot with breaks

We could plot the tracks with breaks by setting multiple genomic ranges.
```{r 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)
```

## Generate and edit an interactive plot using the browseTracks function
As shown above, figures produced by trackViewer are highly customizable, allowing users to alter the label, symbol, color, and size with various functions.

For users who prefer to modify the look and feel of a figure interactively, they can use the function `browseTracks` to draw interactive tracks, leveraging the htmlwidgets package. 

```{r browseTrack,fig.cap='interactive tracks',fig.width=6,fig.height=4}
browseTracks(trackList, gr=gr)
```
The videos at https://youtu.be/lSmeTu4WMlc and  https://youtu.be/lvF0tnJiHQI illustrate how to generate and modify an interactive plot. Please note that the interactive feature is only fully implemented in version 1.19.14 or later.  

## Plot multiple genes in one track

```{r 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)
```


# Operators
If you are interested in drawing a combined track from two input tracks, e.g, adding or substractiong one from the other, then you can try one of the operators such as + and - as showing below.

```{r 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="+")
```

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

Alternatively, you can try **GRoperator** before viewing tracks.

```{r 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)
```

# Lolliplot
Lolliplot is for the visualization of the methylation/variant/mutation data. 
```{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 the lolliplot color
### Change the color of the features.
```{r lolliplot2, fig.width=6, fig.height=3}
features$fill <- c("#FF8833", "#51C6E6", "#DFA32D")
lolliplot(sample.gr, features)
```

### Change the color of the 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)
```

### Change the opacity of the lollipop.
```{r lolliplotOpacity, fig.width=6, fig.height=3}
sample.gr$alpha <- sample(100:255, length(SNP), replace = TRUE)/255
lolliplot(sample.gr, features)
```


## Add the 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 <- ifelse(sample.gr$alpha>0.5, "white", "black")
lolliplot(sample.gr, features)
```

## Change the height of the features
```{r 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)
```

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

## Change the height of a lollipop plot
```{r 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)
```
```{r 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
```

## Customize the x-axis label position
```{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 the y-axis label position
```{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] <- "y-axis" 
lolliplot(sample.gr, features, yaxis=yaxis)
```

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

## Add a 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, 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")
```

## Control the labels

Users can control the paramters of labels by naming the metadata start as label.parameter  such as label.parameter.rot or 
label.parameter.gp. The parameter is used for [grid.text](https://www.rdocumentation.org/packages/grid/topics/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 the ylab, please try to set the labels in the ylab.
Please note that **lolliplot** does not support any parameters to set the title and xlab. If you want to 
add the title and xlab, please try to add them by [grid.text](https://www.rdocumentation.org/packages/grid/topics/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("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"))
```


Users can control the labels one by one by setting label.parameter.gp. 
Please note that for each label, the label.parameter.gp must be a list.

```{r 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)
```


## Change the lolliplot type

### Change the shape for "circle" plot

```{r 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")
```

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

### Flag

```{r 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")
```

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


## 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 possible.
```{r 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)
```

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

## EMBL-EBI Proteins API

Following code will show how to use [EBI Proteins REST API](https://www.ebi.ac.uk/proteins/api/doc/) to get annotations of protein domains.

```{r 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. error code: ", http_status(response))
    }
  }else{
    message("Can not get features. error code: ", http_status(response))
  }
},error=function(e){
  message(e)
})
```


## 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"))
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 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 the node size
In the above example, 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 scale of the x-axis (xscale)
In the above examples, some of the nodes are moved too far from the original coordinates.
To rescale, the x-axis could be reset as below.
```{r 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)
```

# Dandelion Plot
Sometimes, there are as many as hundreds of SNPs invoved in one gene. Dandelion plot can be used to depict such dense SNPs. Please note that the height of the dandelion indicates 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, i.e., type "fan". The area of the fan indicates 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 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")
```
```{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=3.5}
## Less dandelions
dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=1/10)
```
```{r fig.width=4.5,fig.height=4}
## More dandelions
dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=1/100)
```
Users can also specity the maximum distance between neighboring dandelions by settimg the maxgaps as a GRanges object.

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

## Add y-axis (yaxis)
Set yaxis to TRUE to add y-axis, and set heightMethod=mean to use the mean score as the height.
```{r 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')
```
```{r 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')
```

# Plot the lollipop plot with the coverage and annotation tracks
```{r 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)
```

# Ideogram Plot

Plot ideograms with a list of chromosomes and a genome.
```{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 genomic interactions data as tracks.

```{r}
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)
```


Different from most of the available tools, plotGInteractions try to plot the data with the 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 the enhancers, promoters or genes. 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)
```

# Web application of **trackViewer**
We created a web application of **trackViewer** (available in 1.19.14 or later) by leveraging the R package **Shiny**. The web application of trackViewer and sample data are available at  https://github.com/jianhong/trackViewer.documentation/tree/master/trackViewerShinyApp. Here is a demo on how to to use the web application at https://www.nature.com/articles/s41592-019-0430-y#Sec2 Supplementary Video 5.

# Session Info
```{r sessionInfo, results='asis'}
sessionInfo()
```