## ---- echo=FALSE, results="hide", warning=FALSE-------------------------------
suppressPackageStartupMessages({
library(trackViewer)
library(rtracklayer)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
library(VariantAnnotation)
library(httr)
})
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
## ----lolliplot1, fig.width=6, fig.height=3------------------------------------
library(trackViewer)
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)
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$node.label <- as.character(1:length(sample.gr))
sample.gr$node.label.col <-
ifelse(sample.gr$alpha>0.5 | sample.gr$color==1, "white", "black")
sample.gr$node.label.cex <- sample.int(3, length(sample.gr), replace = TRUE)/2
lolliplot(sample.gr, features)
sample.gr$node.label.cex <- 1 ## change it back for pretty showcase.
## ----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) ## define the position
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.col <- "brown"
lolliplot(sample.gr.rot, features, legend=legend)
## try different colors
sample.gr.rot$label.parameter.col <- sample.int(7,
length(sample.gr),
replace = TRUE)
sample.gr.rot$label.parameter.draw <- TRUE
sample.gr.rot$label.parameter.draw[[1]] <- FALSE ## remove the first label
lolliplot(sample.gr.rot, features, legend=legend)
## ----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)
## ----lolliplot.labels.on.feature, fig.width=6, fig.height=5-------------------
lolliplot(sample.gr, features, label_on_feature=TRUE)
## ----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.motiflogo, fig.width=6, fig.height=5.2-------------------------
library(motifStack)
pcms<-readPCM(file.path(find.package("motifStack"), "extdata"),"pcm$")
sample.gr.rot$label.parameter.pfm <- pcms[sample(seq_along(pcms),
length(sample.gr.rot),
replace = TRUE)]
lolliplot(sample.gr.rot, features, legend=legend)
## ----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$node.label.rot <- 15
sample.gr.flag$node.label.rot[c(2, 5)] <- c(60, -15)
sample.gr.flag$label[7] <- "I have a long name"
sample.gr.flag$node.label.cex <- 1
sample.gr.flag$node.label.cex[4] <- 2
lolliplot(sample.gr.flag, features,
ranges=GRanges("chr1", IRanges(0, 1600)),## use ranges to leave more space on the right margin.
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$node.label.col <- NULL
x <- sample.int(100, length(SNP))
sample.gr$value1 <- x
sample.gr$value2 <- 100 - x # for pie plot, 2 value columns are required.
## 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))
## remove one mutation label
sample.gr.mul.patient$label.parameter.draw <- TRUE
sample.gr.mul.patient$label.parameter.draw[
names(sample.gr.mul.patient)==
sample(unique(names(sample.gr.mul.patient)), 1)] <-
FALSE
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, eval=FALSE, echo=TRUE------------
# 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"
if(!grepl("chr", seqlevels(mutation.frequency)[1])){
seqlevels(mutation.frequency) <-
paste0("chr", seqlevels(mutation.frequency))
}
}
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,
# table="wgEncodeHaibMethylRrbs",
# 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=2.5-----------------------------------------------
lolliplot(methy, features, ranges=gr, rescale = "exon")
# exon region occupy 99% of the plot region.
lolliplot(methy, features, ranges=gr, rescale = "exon_99")
lolliplot(methy, features, ranges=gr, rescale = "intron")
## ----fig.width=9,fig.height=8-------------------------------------------------
grSplited <- tile(gr, n=2)
lolliplot(methy, features, ranges=grSplited, type="pin")
## ----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
extdata <- system.file("extdata", package="trackViewer",
mustWork=TRUE)
repA <- importScore(file.path(extdata, "cpsf160.repA_-.wig"),
file.path(extdata, "cpsf160.repA_+.wig"),
format="WIG")
fox2 <- importScore(file.path(extdata, "fox2.bed"), format="BED",
ranges=GRanges("chr11", IRanges(122830799, 123116707)))
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)
## ----sessionInfo, results='asis'----------------------------------------------
sessionInfo()