trackViewer 氨基酸位点变异位置图谱展示

http://www.bioconductor.org/packages/release/bioc/vignettes/trackViewer/inst/doc/trackViewer.html

 

Contents

1Introduction

There are two packages available in Bioconductor for visualizing genomic data: rtracklayerand Gviz. rtracklayer provides an interface to genome browsers and associated annotation tracks. Gviz plots data and annotation information along genomic coordinates. TrackViewer is a light-weighted visualization tool for generating neat figures for publication. It utilizes Gviz, is easy to use, and has a low memory and cpu consumption.

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="-")
Plot data with **Gviz** and **trackViewer**. Note that **trackViewer** can generate similar figure as **Gviz** with several lines of simple codes.

Figure 1: Plot data with Gviz and trackViewer
Note that trackViewer can generate similar figure as Gviz with several lines of simple codes.

TrackViewer not only has the functionalities to plot the figures generated by Gviz, as shown in Figure above, but also provides additional plotting styles as shown in Figure below. The mimimalist design requires minimum input from users while retaining the flexibility to change output style easily.

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

Figure 2: 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.

It requires huge memory space to handle big wig files. To solve this problem, trackViewerrewrote the import function to import whole file first and parse it later when plot. trackViewer provides higher import speed (21 min vs. over 180 min) and acceptable memory cost (5.32G vs. over 10G) for a half giga wig file (GSM917672) comparing to Gviz.

2Browse tracks

2.1Steps of using

2.1.1step1 import data

Function importScore is used to import BED, WIG, bedGraph or BigWig files. Function importBam is employed to import bam file. Here is the example.

library(trackViewer)
extdata <- system.file("extdata", package="trackViewer",
                       mustWork=TRUE)
repA <- importScore(file.path(extdata, "cpsf160.repA_-.wig"),
                    file.path(extdata, "cpsf160.repA_+.wig"),
                    format="WIG")
## because the wig file does not contain strand info, 
## we need to set it manually
strand(repA$dat) <- "-"
strand(repA$dat2) <- "+"

Function coverageGR could be used to calculate coverage after importing if needed.

fox2 <- importScore(file.path(extdata, "fox2.bed"), format="BED",
                    ranges=GRanges("chr11", IRanges(122929000, 122931000)))
dat <- coverageGR(fox2$dat)
## we can split the data by strand into two different track channels
## here we set the dat2 slot to save the negative strand info, 
## reverse order as previous. 
fox2$dat <- dat[strand(dat)=="+"]
fox2$dat2 <- dat[strand(dat)=="-"]

2.1.2step2 build gene model

The gene model can be built for a given genomic range using geneModelFromTxdbfunction which uses TranscriptDb object as input.

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)

2.1.3step3 view the tracks

Use viewTracks function to plot data and annotation information along genomic coordinates. addGuideLine or addArrowMark can be used to highlight the peaks.

viewerStyle <- trackViewerStyle()
setTrackViewerStyleParam(viewerStyle, "margin", c(.1, .05, .02, .02))
vp <- viewTracks(trackList(repA, fox2, trs), 
                 gr=gr, viewerStyle=viewerStyle, 
                 autoOptimizeStyle=TRUE)
addGuideLine(c(122929767, 122929969), vp=vp)
addArrowMark(list(x=122929650, 
                  y=2), # 2 means track 2 from bottom.
             label="label",
             col="blue",
             vp=vp)
plot data and annotation information along genomic coordinates

Figure 3: plot data and annotation information along genomic coordinates

2.2Adjust the styles

2.2.1adjust x axis or x-scale

In most cases, researchers are interested in the relative position of peaks in the gene. Sometimes, margin needs to be adjusted to be able to show the entire gene model. Figure below shows how to add an x-scale and remove x-axis using addGuideLine Function .

optSty <- optimizeStyle(trackList(repA, fox2, trs))
trackList <- optSty$tracks
viewerStyle <- optSty$style
setTrackViewerStyleParam(viewerStyle, "xaxis", FALSE)
setTrackViewerStyleParam(viewerStyle, "margin", c(.01, .05, .01, .01))
setTrackXscaleParam(trackList[[1]], "draw", TRUE)
setTrackXscaleParam(trackList[[1]], "gp", list(cex=.5))
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)
plot data with x-scale

Figure 4: plot data with x-scale

2.2.2adjust y axis

y-axis can be put to right side of the track by setting main slot to FALSE in y-axis slot of each track. And ylim can be set by setTrackStyleParam.

setTrackViewerStyleParam(viewerStyle, "margin", c(.01, .05, .01, .05))
for(i in 1:2){
    setTrackYaxisParam(trackList[[i]], "main", FALSE)
}
## adjust y scale
setTrackStyleParam(trackList[[1]], "ylim", c(0, 25))
setTrackStyleParam(trackList[[2]], "ylim", c(-25, 0))
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)
plot data with y-axis in right side

Figure 5: plot data with y-axis in right side

2.2.3adjust y label

Y label style can be changed by setting the ylabgp slot in style of each track.

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)
plot data with adjusted color and size of y label

Figure 6: plot data with adjusted color and size of y label

Y label can be also put to top or bottom of each track.

setTrackStyleParam(trackList[[1]], "ylabpos", "bottomleft")
setTrackStyleParam(trackList[[2]], "ylabpos", "topright")
setTrackStyleParam(trackList[[2]], "marginTop", .2)
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)
plot data with adjusted y label position

Figure 7: plot data with adjusted y label position

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

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)
plot data with adjusted transcripts name position

Figure 8: plot data with adjusted transcripts name position

2.2.4adjust track color

The track color can be changed by setting the color slot in style of each track. The first color is for dat slot of track and seconde color is for dat2 slot.

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)
plot data with adjusted track color

Figure 9: plot data with adjusted track color

2.2.5adjust track height

The track height can be changed by setting the height slot in style of each track. However, the total height for all the tracks should be 1.

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)
plot data with adjusted track height

Figure 10: plot data with adjusted track height

2.2.6change track names

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

names(trackList) <- c("cpsf160", "fox2", rep("HSPA8", 5))
viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)
change the track names

Figure 11: change the track names

2.2.7show paired data in the same track

trackViewer can be used to show to-be-compared data in the same track side by side.

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)
show two data in the same track

Figure 12: show two data in the same track

2.2.8flip the x-axis

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

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)
show data in the flipped track

Figure 13: show data in the flipped track

2.2.9optimize with theme

We support two themes now: bw and col.

optSty <- optimizeStyle(trackList(repA, fox2, trs), theme="bw")
trackList <- optSty$tracks
viewerStyle <- optSty$style
vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)
theme_bw

(#fig:theme_bw)theme_bw

optSty <- optimizeStyle(trackList(repA, fox2, trs), theme="col")
trackList <- optSty$tracks
viewerStyle <- optSty$style
vp <- viewTracks(trackList, gr=gr, viewerStyle=viewerStyle)
theme_col

(#fig:theme_col)theme_col

2.2.10plot with breaks

We could plot the tracks with breaks by set multiple genomic ranges.

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)
axis.break

(#fig:axis.break)axis.break

2.3browse tracks by htmlwidgets

To change the parameter may be complicated for users. Users may want to have a try function browseTracks.

browseTracks(trackList, gr=gr)
0.05.5repA-220fox2uc001pyo.3uc001pyp.3uc009zbc.3uc009zbd.2uc010rzu.2122,929,400122,929,600122,929,800122,930,000

Figure 14: browseTrack

3Operators

If there are two tracks and we want to draw the two track by adding or substract one from another, we can try operators.

newtrack <- repA
## must keep same format for dat and dat2
newtrack <- parseWIG(newtrack, "chr11", 122929275, 122930122)
newtrack$dat2 <- newtrack$dat
newtrack$dat <- fox2$dat2
setTrackStyleParam(newtrack, "color", c("blue", "red"))
viewTracks(trackList(newtrack, trs), 
           gr=gr, viewerStyle=viewerStyle, operator="+")
show data with operator "+"

Figure 15: show data with operator “+”

viewTracks(trackList(newtrack, trs), gr=gr, viewerStyle=viewerStyle, operator="-")
show data with operator "-"

Figure 16: show data with operator “-”

Or try GRoperator before view tracks.

newtrack$dat <- GRoperator(newtrack$dat, newtrack$dat2, col="score", operator="-")
newtrack$dat2 <- GRanges()
viewTracks(trackList(newtrack, trs), gr=gr, viewerStyle=viewerStyle)
show data with operator "-"

Figure 17: show data with operator “-”

4Lolliplot

Lolliplot is a mutation distribution graphics tool.

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

4.1Change lolliplot colors

4.1.1Change the color of features.

features$fill <- c("#FF8833", "#51C6E6", "#DFA32D")
lolliplot(sample.gr, features)

4.1.2Change the color of lollipop.

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)

4.2Add index labels in the node

sample.gr$label <- as.character(1:length(sample.gr))
sample.gr$label.col <- "white"
lolliplot(sample.gr, features)

4.3Change the height of features

features$height <- c(0.02, 0.05, 0.08)
lolliplot(sample.gr, features)

## keep the height by giving the unit
features$height <- list(unit(1/16, "inches"), 
                     unit(3, "mm"), 
                     unit(12, "points"))
lolliplot(sample.gr, features)

4.4Multiple transcripts in the features

The metadata ‘featureLayerID’ are used for drawing features in different layers.

features.mul <- rep(features, 2)
features.mul$height[4:6] <- list(unit(1/8, "inches"),
                                 unit(0.5, "lines"),
                                 unit(.2, "char"))
features.mul$fill <- c("#FF8833", "#F9712A", "#DFA32D", 
                       "#51C6E6", "#009DDA", "#4B9CDF")
end(features.mul)[5] <- end(features.mul[5])+50
features.mul$featureLayerID <- 
    paste("tx", rep(1:2, each=length(features)), sep="_")
names(features.mul) <- 
    paste(features.mul$featureLayerID, 
          rep(1:length(features), 2), sep="_")
lolliplot(sample.gr, features.mul)

## one name per transcripts
names(features.mul) <- features.mul$featureLayerID
lolliplot(sample.gr, features.mul)

4.5Change the height of lolliplot

#Note: the score value is integer less than 10
sample.gr$score <- sample.int(5, length(sample.gr), replace = TRUE)
lolliplot(sample.gr, features)

##remove yaxis
lolliplot(sample.gr, features, yaxis=FALSE)

#Try score value greater than 10
sample.gr$score <- sample.int(20, length(sample.gr), replace=TRUE)
lolliplot(sample.gr, features)

#Try float numeric score
sample.gr$score <- runif(length(sample.gr))*10
lolliplot(sample.gr, features)

# score should not be smaller than 1

4.6Customize xaxis label positions

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)

4.7Customize yaxis label positions

yaxis <- c(0, 5) ## define the position
lolliplot(sample.gr, features, yaxis=yaxis)

yaxis <- c(0, 5, 10, 15)
names(yaxis) <- yaxis # define the labels
names(yaxis)[3] <- "yaxis" 
lolliplot(sample.gr, features, yaxis=yaxis)

4.8Jitter the label only

sample.gr$dashline.col <- sample.gr$color
lolliplot(sample.gr, features, jitter="label")

4.9Add legend

legend <- 1:6 ## legend fill color
names(legend) <- paste0("legend", letters[1:6]) ## legend labels
lolliplot(sample.gr, features, legend=legend)

## use list to define more attributes. see ?grid::gpar to get more details.
legend <- list(labels=paste0("legend", LETTERS[1:6]), 
               col=palette()[6:1], 
               fill=palette()[legend])
lolliplot(sample.gr, features, legend=legend)

## if you have multiple tracks, try to set the legend by list.
## see more in section [Plot multiple samples](#plot-multiple-samples)
legend <- list(legend)
lolliplot(sample.gr, features, legend=legend)

4.10Control the labels

Use can control the paramters of labels by name the metadata start as label.parameter.such as label.parameter.rot, label.parameter.gp. The parameter is used for .

sample.gr.rot <- sample.gr
sample.gr.rot$label.parameter.rot <- 45
lolliplot(sample.gr.rot, features, legend=legend)

sample.gr.rot$label.parameter.rot <- 60
sample.gr.rot$label.parameter.gp <- gpar(col="brown")
lolliplot(sample.gr.rot, features, legend=legend)

If you want to change the text in ylab, please try to set the labels in ylab. lolliplot does not support the parameters for title and xlab. If you want to add title and xlab, please try to add them by .

lolliplot(sample.gr.rot, features, legend=legend, ylab="y label here")
grid.text("x label here", x=.5, y=.01, just="bottom")
grid.text("title here", x=.5, y=.98, just="top", 
          gp=gpar(cex=1.5, fontface="bold"))

4.11Change the type of lolliplot

4.11.1Google pin

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

4.11.2Pie plot

sample.gr$score <- NULL ## must be removed, because pie will consider all the numeric columns except column "color", "fill", "lwd", "id" and "id.col".
sample.gr$label <- NULL
sample.gr$label.col <- NULL
x <- sample.int(100, length(SNP))
sample.gr$value1 <- x
sample.gr$value2 <- 100 - x
## the length of color should be no less than the values number
sample.gr$color <- rep(list(c("#87CEFA", '#98CE31')), length(SNP))
sample.gr$border <- "gray30"
lolliplot(sample.gr, features, type="pie")

4.12Plot multiple samples

4.12.1multiple layers

SNP2 <- sample(4000:8000, 30)
x2 <- sample.int(100, length(SNP2), replace=TRUE)
sample2.gr <- GRanges("chr3", IRanges(SNP2, width=1, names=paste0("snp", SNP2)), 
             value1=x2, value2=100-x2)
sample2.gr$color <- rep(list(c('#DB7575', '#FFD700')), length(SNP2))
sample2.gr$border <- "gray30"

features2 <- GRanges("chr3", IRanges(c(5001, 5801, 7001), 
                                    width=c(500, 500, 405),
                                    names=paste0("block", 4:6)),
                    fill=c("orange", "gray30", "lightblue"),
                    height=unit(c(0.5, 0.3, 0.8), "cm"))
legends <- list(list(labels=c("WT", "MUT"), fill=c("#87CEFA", '#98CE31')), 
                list(labels=c("WT", "MUT"), fill=c('#DB7575', '#FFD700')))
lolliplot(list(A=sample.gr, B=sample2.gr), 
          list(x=features, y=features2), 
          type="pie", legend=legends)

Different layouts are also be possible.

sample2.gr$score <- sample2.gr$value1 ## circle layout need score column 
lolliplot(list(A=sample.gr, B=sample2.gr), 
          list(x=features, y=features2), 
          type=c("pie", "circle"), legend=legends)

4.12.2pie.stack layout

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

4.12.3caterpillar layout

Metadata SNPsideID is used to trigger caterpillar layout. SNPsideID must be ‘top’ or ‘bottom’.

sample.gr$SNPsideID <- sample(c("top", "bottom"), 
                               length(sample.gr),
                               replace=TRUE)
lolliplot(sample.gr, features, type="pie", 
          legend=legends[[1]])

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

4.13Variant Call Format (VCF) data

library(VariantAnnotation)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
gr <- GRanges("22", IRanges(50968014, 50970514, names="TYMP"))
tab <- TabixFile(fl)
vcf <- readVcf(fl, "hg19", param=gr)
mutation.frequency <- rowRanges(vcf)
mcols(mutation.frequency) <- cbind(mcols(mutation.frequency), 
                                   VariantAnnotation::info(vcf))
mutation.frequency$border <- "gray30"
mutation.frequency$color <- 
    ifelse(grepl("^rs", names(mutation.frequency)), "lightcyan", "lavender")
## plot Global Allele Frequency based on AC/AN
mutation.frequency$score <- mutation.frequency$AF*100
seqlevelsStyle(gr) <- seqlevelsStyle(mutation.frequency) <- "UCSC" 
trs <- geneModelFromTxdb(TxDb.Hsapiens.UCSC.hg19.knownGene,
                         org.Hs.eg.db,
                         gr=gr)
features <- c(range(trs[[1]]$dat), range(trs[[5]]$dat))
names(features) <- c(trs[[1]]$name, trs[[5]]$name)
features$fill <- c("lightblue", "mistyrose")
features$height <- c(.02, .04)
lolliplot(mutation.frequency, features, ranges=gr)

4.14Methylation data

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)
lolliplot(methy, features, ranges=gr, type="pin")

4.15Change node size

In above sample, some of the nodes overlap each other. To change the node size, cex prameter could be used.

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

4.16Change the xscale

In above samples, some of the nodes are moved too far from the original coordinates. To rescale the xaxis could be used.

methy$cex <- 1
lolliplot(methy, features, ranges=gr, rescale = TRUE)

rescale <- data.frame(
  from.start = c(50968014, 50968515, 50968838), 
  from.end   = c(50968514, 50968837, 50970514),
  to.start   = c(50968014, 50968838, 50969501), 
  to.end     = c(50968837, 50969500, 50970514)
)
xaxis <- c(50968014, 50968514, 50968710, 50968838, 50970514)
lolliplot(methy, features, ranges=gr, type="pin",
          rescale = rescale, xaxis = xaxis)

5Dandelion Plot

Sometimes, there are too many SNPs invoved in one gene. If you want to plot such as more than handred SNPs, Dandelion plot may be the good chioce. Please note, the height of the dandelion means the desity of the SNPs.

dandelion.plot(methy, features, ranges=gr, type="pin")

5.1change the type of Dandelion plot

There are one more type for dandelion plot, type “fan”. The area of the fan indicate the percentage of methylation or rate of mutation.

methy$color <- 3
methy$border <- "gray"
## score info is required, the score must be a number in [0, 1]
m <- max(methy$score)
methy$score <- methy$score/m
dandelion.plot(methy, features, ranges=gr, type="fan")

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)

5.2change the number of dandelions

## want less
dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=1/10)

## want more
dandelion.plot(methy, features, ranges=gr, type="circle", maxgaps=1/100)

6Plot lollipop plot with tracks

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
SNPs$cex <- 1
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)

7Ideogram Plot

Plot ideograms with data.

ideo <- loadIdeogram("hg38")
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")))

8Plot genomic interactions data

Plot GInteractions. Different from most of the available tools, plotGInteractions try to plot the data with 2D structure. The nodes indicate the region with interactions and the edges indicates the interactions. The size of the nodes are relative to the width of the region. The features could be enhancer, promoter or gene. The enhancer and promoter are shown as points with symbol 11 and 13.

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)

posted @ 2018-07-12 08:10  nkwy2012  阅读(3138)  评论(0)  编辑  收藏  举报